C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
srmatrix.hpp
00001 /*
00002 **  CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4)
00003 **
00004 **  Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik,
00005 **                          Universitaet Karlsruhe, Germany
00006 **            (C) 2000-2014 Wiss. Rechnen/Softwaretechnologie
00007 **                          Universitaet Wuppertal, Germany   
00008 **
00009 **  This library is free software; you can redistribute it and/or
00010 **  modify it under the terms of the GNU Library General Public
00011 **  License as published by the Free Software Foundation; either
00012 **  version 2 of the License, or (at your option) any later version.
00013 **
00014 **  This library is distributed in the hope that it will be useful,
00015 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
00016 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017 **  Library General Public License for more details.
00018 **
00019 **  You should have received a copy of the GNU Library General Public
00020 **  License along with this library; if not, write to the Free
00021 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 */
00023 
00024 /* CVS $Id: srmatrix.hpp,v 1.21 2014/01/30 17:23:49 cxsc Exp $ */
00025 
00026 #ifndef _CXSC_SRMATRIX_HPP_INCLUDED
00027 #define _CXSC_SRMATRIX_HPP_INCLUDED
00028 
00029 #include <real.hpp>
00030 #include <rmatrix.hpp>
00031 #include <srvector.hpp>
00032 #include <cidot.hpp>
00033 #include <vector>
00034 #include <algorithm>
00035 #include <iostream>
00036 #include <sparsedot.hpp>
00037 #include <sparsematrix.hpp>
00038 
00039 namespace cxsc {
00040 
00042 enum STORAGE_TYPE{triplet,compressed_row,compressed_column};
00043 
00044 class srmatrix_slice;
00045 class srmatrix_subv;
00046 class scmatrix;
00047 class scmatrix_slice;
00048 class scmatrix_subv;
00049 class simatrix;
00050 class simatrix_slice;
00051 class simatrix_subv;
00052 class sivector_slice;
00053 class scimatrix;
00054 class scimatrix_slice;
00055 class scimatrix_subv;
00056 class scivector_slice;
00057 
00058 
00059 inline bool comp_pair(std::pair<int,real> p1, std::pair<int,real> p2) {
00060   return p1.first < p2.first;
00061 }
00062 
00063 
00065 
00077 class srmatrix {
00078 
00079   private:
00080     std::vector<int> p;
00081     std::vector<int> ind;
00082     std::vector<real> x;
00083     int m;
00084     int n;
00085     int lb1,ub1,lb2,ub2;
00086 
00087   public:
00088 
00090     std::vector<int>& column_pointers() {
00091       return p;
00092     }
00093 
00095     std::vector<int>& row_indices() {
00096       return ind;
00097     }
00098 
00100     std::vector<real>& values() {
00101       return x;
00102     }
00103 
00105     const std::vector<int>& column_pointers() const {
00106       return p;
00107     }
00108 
00110     const std::vector<int>& row_indices() const {
00111       return ind;
00112     }
00113 
00115     const std::vector<real>& values() const {
00116       return x;
00117     }
00118 
00120     srmatrix() {
00121       p.push_back(0);
00122       m = n = 0;
00123       lb1 = lb2 = ub1 = ub2 = 0;
00124     }
00125 
00127     srmatrix(const int r, const int c) : m(r),n(c),lb1(1),ub1(r),lb2(1),ub2(c) {
00128       p = std::vector<int>((n>0) ? n+1 : 1, 0);
00129       ind.reserve(2*(m+n));
00130       x.reserve(2*(m+n));
00131 
00132       p[0] = 0;
00133     }
00134 
00136     srmatrix(const int r, const int c, const int e) : m(r),n(c),lb1(1),ub1(r),lb2(1),ub2(c) {
00137       p = std::vector<int>((n>0) ? n+1 : 1, 0);
00138       ind.reserve(e);
00139       x.reserve(e);
00140 
00141       p[0] = 0;
00142     }
00143 
00145 
00151     srmatrix(const int m, const int n, const int nnz, const intvector& rows, const intvector& cols, const rvector& values, const enum STORAGE_TYPE t = triplet) {
00152       if(t == triplet) {
00153          this->m = m;
00154          this->n = n;
00155          p = std::vector<int>(n+1,0);
00156          ind.reserve(nnz);
00157          x.reserve(nnz);
00158          lb1 = lb2 = 1;
00159          ub1 = m; ub2 = n;
00160 
00161          std::vector<triplet_store<real> > work;
00162          work.reserve(nnz);
00163 
00164          for(int k=0 ; k<nnz ; k++) {
00165            work.push_back(triplet_store<real>(rows[Lb(rows)+k],cols[Lb(cols)+k],values[Lb(values)+k]));
00166          }
00167 
00168          sort(work.begin(), work.end());
00169 
00170          int i=0;
00171 
00172          for(int j=0 ; j<n ; j++) {        
00173 
00174            while((unsigned int)i < work.size() && work[i].col == j ) {
00175                ind.push_back(work[i].row);
00176                x.push_back(work[i].val);
00177                i++;
00178            }
00179 
00180            p[j+1] = i;
00181          }
00182          
00183       } else if(t == compressed_row) {
00184 
00185          this->m = m;
00186          this->n = n;
00187          p = std::vector<int>(n+1,0);
00188          ind.reserve(nnz);
00189          x.reserve(nnz);
00190          lb1 = lb2 = 1;
00191          ub1 = m; ub2 = n;
00192 
00193          for(int i=0 ; i<n+1 ; i++)
00194            p[i] = rows[Lb(rows)+i];
00195 
00196          std::vector<triplet_store<real> > work;
00197          work.reserve(nnz);
00198 
00199          for(int j=0 ; j<n ; j++) {
00200            for(int k=p[j] ; k<p[j+1] ; k++) {
00201              work.push_back(triplet_store<real>(j,cols[Lb(cols)+k],values[Lb(values)+k]));
00202            }
00203          }
00204 
00205          sort(work.begin(), work.end());
00206 
00207          int i=0;
00208 
00209          for(int j=0 ; j<n ; j++) {        
00210 
00211            while((unsigned int)i < work.size() && work[i].col == j ) {
00212                ind.push_back(work[i].row);
00213                x.push_back(work[i].val);
00214                i++;
00215            }
00216 
00217            p[j+1] = i;
00218          }
00219     
00220       } else if(t == compressed_column) {
00221          this->m = m;
00222          this->n = n;
00223          p = std::vector<int>(n+1,0);
00224          ind.reserve(nnz);
00225          x.reserve(nnz);
00226          lb1 = lb2 = 1;
00227          ub1 = m; ub2 = n;
00228 
00229          for(int i=0 ; i<n+1 ; i++)
00230            p[i] = rows[Lb(rows)+i];
00231 
00232          std::vector<std::pair<int,real> > work;
00233          work.reserve(n);
00234 
00235          for(int j=0 ; j<n ; j++) {
00236            work.clear();
00237 
00238            for(int k=p[j] ; k<p[j+1] ; k++) {
00239              work.push_back(std::make_pair(cols[Lb(cols)+k],values[Lb(values)+k]));
00240            }
00241 
00242            std::sort(work.begin(),work.end(),comp_pair);
00243 
00244            for(unsigned int i=0 ; i<work.size() ; i++) {
00245              ind.push_back(work[i].first);
00246              x.push_back(work[i].second);
00247            }
00248          }
00249 
00250       }
00251 
00252     }
00253 
00254     
00256 
00263     srmatrix(const int m, const int n, const int nnz, const int* rows, const int* cols, const real* values, const enum STORAGE_TYPE t = triplet) {
00264       if(t == triplet) {
00265          this->m = m;
00266          this->n = n;
00267          p = std::vector<int>(n+1,0);
00268          ind.reserve(nnz);
00269          x.reserve(nnz);
00270          lb1 = lb2 = 1;
00271          ub1 = m; ub2 = n;
00272 
00273          std::vector<triplet_store<real> > work;
00274          work.reserve(nnz);
00275 
00276          for(int k=0 ; k<nnz ; k++) {
00277            work.push_back(triplet_store<real>(rows[k],cols[k],values[k]));
00278          }
00279 
00280          sort(work.begin(), work.end());
00281 
00282          int i=0;
00283 
00284          for(int j=0 ; j<n ; j++) {        
00285 
00286            while((unsigned int)i < work.size() && work[i].col == j ) {
00287                ind.push_back(work[i].row);
00288                x.push_back(work[i].val);
00289                i++;
00290            }
00291 
00292            p[j+1] = i;
00293          }
00294          
00295       } else if(t == compressed_row) {
00296 
00297          this->m = m;
00298          this->n = n;
00299          p = std::vector<int>(n+1,0);
00300          ind.reserve(nnz);
00301          x.reserve(nnz);
00302          lb1 = lb2 = 1;
00303          ub1 = m; ub2 = n;
00304 
00305          for(int i=0 ; i<n+1 ; i++)
00306            p[i] = rows[i];
00307 
00308          std::vector<triplet_store<real> > work;
00309          work.reserve(nnz);
00310 
00311          for(int j=0 ; j<n ; j++) {
00312            for(int k=p[j] ; k<p[j+1] ; k++) {
00313              work.push_back(triplet_store<real>(j,cols[k],values[k]));
00314            }
00315          }
00316 
00317          sort(work.begin(), work.end());
00318 
00319          int i=0;
00320 
00321          for(int j=0 ; j<n ; j++) {        
00322 
00323            while((unsigned int)i < work.size() && work[i].col == j ) {
00324                ind.push_back(work[i].row);
00325                x.push_back(work[i].val);
00326                i++;
00327            }
00328 
00329            p[j+1] = i;
00330          }
00331     
00332       } else if(t == compressed_column) {
00333          this->m = m;
00334          this->n = n;
00335          p = std::vector<int>(n+1,0);
00336          ind.reserve(nnz);
00337          x.reserve(nnz);
00338          lb1 = lb2 = 1;
00339          ub1 = m; ub2 = n;
00340 
00341          for(int i=0 ; i<n+1 ; i++)
00342            p[i] = rows[i];
00343 
00344          std::vector<std::pair<int,real> > work;
00345          work.reserve(n);
00346 
00347          for(int j=0 ; j<n ; j++) {
00348            work.clear();
00349 
00350            for(int k=p[j] ; k<p[j+1] ; k++) {
00351              work.push_back(std::make_pair(cols[k],values[k]));
00352            }
00353 
00354            std::sort(work.begin(),work.end(),comp_pair);
00355 
00356            for(unsigned int i=0 ; i<work.size() ; i++) {
00357              ind.push_back(work[i].first);
00358              x.push_back(work[i].second);
00359            }
00360          }
00361 
00362       }
00363 
00364     }
00365 
00367     srmatrix(const rmatrix& A) : m(ColLen(A)),n(RowLen(A)),lb1(Lb(A,1)),ub1(Ub(A,1)),lb2(Lb(A,2)),ub2(Ub(A,2)) {
00368       p = std::vector<int>((n>0) ? n+1 : 1, 0);
00369       ind.reserve((m*n*0.1 < 2*m) ? (int)(m*n*0.1) : 2*m);
00370       x.reserve((m*n*0.1 < 2*m) ? (int)(m*n*0.1) : 2*m);
00371 
00372       p[0] = 0;
00373       int nnz = 0;
00374 
00375       for(int j=0 ; j<n ; j++) {
00376         for(int i=0 ; i<m ; i++) {
00377           if(A[i+lb1][j+lb2] != 0.0) {
00378              ind.push_back(i);
00379              x.push_back(A[i+lb1][j+lb2]);
00380              nnz++;
00381           }
00382         }
00383           
00384         p[j+1] = nnz;
00385       }
00386 
00387     }
00388 
00389 
00391 
00394     srmatrix(const int ms, const int ns, const rmatrix& A) : m(ms), n(ns), lb1(1), ub1(ms), lb2(1), ub2(ns)  {
00395       //Banded matrix constructor
00396       int nnz = RowLen(A)*ColLen(A);
00397       p = std::vector<int>((n>0) ? n+1 : 1, 0);
00398       ind.reserve(nnz);
00399       x.reserve(nnz);
00400 
00401       std::vector<triplet_store<real> > work;
00402       work.reserve(nnz);
00403 
00404       
00405       for(int i=0 ; i<ColLen(A) ; i++) {
00406         for(int j=Lb(A,2) ; j<=Ub(A,2) ; j++) {
00407           if(i+j >=0  &&  i+j < n) {
00408             work.push_back(triplet_store<real>(i,i+j,A[i+Lb(A,1)][j]));
00409           }
00410         }
00411       }
00412 
00413       sort(work.begin(), work.end());
00414 
00415       int i=0;
00416 
00417       for(int j=0 ; j<n ; j++) {        
00418 
00419         while((unsigned int)i < work.size() && work[i].col == j ) {
00420           ind.push_back(work[i].row);
00421           x.push_back(work[i].val);
00422           i++;
00423         }
00424 
00425         p[j+1] = i;
00426       }
00427 
00428     }
00429 
00431     srmatrix(const srmatrix_slice&);
00432 
00434     void full(rmatrix& A) const {
00435        A = rmatrix(lb1,ub1,lb2,ub2);
00436        A = 0.0;
00437        for(int j=0 ; j<n ; j++) {
00438           for(int k=p[j] ; k<p[j+1] ; k++) {
00439              A[ind[k]+lb1][j+lb2] = x[k];
00440           }
00441        }
00442     }
00443 
00445 
00449     void dropzeros() {
00450       std::vector<int> pnew(n+1,0);
00451       std::vector<int> indnew;
00452       std::vector<real> xnew;
00453       int nnznew = 0;
00454 
00455       for(int j=0 ; j<n ; j++) {
00456         for(int k=p[j] ; k<p[j+1] ; k++) {
00457           if(x[k] != 0.0) {
00458             xnew.push_back(x[k]);
00459             indnew.push_back(ind[k]);
00460             nnznew++;
00461           }
00462         }
00463         pnew[j+1] = nnznew;
00464       }
00465 
00466       p = pnew;
00467       ind = indnew;
00468       x = xnew;
00469     }
00470 
00471 
00473     srmatrix& operator=(const real& A) {
00474       return sp_ms_assign<srmatrix,real,real>(*this,A);
00475     }
00476 
00478     srmatrix& operator=(const rmatrix& A) {
00479       return spf_mm_assign<srmatrix,rmatrix,real>(*this,A);
00480     }
00481 
00483     srmatrix& operator=(const rmatrix_slice& A) {
00484       return spf_mm_assign<srmatrix,rmatrix,real>(*this,A);
00485     }
00486 
00487 /*    srmatrix& operator=(const srmatrix& A) {
00488       p = A.p;
00489       ind = A.ind;
00490       x = A.x;
00491       return *this;
00492     }*/
00493 
00495     srmatrix& operator=(const srmatrix_slice&);
00496 
00498 
00504     const real operator()(int i, int j) const {
00505 #if(CXSC_INDEX_CHECK)
00506       if(i<lb1 || i>ub1 || j<lb2 || j>ub2)
00507         cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix::operator()(int, int)"));
00508 #endif
00509       real r = 0.0;
00510       for(int k=p[j-lb2] ; k<p[j-lb2+1] && ind[k]<=i-lb1 ; k++) {
00511         if(ind[k] == i-lb1)  r = x[k];
00512       }
00513       return r;
00514     }
00515 
00517 
00525     real& element(int i, int j) {
00526 #if(CXSC_INDEX_CHECK)
00527       if(i<lb1 || i>ub1 || j<lb2 || j>ub2)
00528         cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix::element(int, int)"));
00529 #endif
00530       int k;
00531       for(k=p[j-lb2] ; k<p[j-lb2+1] && ind[k]<=i-lb1 ; k++) {
00532         if(ind[k] == i-lb1)  return x[k];
00533       }
00534 
00535       //Nicht gefunden, Element muss angelegt werden, da Schreibzugriff moeglich
00536       std::vector<int>::iterator ind_it = ind.begin() + k;
00537       std::vector<real>::iterator x_it  = x.begin() + k;
00538       ind.insert(ind_it, i-lb1);
00539       x_it = x.insert(x_it, 0.0);
00540       for(k=j-lb2+1 ; k<(int)p.size() ; k++)
00541         p[k]++;
00542 
00543       return *x_it;
00544     }
00545 
00547     srmatrix_subv operator[](const cxscmatrix_column&);
00549     srmatrix_subv operator[](const int);
00551     const srmatrix_subv operator[](const cxscmatrix_column&) const;
00553     const srmatrix_subv operator[](const int) const;
00554 
00556     srmatrix_slice operator()(const int, const int , const int, const int);
00558     const srmatrix_slice operator()(const int, const int , const int, const int) const;
00559 
00561     srmatrix operator()(const intvector& pervec, const intvector& q) {
00562       srmatrix A(m,n,get_nnz());
00563       intvector per = perminv(pervec);
00564 
00565       int nnz=0;
00566       for(int k=0 ; k<n ; k++) {
00567         A.p[k] = nnz;
00568 
00569         std::map<int,real> work;
00570         for(int j=p[q[Lb(q)+k]] ; j<p[q[Lb(q)+k]+1] ; j++) 
00571            work.insert(std::make_pair(per[Lb(per)+ind[j]], x[j]));
00572         
00573         for(std::map<int,real>::iterator it = work.begin() ; it != work.end() ; it++) {
00574            A.ind.push_back(it->first);
00575            A.x.push_back(it->second);
00576         }
00577 
00578         nnz += work.size();
00579  
00580       }
00581 
00582       A.p[n] = nnz;
00583 
00584       return A;
00585     }
00586 
00588     srmatrix operator()(const intvector& pervec) {
00589       srmatrix A(m,n,get_nnz());
00590       intvector per = perminv(pervec);
00591 
00592       for(int k=0 ; k<n ; k++) {
00593         A.p[k] = p[k];
00594 
00595         std::map<int,real> work;
00596         for(int j=p[k] ; j<p[k+1] ; j++) 
00597            work.insert(std::make_pair(per[Lb(per)+ind[j]], x[j]));
00598         
00599         for(std::map<int,real>::iterator it = work.begin() ; it != work.end() ; it++) {
00600            A.ind.push_back(it->first);
00601            A.x.push_back(it->second);
00602         }
00603  
00604       }
00605 
00606       A.p[n] = p[n];
00607 
00608       return A;
00609     }
00610 
00612     srmatrix operator()(const intmatrix& P, const intmatrix& Q) {
00613       intvector p = permvec(P);
00614       intvector q = perminv(permvec(Q));
00615       return (*this)(p,q);
00616     }
00617 
00619     srmatrix operator()(const intmatrix& P) {
00620       intvector p = permvec(P);
00621       return (*this)(p);
00622     }
00623 
00625     real density() const {
00626       return p[n]/((double)m*n);
00627     }
00628 
00630     int get_nnz() const {
00631       return p[n];
00632     }
00633 
00635     srmatrix& operator+=(const rmatrix& B) {
00636       return spf_mm_addassign<srmatrix,rmatrix,rmatrix>(*this,B);
00637     }
00638 
00640     srmatrix& operator+=(const rmatrix_slice& B) {
00641       return spf_mm_addassign<srmatrix,rmatrix_slice,rmatrix>(*this,B);
00642     }
00643 
00645     srmatrix& operator+=(const srmatrix& B) {
00646       return spsp_mm_addassign<srmatrix,srmatrix,real>(*this,B);
00647     }
00648 
00650     srmatrix& operator-=(const rmatrix& B) {
00651       return spf_mm_subassign<srmatrix,rmatrix,rmatrix>(*this,B);
00652     }
00653 
00655     srmatrix& operator-=(const rmatrix_slice& B) {
00656       return spf_mm_subassign<srmatrix,rmatrix_slice,rmatrix>(*this,B);
00657     }
00658 
00660     srmatrix& operator-=(const srmatrix& B) {
00661       return spsp_mm_subassign<srmatrix,srmatrix,real>(*this,B);
00662     }
00663 
00665     srmatrix& operator*=(const rmatrix& B) {
00666       return spf_mm_multassign<srmatrix,rmatrix,sparse_dot,rmatrix>(*this,B);
00667     }
00668 
00670     srmatrix& operator*=(const rmatrix_slice& B) {
00671       return spf_mm_multassign<srmatrix,rmatrix,sparse_dot,rmatrix>(*this,B);
00672     }
00673 
00675     srmatrix& operator*=(const srmatrix& B) {
00676       return spsp_mm_multassign<srmatrix,srmatrix,sparse_dot,real>(*this,B);
00677     }
00678 
00680     srmatrix& operator*=(const real& r) {
00681       return sp_ms_multassign(*this,r);
00682     }
00683 
00685     srmatrix& operator/=(const real& r) {
00686       return sp_ms_divassign(*this,r);
00687     }
00688 
00689     friend int Lb(const srmatrix&, int);
00690     friend int Ub(const srmatrix&, int);
00691     friend void SetLb(srmatrix&, const int, const int);
00692     friend void SetUb(srmatrix&, const int, const int);
00693     friend int RowLen(const srmatrix&);
00694     friend int ColLen(const srmatrix&);
00695     friend srmatrix Re(const scmatrix&);
00696     friend srmatrix Im(const scmatrix&);
00697     friend srmatrix Inf(const simatrix&);
00698     friend srmatrix Sup(const simatrix&);
00699     friend srmatrix InfRe(const scimatrix&);
00700     friend srmatrix SupRe(const scimatrix&);
00701     friend srmatrix InfIm(const scimatrix&);
00702     friend srmatrix SupIm(const scimatrix&);
00703     friend srmatrix mid(const simatrix&);
00704     friend srmatrix diam(const simatrix&);
00705     friend srmatrix absmin(const simatrix&);
00706     friend srmatrix absmax(const simatrix&);
00707     friend srmatrix abs(const srmatrix&);
00708     friend srmatrix abs(const scmatrix&);    
00709 
00710     friend srmatrix CompMat(const simatrix&);
00711     friend srmatrix CompMat(const srmatrix&);
00712     friend srmatrix CompMat(const scmatrix&);
00713     friend srmatrix CompMat(const scimatrix&);    
00714     friend srmatrix transp(const srmatrix&);
00715     friend srmatrix Id(const srmatrix&);
00716 
00717     friend std::istream& operator>>(std::istream&, srmatrix_slice&);
00718     friend std::istream& operator>>(std::istream&, srmatrix_subv&);
00719 
00720     friend class srmatrix_slice;
00721     friend class srmatrix_subv;
00722     friend class srvector;
00723     friend class scmatrix;
00724     friend class scmatrix_slice;
00725     friend class scmatrix_subv;
00726     friend class simatrix;
00727     friend class simatrix_slice;
00728     friend class simatrix_subv;
00729     friend class scivector;
00730     friend class scimatrix;
00731     friend class scimatrix_slice;
00732     friend class scimatrix_subv;
00733     friend class rmatrix;
00734     friend class imatrix;
00735     friend class cmatrix;
00736     friend class cimatrix;
00737 
00738 #include "matrix_friend_declarations.inl"
00739 };
00740 
00741 inline rmatrix::rmatrix(const srmatrix& A) {
00742   dat = new real[A.m*A.n];
00743   lb1 = A.lb1; lb2 = A.lb2; ub1 = A.ub1; ub2 = A.ub2;
00744   xsize = A.n;
00745   ysize = A.m;
00746   *this = 0.0;
00747 
00748   for(int j=0 ; j<A.n ; j++) {
00749      for(int k=A.p[j] ; k<A.p[j+1] ; k++) {
00750         dat[A.ind[k]*A.n+j] = A.x[k];
00751      }
00752   }
00753 }
00754 
00756 inline srmatrix Id(const srmatrix& A) {
00757   srmatrix I(A.m, A.n, (A.m>A.n) ? A.m : A.n);
00758   I.lb1 = A.lb1; I.lb2 = A.lb2;
00759   I.ub1 = A.ub1; I.ub2 = A.ub2;
00760 
00761   if(A.m < A.n) {
00762     for(int i=0 ; i<A.m ; i++) {
00763       I.p[i+1] = I.p[i] + 1;
00764       I.ind.push_back(i);
00765       I.x.push_back(1.0);
00766     }
00767   } else {
00768     for(int i=0 ; i<A.n ; i++) {
00769       I.p[i+1] = I.p[i] + 1;
00770       I.ind.push_back(i);
00771       I.x.push_back(1.0);
00772     }
00773   }
00774 
00775   return I;
00776 }
00777 
00779 inline srmatrix transp(const srmatrix& A) {
00780   srmatrix B(A.n, A.m, A.get_nnz());
00781 
00782   //NIchtnullen pro Zeile bestimmen
00783   std::vector<int> w(A.m,0);
00784   for(unsigned int i=0 ; i<A.ind.size() ; i++) 
00785     w[A.ind[i]]++;
00786 
00787   //Spalten"pointer" setzen
00788   B.p.resize(A.m+1);
00789   B.p[0] = 0;
00790   for(unsigned int i=1 ; i<B.p.size() ; i++)
00791     B.p[i] = w[i-1] + B.p[i-1];
00792 
00793   //w vorbereiten
00794   w.insert(w.begin(), 0); 
00795   for(unsigned int i=1 ; i<w.size() ; i++) {
00796     w[i] += w[i-1];
00797   }
00798 
00799   //neuer zeilenindex und wert wird gesetzt
00800   int q;
00801   B.ind.resize(A.get_nnz());
00802   B.x.resize(A.get_nnz());
00803   for(int j=0 ; j<A.n ; j++) {
00804     for(int k=A.p[j] ; k<A.p[j+1] ; k++) {
00805       q = w[A.ind[k]]++;
00806       B.ind[q] = j;
00807       B.x[q] = A.x[k];
00808     }
00809   }
00810 
00811   return B;
00812 }
00813 
00815 inline srmatrix abs(const srmatrix& A) {
00816   srmatrix ret(A);
00817   for(unsigned int i=0 ; i<ret.x.size() ; i++) 
00818     ret.x[i] = abs(ret.x[i]);
00819   return ret;
00820 }
00821 
00823 inline srmatrix CompMat(const srmatrix& A) {
00824   srmatrix res(A.m,A.n,A.get_nnz());
00825   res.lb1 = A.lb1;
00826   res.lb2 = A.lb2;
00827   res.ub1 = A.ub1;
00828   res.ub2 = A.ub2;
00829   res.p   = A.p;
00830   res.ind = A.ind;
00831   res.p[A.n] = A.p[A.n];
00832 
00833   for(int j=0 ; j<res.n ; j++) {
00834     for(int k=A.p[j] ; k<A.p[j+1] ; k++) {
00835       if(A.ind[k] == j)
00836         res.x.push_back(abs(A.x[k]));
00837       else
00838         res.x.push_back(-abs(A.x[k]));
00839     }
00840   }
00841 
00842   res.dropzeros();
00843 
00844   return res; 
00845 }
00846 
00848 
00852 inline void SetLb(srmatrix& A, const int i, const int j) {
00853   if(i==1) {
00854     A.lb1 = j;
00855     A.ub1 = j + A.m - 1;
00856   } else if(i==2) {
00857     A.lb2 = j;
00858     A.ub2 = j + A.n - 1;
00859   }
00860 }
00861 
00863 
00867 inline void SetUb(srmatrix& A, const int i, const int j) {
00868   if(i==1) {
00869     A.ub1 = j;
00870     A.lb1 = j - A.m + 1;
00871   } else if(i==2) {
00872     A.ub2 = j;
00873     A.lb2 = j - A.n + 1;
00874   }
00875 }
00876 
00878 
00881 inline int Lb(const srmatrix& A, int i) {
00882   if(i==1) 
00883     return A.lb1;
00884   else if(i==2)
00885     return A.lb2;
00886   else
00887     return 1;
00888 }
00889 
00891 
00894 inline int Ub(const srmatrix& A, int i) {
00895   if(i==1) 
00896     return A.ub1;
00897   else if(i==2)
00898     return A.ub2;
00899   else
00900     return 1;
00901 }
00902 
00904 inline int RowLen(const srmatrix& A) {
00905   return A.n;
00906 }
00907 
00909 inline int ColLen(const srmatrix& A) {
00910   return A.m;
00911 }
00912 
00914 inline void Resize(srmatrix& A) {
00915   sp_m_resize(A);
00916 }
00917 
00919 inline void Resize(srmatrix& A, const int m, const int n) {
00920   sp_m_resize(A,m,n);
00921 }
00922 
00924 inline void Resize(srmatrix& A, const int l1, const int u1, const int l2, const int u2) {
00925   sp_m_resize(A,l1,u1,l2,u2);
00926 }
00927 
00929 
00935 inline rmatrix operator*(const rmatrix& A, const srmatrix& B) {
00936   return fsp_mm_mult<rmatrix,srmatrix,rmatrix,sparse_dot>(A,B);
00937 }
00938 
00940 
00946 inline rmatrix operator*(const srmatrix& A, const rmatrix& B) {
00947   return spf_mm_mult<srmatrix,rmatrix,rmatrix,sparse_dot>(A,B);
00948 }
00949 
00951 
00957 inline rmatrix operator*(const rmatrix_slice& A, const srmatrix& B) {
00958   return fsp_mm_mult<rmatrix_slice,srmatrix,rmatrix,sparse_dot>(A,B);
00959 }
00960 
00962 
00968 inline rmatrix operator*(const srmatrix& A, const rmatrix_slice& B) {
00969   return spf_mm_mult<srmatrix,rmatrix_slice,rmatrix,sparse_dot>(A,B);
00970 }
00971 
00973 
00979 inline srmatrix operator*(const srmatrix& A, const srmatrix& B) {
00980   return spsp_mm_mult<srmatrix,srmatrix,srmatrix,sparse_dot,real>(A,B);
00981 }
00982 
00984 inline srmatrix operator*(const srmatrix& A, const real& r) {
00985   return sp_ms_mult<srmatrix,real,srmatrix>(A,r);
00986 }
00987 
00989 inline srmatrix operator/(const srmatrix& A, const real& r) {
00990   return sp_ms_div<srmatrix,real,srmatrix>(A,r);
00991 }
00992 
00994 inline srmatrix operator*(const real& r, const srmatrix& A) {
00995   return sp_sm_mult<real,srmatrix,srmatrix>(r,A);
00996 }
00997 
00999 
01005 inline rvector operator*(const srmatrix& A, const rvector& v) {
01006   return spf_mv_mult<srmatrix,rvector,rvector,sparse_dot>(A,v);
01007 }
01008 
01010 
01016 inline rvector operator*(const srmatrix& A, const rvector_slice& v) {
01017   return spf_mv_mult<srmatrix,rvector_slice,rvector,sparse_dot>(A,v);
01018 }
01019 
01021 
01027 inline srvector operator*(const srmatrix& A, const srvector& v) {
01028   return spsp_mv_mult<srmatrix,srvector,srvector,sparse_dot,real>(A,v);
01029 }
01030 
01032 
01038 inline srvector operator*(const srmatrix& A, const srvector_slice& v) {
01039   return spsl_mv_mult<srmatrix,srvector_slice,srvector,sparse_dot,real>(A,v);
01040 }
01041 
01043 
01049 inline rvector operator*(const rmatrix& A, const srvector& v) {
01050   return fsp_mv_mult<rmatrix,srvector,rvector,sparse_dot>(A,v);
01051 }
01052 
01054 
01060 inline rvector operator*(const rmatrix_slice& A, const srvector& v) {
01061   return fsp_mv_mult<rmatrix_slice,srvector,rvector,sparse_dot>(A,v);
01062 }
01063 
01065 
01071 inline rvector operator*(const rmatrix& A, const srvector_slice& v) {
01072   return fsl_mv_mult<rmatrix,srvector_slice,rvector,sparse_dot>(A,v);
01073 }
01074 
01076 
01082 inline rvector operator*(const rmatrix_slice& A, const srvector_slice& v) {
01083   return fsl_mv_mult<rmatrix_slice,srvector_slice,rvector,sparse_dot>(A,v);
01084 }
01085 
01087 inline rmatrix operator+(const rmatrix& A, const srmatrix& B) {
01088   return fsp_mm_add<rmatrix,srmatrix,rmatrix>(A,B);
01089 }
01090 
01092 inline rmatrix operator+(const srmatrix& A, const rmatrix& B) {
01093   return spf_mm_add<srmatrix,rmatrix,rmatrix>(A,B);
01094 }
01095 
01097 inline rmatrix operator+(const rmatrix_slice& A, const srmatrix& B) {
01098   return fsp_mm_add<rmatrix_slice,srmatrix,rmatrix>(A,B);
01099 }
01100 
01102 inline rmatrix operator+(const srmatrix& A, const rmatrix_slice& B) {
01103   return spf_mm_add<srmatrix,rmatrix_slice,rmatrix>(A,B);
01104 }
01105 
01107 inline srmatrix operator+(const srmatrix& A, const srmatrix& B) {
01108   return spsp_mm_add<srmatrix,srmatrix,srmatrix,real>(A,B);
01109 }
01110 
01112 inline rmatrix operator-(const rmatrix& A, const srmatrix& B) {
01113   return fsp_mm_sub<rmatrix,srmatrix,rmatrix>(A,B);
01114 }
01115 
01117 inline rmatrix operator-(const srmatrix& A, const rmatrix& B) {
01118   return spf_mm_sub<srmatrix,rmatrix,rmatrix>(A,B);
01119 }
01120 
01122 inline rmatrix operator-(const rmatrix_slice& A, const srmatrix& B) {
01123   return fsp_mm_sub<rmatrix_slice,srmatrix,rmatrix>(A,B);
01124 }
01125 
01127 inline rmatrix operator-(const srmatrix& A, const rmatrix_slice& B) {
01128   return spf_mm_sub<srmatrix,rmatrix_slice,rmatrix>(A,B);
01129 }
01130 
01132 inline srmatrix operator-(const srmatrix& A, const srmatrix& B) {
01133   return spsp_mm_sub<srmatrix,srmatrix,srmatrix,real>(A,B);
01134 }
01135 
01137 inline srmatrix operator-(const srmatrix& M) {
01138   return sp_m_negative<srmatrix,srmatrix>(M);
01139 }
01140 
01142 inline srmatrix& operator+(srmatrix& A) {
01143   return A;
01144 }
01145 
01146 inline rmatrix& rmatrix::operator=(const srmatrix& B) {
01147   *this = rmatrix(B);
01148   return *this;
01149 }
01150 
01151 inline rmatrix_slice& rmatrix_slice::operator=(const srmatrix& B) {
01152   *this = rmatrix(B);
01153   return *this;
01154 }
01155 
01156 inline rmatrix& rmatrix::operator+=(const srmatrix& B) {
01157   return fsp_mm_addassign(*this,B);
01158 }
01159 
01160 inline rmatrix_slice& rmatrix_slice::operator+=(const srmatrix& B) {
01161   return fsp_mm_addassign(*this,B);
01162 }
01163 
01164 inline rmatrix& rmatrix::operator-=(const srmatrix& B) {
01165   return fsp_mm_subassign(*this,B);
01166 }
01167 
01168 inline rmatrix_slice& rmatrix_slice::operator-=(const srmatrix& B) {
01169   return fsp_mm_subassign(*this,B);
01170 }
01171 
01172 inline rmatrix& rmatrix::operator*=(const srmatrix& B) {
01173   return fsp_mm_multassign<rmatrix,srmatrix,sparse_dot,rmatrix>(*this,B);
01174 }
01175 
01176 inline rmatrix_slice& rmatrix_slice::operator*=(const srmatrix& B) {
01177   return fsp_mm_multassign<rmatrix_slice,srmatrix,sparse_dot,rmatrix>(*this,B);
01178 }
01179 
01181 inline bool operator==(const srmatrix& A, const srmatrix& B) {
01182   return spsp_mm_comp(A,B);
01183 }
01184 
01186 inline bool operator==(const srmatrix& A, const rmatrix& B) {
01187   return spf_mm_comp(A,B);
01188 }
01189 
01191 inline bool operator==(const rmatrix& A, const srmatrix& B) {
01192   return fsp_mm_comp(A,B);
01193 }
01194 
01196 inline bool operator==(const rmatrix_slice& A, const srmatrix& B) {
01197   return fsp_mm_comp(A,B);
01198 }
01199 
01201 inline bool operator==(const srmatrix& A, const rmatrix_slice& B) {
01202   return spf_mm_comp(A,B);
01203 }
01204 
01206 inline bool operator!=(const srmatrix& A, const srmatrix& B) {
01207   return !spsp_mm_comp(A,B);
01208 }
01209 
01211 inline bool operator!=(const srmatrix& A, const rmatrix& B) {
01212   return !spf_mm_comp(A,B);
01213 }
01214 
01216 inline bool operator!=(const rmatrix& A, const srmatrix& B) {
01217   return !fsp_mm_comp(A,B);
01218 }
01219 
01221 inline bool operator!=(const rmatrix_slice& A, const srmatrix& B) {
01222   return !fsp_mm_comp(A,B);
01223 }
01224 
01226 inline bool operator!=(const srmatrix& A, const rmatrix_slice& B) {
01227   return !spf_mm_comp(A,B);
01228 }
01229 
01231 inline bool operator<(const srmatrix& A, const srmatrix& B) {
01232   return spsp_mm_less<srmatrix,srmatrix,real>(A,B);
01233 }
01234 
01236 inline bool operator<(const srmatrix& A, const rmatrix& B) {
01237   return spf_mm_less<srmatrix,rmatrix,real>(A,B);
01238 }
01239 
01241 inline bool operator<(const rmatrix& A, const srmatrix& B) {
01242   return fsp_mm_less<rmatrix,srmatrix,real>(A,B);
01243 }
01244 
01246 inline bool operator<(const rmatrix_slice& A, const srmatrix& B) {
01247   return fsp_mm_less<rmatrix_slice,srmatrix,real>(A,B);
01248 }
01249 
01251 inline bool operator<(const srmatrix& A, const rmatrix_slice& B) {
01252   return spf_mm_less<srmatrix,rmatrix_slice,real>(A,B);
01253 }
01254 
01256 inline bool operator<=(const srmatrix& A, const srmatrix& B) {
01257   return spsp_mm_leq<srmatrix,srmatrix,real>(A,B);
01258 }
01259 
01261 inline bool operator<=(const srmatrix& A, const rmatrix& B) {
01262   return spf_mm_leq<srmatrix,rmatrix,real>(A,B);
01263 }
01264 
01266 inline bool operator<=(const rmatrix& A, const srmatrix& B) {
01267   return fsp_mm_leq<rmatrix,srmatrix,real>(A,B);
01268 }
01269 
01271 inline bool operator<=(const rmatrix_slice& A, const srmatrix& B) {
01272   return fsp_mm_leq<rmatrix_slice,srmatrix,real>(A,B);
01273 }
01274 
01276 inline bool operator<=(const srmatrix& A, const rmatrix_slice& B) {
01277   return spf_mm_leq<srmatrix,rmatrix_slice,real>(A,B);
01278 }
01279 
01281 inline bool operator>(const srmatrix& A, const srmatrix& B) {
01282   return spsp_mm_greater<srmatrix,srmatrix,real>(A,B);
01283 }
01284 
01286 inline bool operator>(const srmatrix& A, const rmatrix& B) {
01287   return spf_mm_greater<srmatrix,rmatrix,real>(A,B);
01288 }
01289 
01291 inline bool operator>(const rmatrix& A, const srmatrix& B) {
01292   return fsp_mm_greater<rmatrix,srmatrix,real>(A,B);
01293 }
01294 
01296 inline bool operator>(const rmatrix_slice& A, const srmatrix& B) {
01297   return fsp_mm_greater<rmatrix_slice,srmatrix,real>(A,B);
01298 }
01299 
01301 inline bool operator>(const srmatrix& A, const rmatrix_slice& B) {
01302   return spf_mm_greater<srmatrix,rmatrix_slice,real>(A,B);
01303 }
01304 
01306 inline bool operator>=(const srmatrix& A, const srmatrix& B) {
01307   return spsp_mm_geq<srmatrix,srmatrix,real>(A,B);
01308 }
01309 
01311 inline bool operator>=(const srmatrix& A, const rmatrix& B) {
01312   return spf_mm_geq<srmatrix,rmatrix,real>(A,B);
01313 }
01314 
01316 inline bool operator>=(const rmatrix& A, const srmatrix& B) {
01317   return fsp_mm_geq<rmatrix,srmatrix,real>(A,B);
01318 }
01319 
01321 inline bool operator>=(const rmatrix_slice& A, const srmatrix& B) {
01322   return fsp_mm_geq<rmatrix_slice,srmatrix,real>(A,B);
01323 }
01324 
01326 inline bool operator>=(const srmatrix& A, const rmatrix_slice& B) {
01327   return spf_mm_geq<srmatrix,rmatrix_slice,real>(A,B);
01328 }
01329 
01331 inline bool operator!(const srmatrix& A) {
01332   return sp_m_not(A);
01333 }
01334 
01336 
01341 inline std::ostream& operator<<(std::ostream& os, const srmatrix& A) {
01342   return sp_m_output<srmatrix,real>(os,A);
01343 }
01344 
01346 
01351 inline std::istream& operator>>(std::istream& is, srmatrix& A) {
01352   return sp_m_input<srmatrix,real>(is,A);
01353 }
01354 
01356 
01360 class srmatrix_slice {
01361   public:
01362     srmatrix  A;
01363     srmatrix* M; //Originalmatrix
01364 
01365   private:
01366     srmatrix_slice(srmatrix& Mat, int sl1l, int sl1u, int sl2l, int sl2u) {    
01367         A.lb1 = sl1l;
01368         A.lb2 = sl2l;
01369         A.ub1 = sl1u;
01370         A.ub2 = sl2u;
01371         A.m   = sl1u-sl1l+1;
01372         A.n   = sl2u-sl2l+1;
01373 
01374         //Kopieren der Werte aus A
01375         A.p = std::vector<int>(A.n+1, 0);
01376         A.ind.reserve(A.m + A.n);
01377         A.x.reserve(A.m + A.n);
01378 
01379         for(int i=0 ; i<A.n ; i++) {
01380            A.p[i+1] = A.p[i];
01381            for(int j=Mat.p[sl2l-Mat.lb2+i] ; j<Mat.p[sl2l-Mat.lb2+i+1] ; j++) {
01382               if(Mat.ind[j] >= sl1l-Mat.lb1  &&  Mat.ind[j] <= sl1u-Mat.lb1) {
01383                 A.ind.push_back(Mat.ind[j]-(sl1l-Mat.lb1));
01384                 A.x.push_back(Mat.x[j]);
01385                 A.p[i+1]++;
01386               }
01387            }
01388         }
01389 
01390         //Zeiger auf A fuer Datenmanipulationen
01391         M = &Mat;
01392     }
01393 
01394     srmatrix_slice(const srmatrix& Mat, int sl1l, int sl1u, int sl2l, int sl2u) {    
01395         A.lb1 = sl1l;
01396         A.lb2 = sl2l;
01397         A.ub1 = sl1u;
01398         A.ub2 = sl2u;
01399         A.m   = sl1u-sl1l+1;
01400         A.n   = sl2u-sl2l+1;
01401 
01402         //Kopieren der Werte aus A
01403         A.p = std::vector<int>(A.n+1, 0);
01404         A.ind.reserve(A.m + A.n);
01405         A.x.reserve(A.m + A.n);
01406 
01407         for(int i=0 ; i<A.n ; i++) {
01408            A.p[i+1] = A.p[i];
01409            for(int j=Mat.p[sl2l-Mat.lb2+i] ; j<Mat.p[sl2l-Mat.lb2+i+1] ; j++) {
01410               if(Mat.ind[j] >= sl1l-Mat.lb1  &&  Mat.ind[j] <= sl1u-Mat.lb1) {
01411                 A.ind.push_back(Mat.ind[j]-(sl1l-Mat.lb1));
01412                 A.x.push_back(Mat.x[j]);
01413                 A.p[i+1]++;
01414               }
01415            }
01416         }
01417 
01418         //Zeiger auf A fuer Datenmanipulationen
01419         M = const_cast<srmatrix*>(&Mat);
01420     }
01421 
01422   public:
01424     srmatrix_slice& operator=(const real& C) {
01425       return sl_ms_assign<srmatrix_slice, real, std::vector<real>::iterator, real>(*this,C);
01426     }
01427 
01429     srmatrix_slice& operator=(const srmatrix& C) {
01430       return slsp_mm_assign<srmatrix_slice, srmatrix, std::vector<real>::iterator>(*this,C);
01431     }
01432 
01434     srmatrix_slice& operator=(const rmatrix& C) {
01435       return slf_mm_assign<srmatrix_slice, rmatrix, std::vector<real>::iterator, real>(*this,C);
01436     }
01437 
01439     srmatrix_slice& operator=(const rmatrix_slice& C) {
01440       return slf_mm_assign<srmatrix_slice, rmatrix_slice, std::vector<real>::iterator, real>(*this,C);
01441     }
01442 
01444     srmatrix_slice& operator=(const srmatrix_slice& C) {
01445       *this = C.A;
01446       return *this;
01447     }
01448 
01450     srmatrix_slice& operator*=(const srmatrix_slice& M) {
01451       *this = A*M.A;
01452       return *this;
01453     }
01454 
01456     srmatrix_slice& operator*=(const srmatrix& M) {
01457       *this = A*M;
01458       return *this;
01459     }
01460 
01462     srmatrix_slice& operator*=(const rmatrix& M) {
01463       *this = A*M;
01464       return *this;
01465     }
01466 
01468     srmatrix_slice& operator*=(const rmatrix_slice& M) {
01469       *this = A*M;
01470       return *this;
01471     }
01472 
01474     srmatrix_slice& operator*=(const real& r) {
01475       *this = A*r;
01476       return *this;
01477     }
01478 
01480     srmatrix_slice& operator/=(const real& r) {
01481       *this = A/r;
01482       return *this;
01483     }
01484 
01486     srmatrix_slice& operator+=(const srmatrix_slice& M) {
01487       *this = A+M.A;
01488       return *this;
01489     } 
01490 
01492     srmatrix_slice& operator+=(const srmatrix& M) {
01493       *this = A+M;
01494       return *this;
01495     } 
01496 
01498     srmatrix_slice& operator+=(const rmatrix& M) {
01499       *this = A+M;
01500       return *this;
01501     } 
01502 
01504     srmatrix_slice& operator+=(const rmatrix_slice& M) {
01505       *this = A+M;
01506       return *this;
01507     } 
01508 
01510     srmatrix_slice& operator-=(const srmatrix_slice& M) {
01511       *this = A-M.A;
01512       return *this;
01513     } 
01514 
01516     srmatrix_slice& operator-=(const srmatrix& M) {
01517       *this = A-M;
01518       return *this;
01519     } 
01520 
01522     srmatrix_slice& operator-=(const rmatrix& M) {
01523       *this = A-M;
01524       return *this;
01525     } 
01526 
01528     srmatrix_slice& operator-=(const rmatrix_slice& M) {
01529       *this = A-M;
01530       return *this;
01531     }
01532 
01534 
01538     const real operator()(const int i, const int j) const {
01539 #if(CXSC_INDEX_CHECK)
01540       if(i<A.lb1 || i>A.ub1 || j<A.lb2 || j>A.ub2)
01541         cxscthrow(ELEMENT_NOT_IN_VEC("srmatrix_slice::operator()(int, int)"));
01542 #endif
01543       real r = A(i,j);
01544       return r;
01545     }
01546 
01548 
01552     real& element(const int i, const int j) {
01553 #if(CXSC_INDEX_CHECK)
01554       if(i<A.lb1 || i>A.ub1 || j<A.lb2 || j>A.ub2)
01555         cxscthrow(ELEMENT_NOT_IN_VEC("srmatrix::element(int, int)"));
01556 #endif
01557       return M->element(i,j);
01558     }
01559 
01561     srmatrix_subv operator[](const int);
01563     srmatrix_subv operator[](const cxscmatrix_column&);
01565     const srmatrix_subv operator[](const int) const;
01567     const srmatrix_subv operator[](const cxscmatrix_column&) const;
01568 
01569     friend int Lb(const srmatrix_slice&, const int);
01570     friend int Ub(const srmatrix_slice&, const int);
01571     friend int RowLen(const srmatrix_slice&);
01572     friend int ColLen(const srmatrix_slice&);
01573 
01574     friend class srmatrix;
01575     friend class srmatrix_subv;
01576     friend class srvector;
01577     friend class scmatrix;
01578     friend class scmatrix_slice;
01579     friend class scmatrix_subv;
01580     friend class scvector;
01581     friend class simatrix;
01582     friend class simatrix_slice;
01583     friend class simatrix_subv;
01584     friend class sivector;
01585     friend class scimatrix;
01586     friend class scimatrix_slice;
01587     friend class scimatrix_subv;
01588     friend class scivector;
01589     friend class rmatrix;
01590     friend class imatrix;
01591     friend class cmatrix;
01592     friend class cimatrix;
01593 
01594 #include "matrix_friend_declarations.inl"    
01595 };
01596 
01597 inline rmatrix::rmatrix(const srmatrix_slice& A) {
01598   dat = new real[A.A.m*A.A.n];
01599   lb1 = A.A.lb1; lb2 = A.A.lb2; ub1 = A.A.ub1; ub2 = A.A.ub2;
01600   xsize = A.A.n;
01601   ysize = A.A.m;
01602   *this = 0.0;
01603   for(int j=0 ; j<A.A.n ; j++) {
01604      for(int k=A.A.p[j] ; k<A.A.p[j+1] ; k++) {
01605         dat[A.A.ind[k]*A.A.n+j] = A.A.x[k];
01606      }
01607   }
01608 }
01609 
01611 inline int Lb(const srmatrix_slice& S, const int i) {
01612   return Lb(S.A, i);
01613 }
01614 
01616 inline int Ub(const srmatrix_slice& S, const int i) {
01617   return Ub(S.A, i);
01618 }
01619 
01621 inline int RowLen(const srmatrix_slice& S) {
01622   return RowLen(S.A);
01623 }
01624 
01626 inline int ColLen(const srmatrix_slice& S) {
01627   return ColLen(S.A);
01628 }
01629 
01630 inline srmatrix_slice srmatrix::operator()(const int i, const int j, const int k, const int l) {
01631 #if(CXSC_INDEX_CHECK)
01632   if(i<lb1 || j>ub1 || k<lb2 || l>ub2)
01633     cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix::operator()(int, int, int, int)"));
01634 #endif
01635   return srmatrix_slice(*this, i, j, k, l);
01636 }
01637 
01638 inline const srmatrix_slice srmatrix::operator()(const int i, const int j, const int k, const int l) const {
01639 #if(CXSC_INDEX_CHECK)
01640   if(i<lb1 || j>ub1 || k<lb2 || l>ub2)
01641     cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix::operator()(int, int, int, int) const"));
01642 #endif
01643   return srmatrix_slice(*this, i, j, k, l);
01644 }
01645 
01646 inline srmatrix& srmatrix::operator=(const srmatrix_slice& S) {
01647   *this = S.A;
01648   return *this;
01649 }
01650 
01651 inline srmatrix::srmatrix(const srmatrix_slice& S) {
01652   m = S.A.m;
01653   n = S.A.n;
01654   lb1 = S.A.lb1;
01655   ub1 = S.A.ub1;
01656   lb2 = S.A.lb2;
01657   ub2 = S.A.ub2;
01658   *this = S.A;
01659 }
01660 
01662 inline srmatrix operator-(const srmatrix_slice& M) {
01663   return sp_m_negative<srmatrix,srmatrix>(M.A);
01664 }
01665 
01667 inline srmatrix operator+(const srmatrix_slice& M) {
01668   return M.A;
01669 }
01670 
01672 
01678 inline srmatrix operator*(const srmatrix_slice& M1, const srmatrix_slice& M2) {
01679   return spsp_mm_mult<srmatrix,srmatrix,srmatrix,sparse_dot,real>(M1.A,M2.A);
01680 }
01681 
01683 
01689 inline srmatrix operator*(const srmatrix_slice& M1, const srmatrix& M2) {
01690   return spsp_mm_mult<srmatrix,srmatrix,srmatrix,sparse_dot,real>(M1.A,M2);
01691 }
01692 
01694 
01700 inline srmatrix operator*(const srmatrix& M1, const srmatrix_slice& M2) {
01701   return spsp_mm_mult<srmatrix,srmatrix,srmatrix,sparse_dot,real>(M1,M2.A);
01702 }
01703 
01705 
01711 inline rmatrix operator*(const srmatrix_slice& M1, const rmatrix& M2) {
01712   return spf_mm_mult<srmatrix,rmatrix,rmatrix,sparse_dot>(M1.A,M2);
01713 }
01714 
01716 
01722 inline rmatrix operator*(const rmatrix& M1, const srmatrix_slice& M2) {
01723   return fsp_mm_mult<rmatrix,srmatrix,rmatrix,sparse_dot>(M1,M2.A);
01724 }
01725 
01727 
01733 inline rmatrix operator*(const srmatrix_slice& M1, const rmatrix_slice& M2) {
01734   return spf_mm_mult<srmatrix,rmatrix_slice,rmatrix,sparse_dot>(M1.A,M2);
01735 }
01736 
01738 
01744 inline rmatrix operator*(const rmatrix_slice& M1, const srmatrix_slice& M2) {
01745   return fsp_mm_mult<rmatrix,srmatrix,rmatrix,sparse_dot>(M1,M2.A);
01746 }
01747 
01749 
01755 inline srvector operator*(const srmatrix_slice& M, const srvector& v) {
01756   return spsp_mv_mult<srmatrix,srvector,srvector,sparse_dot,real>(M.A,v);
01757 }
01758 
01760 
01766 inline srvector operator*(const srmatrix_slice& M, const srvector_slice& v) {
01767   return spsl_mv_mult<srmatrix,srvector_slice,srvector,sparse_dot,real>(M.A,v);
01768 }
01769 
01771 
01777 inline rvector operator*(const srmatrix_slice& M, const rvector& v) {
01778   return spf_mv_mult<srmatrix,rvector,rvector,sparse_dot>(M.A,v);
01779 }
01780 
01782 
01788 inline rvector operator*(const srmatrix_slice& M, const rvector_slice& v) {
01789   return spf_mv_mult<srmatrix,rvector_slice,rvector,sparse_dot>(M.A,v);
01790 }
01791 
01793 inline srmatrix operator*(const srmatrix_slice& M, const real& r) {
01794   return sp_ms_mult<srmatrix,real,srmatrix>(M.A,r);
01795 }
01796 
01798 inline srmatrix operator/(const srmatrix_slice& M, const real& r) {
01799   return sp_ms_div<srmatrix,real,srmatrix>(M.A,r);
01800 }
01801 
01803 inline srmatrix operator*(const real& r, const srmatrix_slice& M) {
01804   return sp_sm_mult<real,srmatrix,srmatrix>(r,M.A);
01805 }
01806 
01808 inline srmatrix operator+(const srmatrix_slice& M1, const srmatrix_slice& M2) {
01809   return spsp_mm_add<srmatrix,srmatrix,srmatrix,real>(M1.A,M2.A);
01810 }
01811 
01813 inline srmatrix operator+(const srmatrix_slice& M1, const srmatrix& M2) {
01814   return spsp_mm_add<srmatrix,srmatrix,srmatrix,real>(M1.A,M2);
01815 }
01816 
01818 inline srmatrix operator+(const srmatrix& M1, const srmatrix_slice& M2) {
01819   return spsp_mm_add<srmatrix,srmatrix,srmatrix,real>(M1,M2.A);
01820 }
01821 
01823 inline rmatrix operator+(const srmatrix_slice& M1, const rmatrix& M2) {
01824   return spf_mm_add<srmatrix,rmatrix,rmatrix>(M1.A,M2);
01825 }
01826 
01828 inline rmatrix operator+(const rmatrix& M1, const srmatrix_slice& M2) {
01829   return fsp_mm_add<rmatrix,srmatrix,rmatrix>(M1,M2.A);
01830 }
01831 
01833 inline rmatrix operator+(const srmatrix_slice& M1, const rmatrix_slice& M2) {
01834   return spf_mm_add<srmatrix,rmatrix_slice,rmatrix>(M1.A,M2);
01835 }
01836 
01838 inline rmatrix operator+(const rmatrix_slice& M1, const srmatrix_slice& M2) {
01839   return fsp_mm_add<rmatrix_slice,srmatrix,rmatrix>(M1,M2.A);
01840 }
01841 
01843 inline srmatrix operator-(const srmatrix_slice& M1, const srmatrix_slice& M2) {
01844   return spsp_mm_sub<srmatrix,srmatrix,srmatrix,real>(M1.A,M2.A);
01845 }
01846 
01848 inline srmatrix operator-(const srmatrix_slice& M1, const srmatrix& M2) {
01849   return spsp_mm_sub<srmatrix,srmatrix,srmatrix,real>(M1.A,M2);
01850 }
01851 
01853 inline srmatrix operator-(const srmatrix& M1, const srmatrix_slice& M2) {
01854   return spsp_mm_sub<srmatrix,srmatrix,srmatrix,real>(M1,M2.A);
01855 }
01856 
01858 inline rmatrix operator-(const srmatrix_slice& M1, const rmatrix& M2) {
01859   return spf_mm_sub<srmatrix,rmatrix,rmatrix>(M1.A,M2);
01860 }
01861 
01863 inline rmatrix operator-(const rmatrix& M1, const srmatrix_slice& M2) {
01864   return fsp_mm_sub<rmatrix,srmatrix,rmatrix>(M1,M2.A);
01865 }
01866 
01868 inline rmatrix operator-(const srmatrix_slice& M1, const rmatrix_slice& M2) {
01869   return spf_mm_sub<srmatrix,rmatrix_slice,rmatrix>(M1.A,M2);
01870 }
01871 
01873 inline rmatrix operator-(const rmatrix_slice& M1, const srmatrix_slice& M2) {
01874   return fsp_mm_sub<rmatrix_slice,srmatrix,rmatrix>(M1,M2.A);
01875 }
01876 
01877 inline rmatrix& rmatrix::operator=(const srmatrix_slice& M) {
01878   *this = rmatrix(M);
01879   return *this;
01880 }
01881 
01882 inline rmatrix& rmatrix::operator+=(const srmatrix_slice& M) {
01883   *this += M.A;
01884   return *this;
01885 }
01886 
01887 inline rmatrix_slice& rmatrix_slice::operator+=(const srmatrix_slice& M) {
01888   *this += M.A;
01889   return *this;
01890 }
01891 
01892 inline rmatrix& rmatrix::operator-=(const srmatrix_slice& M) {
01893   *this -= M.A;
01894   return *this;
01895 }
01896 
01897 inline rmatrix_slice& rmatrix_slice::operator-=(const srmatrix_slice& M) {
01898   *this -= M.A;
01899   return *this;
01900 }
01901 
01902 inline rmatrix& rmatrix::operator*=(const srmatrix_slice& M) {
01903   *this *= M.A;
01904   return *this;
01905 }
01906 
01907 inline rmatrix_slice& rmatrix_slice::operator*=(const srmatrix_slice& M) {
01908   *this *= M.A;
01909   return *this;
01910 }
01911 
01913 inline bool operator==(const srmatrix_slice& M1, const srmatrix_slice& M2) {
01914   return spsp_mm_comp(M1.A,M2.A);
01915 }
01916 
01918 inline bool operator==(const srmatrix_slice& M1, const srmatrix& M2) {
01919   return spsp_mm_comp(M1.A,M2);
01920 }
01921 
01923 inline bool operator==(const srmatrix& M1, const srmatrix_slice& M2) {
01924   return spsp_mm_comp(M1,M2.A);
01925 }
01926 
01928 inline bool operator==(const srmatrix_slice& M1, const rmatrix& M2) {
01929   return spf_mm_comp(M1.A,M2);
01930 }
01931 
01933 inline bool operator==(const rmatrix& M1, const srmatrix_slice& M2) {
01934   return fsp_mm_comp(M1,M2.A);
01935 }
01936 
01938 inline bool operator==(const rmatrix_slice& M1, const srmatrix_slice& M2) {
01939   return fsp_mm_comp(M1,M2.A);
01940 }
01941 
01943 inline bool operator==(const srmatrix_slice& M1, const rmatrix_slice& M2) {
01944   return spf_mm_comp(M1.A,M2);
01945 }
01946 
01948 inline bool operator!=(const srmatrix_slice& M1, const srmatrix_slice& M2) {
01949   return !spsp_mm_comp(M1.A,M2.A);
01950 }
01951 
01953 inline bool operator!=(const srmatrix_slice& M1, const srmatrix& M2) {
01954   return !spsp_mm_comp(M1.A,M2);
01955 }
01956 
01958 inline bool operator!=(const srmatrix& M1, const srmatrix_slice& M2) {
01959   return !spsp_mm_comp(M1,M2.A);
01960 }
01961 
01963 inline bool operator!=(const srmatrix_slice& M1, const rmatrix& M2) {
01964   return !spf_mm_comp(M1.A,M2);
01965 }
01966 
01968 inline bool operator!=(const rmatrix& M1, const srmatrix_slice& M2) {
01969   return !fsp_mm_comp(M1,M2.A);
01970 }
01971 
01973 inline bool operator!=(const rmatrix_slice& M1, const srmatrix_slice& M2) {
01974   return !fsp_mm_comp(M1,M2.A);
01975 }
01976 
01978 inline bool operator!=(const srmatrix_slice& M1, const rmatrix_slice& M2) {
01979   return !spf_mm_comp(M1.A,M2);
01980 }
01981 
01983 inline bool operator<(const srmatrix_slice& M1, const srmatrix_slice& M2) {
01984   return spsp_mm_less<srmatrix,srmatrix,real>(M1.A,M2.A);
01985 }
01986 
01988 inline bool operator<(const srmatrix_slice& M1, const srmatrix& M2) {
01989   return spsp_mm_less<srmatrix,srmatrix,real>(M1.A,M2);
01990 }
01991 
01993 inline bool operator<(const srmatrix& M1, const srmatrix_slice& M2) {
01994   return spsp_mm_less<srmatrix,srmatrix,real>(M1,M2.A);
01995 }
01996 
01998 inline bool operator<(const srmatrix_slice& M1, const rmatrix& M2) {
01999   return spf_mm_less<srmatrix,rmatrix,real>(M1.A,M2);
02000 }
02001 
02003 inline bool operator<(const rmatrix& M1, const srmatrix_slice& M2) {
02004   return fsp_mm_less<rmatrix,srmatrix,real>(M1,M2.A);
02005 }
02006 
02008 inline bool operator<(const rmatrix_slice& M1, const srmatrix_slice& M2) {
02009   return fsp_mm_less<rmatrix_slice,srmatrix,real>(M1,M2.A);
02010 }
02011 
02013 inline bool operator<(const srmatrix_slice& M1, const rmatrix_slice& M2) {
02014   return spf_mm_less<srmatrix,rmatrix_slice,real>(M1.A,M2);
02015 }
02016 
02018 inline bool operator<=(const srmatrix_slice& M1, const srmatrix_slice& M2) {
02019   return spsp_mm_leq<srmatrix,srmatrix,real>(M1.A,M2.A);
02020 }
02021 
02023 inline bool operator<=(const srmatrix_slice& M1, const srmatrix& M2) {
02024   return spsp_mm_leq<srmatrix,srmatrix,real>(M1.A,M2);
02025 }
02026 
02028 inline bool operator<=(const srmatrix& M1, const srmatrix_slice& M2) {
02029   return spsp_mm_leq<srmatrix,srmatrix,real>(M1,M2.A);
02030 }
02031 
02033 inline bool operator<=(const srmatrix_slice& M1, const rmatrix& M2) {
02034   return spf_mm_leq<srmatrix,rmatrix,real>(M1.A,M2);
02035 }
02036 
02038 inline bool operator<=(const rmatrix& M1, const srmatrix_slice& M2) {
02039   return fsp_mm_leq<rmatrix,srmatrix,real>(M1,M2.A);
02040 }
02041 
02043 inline bool operator<=(const rmatrix_slice& M1, const srmatrix_slice& M2) {
02044   return fsp_mm_leq<rmatrix_slice,srmatrix,real>(M1,M2.A);
02045 }
02046 
02048 inline bool operator<=(const srmatrix_slice& M1, const rmatrix_slice& M2) {
02049   return spf_mm_leq<srmatrix,rmatrix_slice,real>(M1.A,M2);
02050 }
02051 
02053 inline bool operator>(const srmatrix_slice& M1, const srmatrix_slice& M2) {
02054   return spsp_mm_greater<srmatrix,srmatrix,real>(M1.A,M2.A);
02055 }
02056 
02058 inline bool operator>(const srmatrix_slice& M1, const srmatrix& M2) {
02059   return spsp_mm_greater<srmatrix,srmatrix,real>(M1.A,M2);
02060 }
02061 
02063 inline bool operator>(const srmatrix& M1, const srmatrix_slice& M2) {
02064   return spsp_mm_greater<srmatrix,srmatrix,real>(M1,M2.A);
02065 }
02066 
02068 inline bool operator>(const srmatrix_slice& M1, const rmatrix& M2) {
02069   return spf_mm_greater<srmatrix,rmatrix,real>(M1.A,M2);
02070 }
02071 
02073 inline bool operator>(const rmatrix& M1, const srmatrix_slice& M2) {
02074   return fsp_mm_greater<rmatrix,srmatrix,real>(M1,M2.A);
02075 }
02076 
02078 inline bool operator>(const rmatrix_slice& M1, const srmatrix_slice& M2) {
02079   return fsp_mm_greater<rmatrix_slice,srmatrix,real>(M1,M2.A);
02080 }
02081 
02083 inline bool operator>(const srmatrix_slice& M1, const rmatrix_slice& M2) {
02084   return spf_mm_greater<srmatrix,rmatrix_slice,real>(M1.A,M2);
02085 }
02086 
02088 inline bool operator>=(const srmatrix_slice& M1, const srmatrix_slice& M2) {
02089   return spsp_mm_geq<srmatrix,srmatrix,real>(M1.A,M2.A);
02090 }
02091 
02093 inline bool operator>=(const srmatrix_slice& M1, const srmatrix& M2) {
02094   return spsp_mm_geq<srmatrix,srmatrix,real>(M1.A,M2);
02095 }
02096 
02098 inline bool operator>=(const srmatrix& M1, const srmatrix_slice& M2) {
02099   return spsp_mm_geq<srmatrix,srmatrix,real>(M1,M2.A);
02100 }
02101 
02103 inline bool operator>=(const srmatrix_slice& M1, const rmatrix& M2) {
02104   return spf_mm_geq<srmatrix,rmatrix,real>(M1.A,M2);
02105 }
02106 
02108 inline bool operator>=(const rmatrix& M1, const srmatrix_slice& M2) {
02109   return fsp_mm_geq<rmatrix,srmatrix,real>(M1,M2.A);
02110 }
02111 
02113 inline bool operator>=(const rmatrix_slice& M1, const srmatrix_slice& M2) {
02114   return fsp_mm_geq<rmatrix_slice,srmatrix,real>(M1,M2.A);
02115 }
02116 
02118 inline bool operator>=(const srmatrix_slice& M1, const rmatrix_slice& M2) {
02119   return spf_mm_geq<srmatrix,rmatrix_slice,real>(M1.A,M2);
02120 }
02121 
02123 inline bool operator!(const srmatrix_slice& M) {
02124   return sp_m_not(M.A);
02125 }
02126 
02128 
02133 inline std::ostream& operator<<(std::ostream& os, const srmatrix_slice& M) {
02134   return sp_m_output<srmatrix,real>(os, M.A);
02135 }
02136 
02138 
02143 inline std::istream& operator>>(std::istream& is, srmatrix_slice& M) {
02144   srmatrix tmp(M.A.m,M.A.n);
02145   sp_m_input<srmatrix,real>(is, tmp);
02146   M = tmp;
02147   return is;
02148 }
02149 
02150 
02152 
02157 class srmatrix_subv {
02158   private:
02159     srmatrix_slice dat;
02160     bool row;
02161     int index;
02162 
02163 
02164     srmatrix_subv(srmatrix& A, bool r, int i, int j, int k, int l) : dat(A,i,j,k,l), row(r) {
02165        if(row) index=i; else index=k;
02166     }
02167 
02168     srmatrix_subv(const srmatrix& A, bool r, int i, int j, int k, int l) : dat(A,i,j,k,l), row(r){
02169        if(row) index=i; else index=k;
02170     }
02171 
02172   public:
02173     
02175 
02179     real& operator[](const int i) {
02180       if(row) {
02181 #if(CXSC_INDEX_CHECK)
02182         if(i<dat.A.lb2 || i>dat.A.ub2)
02183           cxscthrow(ELEMENT_NOT_IN_VEC("srmatrix_subv::operator[](int)"));
02184 #endif
02185         return dat.element(index,i);
02186       } else {
02187 #if(CXSC_INDEX_CHECK)
02188         if(i<dat.A.lb1 || i>dat.A.ub1)
02189           cxscthrow(ELEMENT_NOT_IN_VEC("srmatrix_subv::operator[](int)"));
02190 #endif
02191         return dat.element(i,index);
02192       }
02193     }
02194 
02196 
02199     const real operator[](const int i) const {
02200       if(row) {
02201 #if(CXSC_INDEX_CHECK)
02202         if(i<dat.A.lb2 || i>dat.A.ub2)
02203           cxscthrow(ELEMENT_NOT_IN_VEC("srmatrix_subv::operator[](int)"));
02204 #endif
02205         return dat(index,i);
02206       } else {
02207 #if(CXSC_INDEX_CHECK)
02208         if(i<dat.A.lb1 || i>dat.A.ub1)
02209           cxscthrow(ELEMENT_NOT_IN_VEC("srmatrix_subv::operator[](int)"));
02210 #endif
02211         return dat(i,index);
02212       }
02213     }
02214 
02216     srmatrix_subv& operator=(const srmatrix_subv& v) {
02217       return svsp_vv_assign(*this,srvector(v));
02218     }
02219 
02221     srmatrix_subv& operator=(const real& v) {
02222       return sv_vs_assign(*this,v);
02223     }
02224 
02226     srmatrix_subv& operator=(const srvector& v) {
02227       return svsp_vv_assign(*this,v);
02228     }
02229 
02231     srmatrix_subv& operator=(const srvector_slice& v) {
02232       return svsl_vv_assign(*this,v);
02233     }
02234 
02236     srmatrix_subv& operator=(const rvector& v) {
02237       return svf_vv_assign(*this,v);
02238     }
02239 
02241     srmatrix_subv& operator=(const rvector_slice& v) {
02242       return svf_vv_assign(*this,v);
02243     }
02244 
02246     srmatrix_subv& operator*=(const real&);
02248     srmatrix_subv& operator/=(const real&);
02250     srmatrix_subv& operator+=(const srvector&);
02252     srmatrix_subv& operator+=(const srvector_slice&);
02254     srmatrix_subv& operator+=(const rvector&);
02256     srmatrix_subv& operator+=(const rvector_slice&);
02258     srmatrix_subv& operator-=(const srvector&);
02260     srmatrix_subv& operator-=(const srvector_slice&);
02262     srmatrix_subv& operator-=(const rvector&);
02264     srmatrix_subv& operator-=(const rvector_slice&);
02265 
02266     friend srvector operator-(const srmatrix_subv&);
02267 
02268     friend std::istream& operator>>(std::istream&, srmatrix_subv&);
02269 
02270     friend int Lb(const srmatrix_subv&);
02271     friend int Ub(const srmatrix_subv&);
02272     friend int VecLen(const srmatrix_subv&);
02273 
02274     friend class srvector;
02275     friend class srmatrix;
02276     friend class srmatrix_slice;
02277     friend class scvector;
02278     friend class scmatrix;
02279     friend class scmatrix_slice;
02280     friend class sivector;
02281     friend class simatrix;
02282     friend class simatrix_slice;
02283     friend class scivector;
02284     friend class scimatrix;
02285     friend class scimatrix_slice;
02286 
02287 #include "vector_friend_declarations.inl"
02288 };
02289 
02291 inline int Lb(const srmatrix_subv& S) {
02292   if(S.row)
02293     return Lb(S.dat, 2);
02294   else
02295     return Lb(S.dat, 1);
02296 }
02297 
02299 inline int Ub(const srmatrix_subv& S) {
02300   if(S.row)
02301     return Ub(S.dat, 2);
02302   else
02303     return Ub(S.dat, 1);
02304 }
02305 
02307 inline int VecLen(const srmatrix_subv& S) {
02308   return Ub(S)-Lb(S)+1;
02309 }
02310 
02312 inline std::ostream& operator<<(std::ostream& os, const srmatrix_subv& v) {
02313   os << srvector(v);
02314   return os;
02315 }
02316 
02318 inline std::istream& operator>>(std::istream& is, srmatrix_subv& v) {
02319   int n = 0;
02320   if(v.row) n=v.dat.A.n; else n=v.dat.A.m;
02321   srvector tmp(n);
02322   is >> tmp;
02323   v = tmp;
02324   return is;
02325 }
02326 
02327 inline srmatrix_subv srmatrix::operator[](const cxscmatrix_column& c) {
02328 #if(CXSC_INDEX_CHECK)
02329   if(c.col()<lb2 || c.col()>ub2)
02330     cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix::operator[](const cxscmatrix_column&)"));
02331 #endif
02332   return srmatrix_subv(*this, false, lb1, ub1, c.col(), c.col());
02333 }
02334 
02335 inline srmatrix_subv srmatrix::operator[](const int i) {
02336 #if(CXSC_INDEX_CHECK)
02337   if(i<lb1 || i>ub1)
02338     cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix::operator[](const int)"));
02339 #endif
02340   return srmatrix_subv(*this, true, i, i, lb2, ub2);
02341 }
02342 
02343 inline const srmatrix_subv srmatrix::operator[](const cxscmatrix_column& c) const {
02344 #if(CXSC_INDEX_CHECK)
02345   if(c.col()<lb2 || c.col()>ub2)
02346     cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix::operator[](const cxscmatrix_column&)"));
02347 #endif
02348   return srmatrix_subv(*this, false, lb1, ub1, c.col(), c.col());
02349 }
02350 
02351 inline const srmatrix_subv srmatrix::operator[](const int i) const {
02352 #if(CXSC_INDEX_CHECK)
02353   if(i<lb1 || i>ub1)
02354     cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix::operator[](const int)"));
02355 #endif
02356   return srmatrix_subv(*this, true, i, i, lb2, ub2);
02357 }
02358 
02359 inline srmatrix_subv srmatrix_slice::operator[](const int i) {
02360 #if(CXSC_INDEX_CHECK)
02361   if(i<A.lb1 || i>A.ub1)
02362     cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix_slice::operator[](const int)"));
02363 #endif
02364   return srmatrix_subv(*M, true, i, i, A.lb2, A.ub2);
02365 }
02366 
02367 inline srmatrix_subv srmatrix_slice::operator[](const cxscmatrix_column& c) {
02368 #if(CXSC_INDEX_CHECK)
02369   if(c.col()<A.lb2 || c.col()>A.ub2)
02370     cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix_slice::operator[](const cxscmatrix_column&)"));
02371 #endif
02372   return srmatrix_subv(*M, false, A.lb1, A.ub1, c.col(), c.col());
02373 }
02374 
02375 inline const srmatrix_subv srmatrix_slice::operator[](const int i) const {
02376 #if(CXSC_INDEX_CHECK)
02377   if(i<A.lb1 || i>A.ub1)
02378     cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix_slice::operator[](const int)"));
02379 #endif
02380   return srmatrix_subv(*M, true, i, i, A.lb2, A.ub2);
02381 }
02382 
02383 inline const srmatrix_subv srmatrix_slice::operator[](const cxscmatrix_column& c) const {
02384 #if(CXSC_INDEX_CHECK)
02385   if(c.col()<A.lb2 || c.col()>A.ub2)
02386     cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix_slice::operator[](const cxscmatrix_column&)"));
02387 #endif
02388   return srmatrix_subv(*M, false, A.lb1, A.ub1, c.col(), c.col());
02389 }
02390 
02391 inline srvector::srvector(const srmatrix_subv& A) {
02392   p.reserve(A.dat.A.get_nnz());
02393   x.reserve(A.dat.A.get_nnz());
02394 
02395   if(A.row) {
02396     lb = A.dat.A.lb2;
02397     ub = A.dat.A.ub2;
02398     n = ub-lb+1; 
02399 
02400     for(int j=0 ; j<n ; j++) {
02401       for(int k=A.dat.A.p[j] ; k<A.dat.A.p[j+1] ; k++) {
02402         p.push_back(j);
02403         x.push_back(A.dat.A.x[k]);
02404       }
02405     }
02406 
02407   } else {
02408     lb = A.dat.A.lb1;
02409     ub = A.dat.A.ub1;
02410     n = ub-lb+1; 
02411 
02412     for(unsigned int k=0 ; k<A.dat.A.ind.size() ; k++) {
02413         p.push_back(A.dat.A.ind[k]);
02414         x.push_back(A.dat.A.x[k]);
02415     }
02416   }
02417 }
02418 
02420 inline srvector operator-(const srmatrix_subv& v) {
02421  srvector s(v);
02422  return -s;
02423 }
02424 
02426 inline srvector operator*(const srmatrix_subv& v1, const real& v2) {
02427   return srvector(v1) * v2;
02428 }
02429 
02431 inline srvector operator/(const srmatrix_subv& v1, const real& v2) {
02432   return srvector(v1) / v2;
02433 }
02434 
02436 inline srvector operator*(const real& v1, const srmatrix_subv& v2) {
02437   return v1 * srvector(v2);
02438 }
02439 
02441 
02447 inline real operator*(const srmatrix_subv& v1, const srvector& v2) {
02448   return srvector(v1) * v2;
02449 }
02450 
02452 
02458 inline real operator*(const srmatrix_subv& v1, const srvector_slice& v2) {
02459   return srvector(v1) * v2;
02460 }
02461 
02463 
02469 inline real operator*(const srmatrix_subv& v1, const rvector& v2) {
02470   return srvector(v1) * v2;
02471 }
02472 
02474 
02480 inline real operator*(const srmatrix_subv& v1, const rvector_slice& v2) {
02481   return srvector(v1) * v2;
02482 }
02483 
02485 
02491 inline real operator*(const srvector& v1, const srmatrix_subv& v2) {
02492   return v1 * srvector(v2);
02493 }
02494 
02496 
02502 inline real operator*(const srvector_slice& v1, const srmatrix_subv& v2) {
02503   return v1 * srvector(v2);
02504 }
02505 
02507 
02513 inline real operator*(const rvector& v1, const srmatrix_subv& v2) {
02514   return v1 * srvector(v2);
02515 }
02516 
02518 
02524 inline real operator*(const rvector_slice& v1, const srmatrix_subv& v2) {
02525   return v1 * srvector(v2);
02526 }
02527 
02529 inline srvector operator+(const srmatrix_subv& v1, const srvector& v2) {
02530   return srvector(v1) + v2;
02531 }
02532 
02534 inline srvector operator+(const srmatrix_subv& v1, const srvector_slice& v2) {
02535   return srvector(v1) + v2;
02536 }
02537 
02539 inline rvector operator+(const srmatrix_subv& v1, const rvector& v2) {
02540   return srvector(v1) + v2;
02541 }
02542 
02544 inline rvector operator+(const srmatrix_subv& v1, const rvector_slice& v2) {
02545   return srvector(v1) + v2;
02546 }
02547 
02549 inline srvector operator+(const srvector& v1, const srmatrix_subv& v2) {
02550   return v1 + srvector(v2);
02551 }
02552 
02554 inline srvector operator+(const srvector_slice& v1, const srmatrix_subv& v2) {
02555   return v1 + srvector(v2);
02556 }
02557 
02559 inline rvector operator+(const rvector& v1, const srmatrix_subv& v2) {
02560   return v1 + srvector(v2);
02561 }
02562 
02564 inline rvector operator+(const rvector_slice& v1, const srmatrix_subv& v2) {
02565   return v1 + srvector(v2);
02566 }
02567 
02569 inline srvector operator-(const srmatrix_subv& v1, const srvector& v2) {
02570   return srvector(v1) - v2;
02571 }
02572 
02574 inline srvector operator-(const srmatrix_subv& v1, const srvector_slice& v2) {
02575   return srvector(v1) - v2;
02576 }
02577 
02579 inline rvector operator-(const srmatrix_subv& v1, const rvector& v2) {
02580   return srvector(v1) - v2;
02581 }
02582 
02584 inline rvector operator-(const srmatrix_subv& v1, const rvector_slice& v2) {
02585   return srvector(v1) - v2;
02586 }
02587 
02589 inline srvector operator-(const srvector& v1, const srmatrix_subv& v2) {
02590   return v1 - srvector(v2);
02591 }
02592 
02594 inline srvector operator-(const srvector_slice& v1, const srmatrix_subv& v2) {
02595   return v1 - srvector(v2);
02596 }
02597 
02599 inline rvector operator-(const rvector& v1, const srmatrix_subv& v2) {
02600   return v1 - srvector(v2);
02601 }
02602 
02604 inline rvector operator-(const rvector_slice& v1, const srmatrix_subv& v2) {
02605   return v1 - srvector(v2);
02606 }
02607 
02608 inline srmatrix_subv& srmatrix_subv::operator*=(const real& v) {
02609   *this = *this * v;
02610   return *this;
02611 }
02612 
02613 inline srmatrix_subv& srmatrix_subv::operator/=(const real& v) {
02614   *this = *this * v;
02615   return *this;
02616 }
02617 
02618 inline srmatrix_subv& srmatrix_subv::operator+=(const srvector& v) {
02619   *this = *this + v;
02620   return *this;
02621 }
02622 
02623 inline srmatrix_subv& srmatrix_subv::operator+=(const srvector_slice& v) {
02624   *this = *this + v;
02625   return *this;
02626 }
02627 
02628 inline srmatrix_subv& srmatrix_subv::operator+=(const rvector& v) {
02629   *this = *this + v;
02630   return *this;
02631 }
02632 
02633 inline srmatrix_subv& srmatrix_subv::operator+=(const rvector_slice& v) {
02634   *this = *this + v;
02635   return *this;
02636 }
02637 
02638 inline srmatrix_subv& srmatrix_subv::operator-=(const srvector& v) {
02639   *this = *this - v;
02640   return *this;
02641 }
02642 
02643 inline srmatrix_subv& srmatrix_subv::operator-=(const srvector_slice& v) {
02644   *this = *this - v;
02645   return *this;
02646 }
02647 
02648 inline srmatrix_subv& srmatrix_subv::operator-=(const rvector& v) {
02649   *this = *this - v;
02650   return *this;
02651 }
02652 
02653 inline srmatrix_subv& srmatrix_subv::operator-=(const rvector_slice& v) {
02654   *this = *this - v;
02655   return *this;
02656 }
02657 
02658 inline rmatrix_subv& rmatrix_subv::operator+=(const srmatrix_subv& v) {
02659   *this += rvector(v);
02660   return *this;
02661 }
02662 
02663 inline rmatrix_subv& rmatrix_subv::operator+=(const srvector& v) {
02664   *this += rvector(v);
02665   return *this;
02666 }
02667 
02668 inline rmatrix_subv& rmatrix_subv::operator+=(const srvector_slice& v) {
02669   *this += rvector(v);
02670   return *this;
02671 }
02672 
02673 inline rmatrix_subv& rmatrix_subv::operator-=(const srmatrix_subv& v) {
02674   *this -= rvector(v);
02675   return *this;
02676 }
02677 
02678 inline rmatrix_subv& rmatrix_subv::operator=(const srvector& v) {
02679   *this = rvector(v);
02680   return *this;
02681 }
02682 
02683 inline rmatrix_subv& rmatrix_subv::operator=(const srvector_slice& v) {
02684   *this = rvector(v);
02685   return *this;
02686 }
02687 
02688 inline rmatrix_subv& rmatrix_subv::operator=(const srmatrix_subv& v) {
02689   *this = rvector(v);
02690   return *this;
02691 }
02692 
02694 inline bool operator==(const srmatrix_subv& v1, const srvector& v2) {
02695   return srvector(v1) == v2;
02696 }
02697 
02699 inline bool operator==(const srmatrix_subv& v1, const srvector_slice& v2) {
02700   return srvector(v1) == v2;
02701 }
02702 
02704 inline bool operator==(const srmatrix_subv& v1, const rvector& v2) {
02705   return srvector(v1) == v2;
02706 }
02707 
02709 inline bool operator==(const srmatrix_subv& v1, const rvector_slice& v2) {
02710   return srvector(v1) == v2;
02711 }
02712 
02714 inline bool operator==(const srvector& v1, const srmatrix_subv& v2) {
02715   return v1 == srvector(v2);
02716 }
02717 
02719 inline bool operator==(const srvector_slice& v1, const srmatrix_subv& v2) {
02720   return v1 == srvector(v2);
02721 }
02722 
02724 inline bool operator==(const rvector& v1, const srmatrix_subv& v2) {
02725   return v1 == srvector(v2);
02726 }
02727 
02729 inline bool operator==(const rvector_slice& v1, const srmatrix_subv& v2) {
02730   return v1 == srvector(v2);
02731 }
02732 
02734 inline bool operator!=(const srmatrix_subv& v1, const srvector& v2) {
02735   return srvector(v1) != v2;
02736 }
02737 
02739 inline bool operator!=(const srmatrix_subv& v1, const srvector_slice& v2) {
02740   return srvector(v1) != v2;
02741 }
02742 
02744 inline bool operator!=(const srmatrix_subv& v1, const rvector& v2) {
02745   return srvector(v1) != v2;
02746 }
02747 
02749 inline bool operator!=(const srmatrix_subv& v1, const rvector_slice& v2) {
02750   return srvector(v1) != v2;
02751 }
02752 
02754 inline bool operator!=(const srvector& v1, const srmatrix_subv& v2) {
02755   return v1 != srvector(v2);
02756 }
02757 
02759 inline bool operator!=(const srvector_slice& v1, const srmatrix_subv& v2) {
02760   return v1 != srvector(v2);
02761 }
02762 
02764 inline bool operator!=(const rvector& v1, const srmatrix_subv& v2) {
02765   return v1 != srvector(v2);
02766 }
02767 
02769 inline bool operator!=(const rvector_slice& v1, const srmatrix_subv& v2) {
02770   return v1 != srvector(v2);
02771 }
02772 
02774 inline bool operator<(const srmatrix_subv& v1, const srvector& v2) {
02775   return srvector(v1) < v2;
02776 }
02777 
02779 inline bool operator<(const srmatrix_subv& v1, const srvector_slice& v2) {
02780   return srvector(v1) < v2;
02781 }
02782 
02784 inline bool operator<(const srmatrix_subv& v1, const rvector& v2) {
02785   return srvector(v1) < v2;
02786 }
02787 
02789 inline bool operator<(const srmatrix_subv& v1, const rvector_slice& v2) {
02790   return srvector(v1) < v2;
02791 }
02792 
02794 inline bool operator<(const srvector& v1, const srmatrix_subv& v2) {
02795   return v1 < srvector(v2);
02796 }
02797 
02799 inline bool operator<(const srvector_slice& v1, const srmatrix_subv& v2) {
02800   return v1 < srvector(v2);
02801 }
02802 
02804 inline bool operator<(const rvector& v1, const srmatrix_subv& v2) {
02805   return v1 < srvector(v2);
02806 }
02807 
02809 inline bool operator<(const rvector_slice& v1, const srmatrix_subv& v2) {
02810   return v1 < srvector(v2);
02811 }
02812 
02814 inline bool operator<=(const srmatrix_subv& v1, const srvector& v2) {
02815   return srvector(v1) <= v2;
02816 }
02817 
02819 inline bool operator<=(const srmatrix_subv& v1, const srvector_slice& v2) {
02820   return srvector(v1) <= v2;
02821 }
02822 
02824 inline bool operator<=(const srmatrix_subv& v1, const rvector& v2) {
02825   return srvector(v1) <= v2;
02826 }
02827 
02829 inline bool operator<=(const srmatrix_subv& v1, const rvector_slice& v2) {
02830   return srvector(v1) <= v2;
02831 }
02832 
02834 inline bool operator<=(const srvector& v1, const srmatrix_subv& v2) {
02835   return v1 <= srvector(v2);
02836 }
02837 
02839 inline bool operator<=(const srvector_slice& v1, const srmatrix_subv& v2) {
02840   return v1 <= srvector(v2);
02841 }
02842 
02844 inline bool operator<=(const rvector& v1, const srmatrix_subv& v2) {
02845   return v1 <= srvector(v2);
02846 }
02847 
02849 inline bool operator<=(const rvector_slice& v1, const srmatrix_subv& v2) {
02850   return v1 <= srvector(v2);
02851 }
02852 
02854 inline bool operator>(const srmatrix_subv& v1, const srvector& v2) {
02855   return srvector(v1) > v2;
02856 }
02857 
02859 inline bool operator>(const srmatrix_subv& v1, const srvector_slice& v2) {
02860   return srvector(v1) > v2;
02861 }
02862 
02864 inline bool operator>(const srmatrix_subv& v1, const rvector& v2) {
02865   return srvector(v1) > v2;
02866 }
02867 
02869 inline bool operator>(const srmatrix_subv& v1, const rvector_slice& v2) {
02870   return srvector(v1) > v2;
02871 }
02872 
02874 inline bool operator>(const srvector& v1, const srmatrix_subv& v2) {
02875   return v1 > srvector(v2);
02876 }
02877 
02879 inline bool operator>(const srvector_slice& v1, const srmatrix_subv& v2) {
02880   return v1 > srvector(v2);
02881 }
02882 
02884 inline bool operator>(const rvector& v1, const srmatrix_subv& v2) {
02885   return v1 > srvector(v2);
02886 }
02887 
02889 inline bool operator>(const rvector_slice& v1, const srmatrix_subv& v2) {
02890   return v1 > srvector(v2);
02891 }
02892 
02894 inline bool operator>=(const srmatrix_subv& v1, const srvector& v2) {
02895   return srvector(v1) >= v2;
02896 }
02897 
02899 inline bool operator>=(const srmatrix_subv& v1, const srvector_slice& v2) {
02900   return srvector(v1) >= v2;
02901 }
02902 
02904 inline bool operator>=(const srmatrix_subv& v1, const rvector& v2) {
02905   return srvector(v1) >= v2;
02906 }
02907 
02909 inline bool operator>=(const srmatrix_subv& v1, const rvector_slice& v2) {
02910   return srvector(v1) >= v2;
02911 }
02912 
02914 inline bool operator>=(const srvector& v1, const srmatrix_subv& v2) {
02915   return v1 >= srvector(v2);
02916 }
02917 
02919 inline bool operator>=(const srvector_slice& v1, const srmatrix_subv& v2) {
02920   return v1 >= srvector(v2);
02921 }
02922 
02924 inline bool operator>=(const rvector& v1, const srmatrix_subv& v2) {
02925   return v1 >= srvector(v2);
02926 }
02927 
02929 inline bool operator>=(const rvector_slice& v1, const srmatrix_subv& v2) {
02930   return v1 >= srvector(v2);
02931 }
02932 
02934 inline bool operator!(const srmatrix_subv& v) {
02935   return sv_v_not(v);
02936 }
02937 
02939 
02942 inline void accumulate(dotprecision& dot, const srmatrix_subv& v1, const srmatrix_subv& v2) {
02943   spsp_vv_accu<dotprecision,srvector,srvector,sparse_dot>(dot, srvector(v1), srvector(v2));
02944 }
02945 
02947 
02950 inline void accumulate(dotprecision& dot, const srmatrix_subv& v1, const srvector& v2) {
02951   spsp_vv_accu<dotprecision,srvector,srvector,sparse_dot>(dot, srvector(v1), v2);
02952 }
02953 
02955 
02958 inline void accumulate(dotprecision& dot, const srmatrix_subv& v1, const srvector_slice& v2) {
02959   spsl_vv_accu<dotprecision,srvector,srvector_slice,sparse_dot>(dot, srvector(v1), v2);
02960 }
02961 
02963 
02966 inline void accumulate(dotprecision& dot, const srmatrix_subv& v1, const rvector& v2) {
02967   spf_vv_accu<dotprecision,srvector,rvector,sparse_dot>(dot, srvector(v1), v2);
02968 }
02969 
02971 
02974 inline void accumulate(dotprecision& dot, const srmatrix_subv& v1, const rvector_slice& v2) {
02975   spf_vv_accu<dotprecision,srvector,rvector_slice,sparse_dot>(dot, srvector(v1), v2);
02976 }
02977 
02979 
02982 inline void accumulate(dotprecision& dot, const srvector& v1, const srmatrix_subv& v2) {
02983   spsp_vv_accu<dotprecision,srvector,srvector,sparse_dot>(dot, v1, srvector(v2));
02984 }
02985 
02987 
02990 inline void accumulate(dotprecision& dot, const srvector_slice& v1, const srmatrix_subv& v2) {
02991   slsp_vv_accu<dotprecision,srvector_slice,srvector,sparse_dot>(dot, v1, srvector(v2));
02992 }
02993 
02995 
02998 inline void accumulate(dotprecision& dot, const rvector& v1, const srmatrix_subv& v2) {
02999   fsp_vv_accu<dotprecision,rvector,srvector,sparse_dot>(dot, v1, srvector(v2));
03000 }
03001 
03003 
03006 inline void accumulate(dotprecision& dot, const rvector_slice& v1, const srmatrix_subv& v2) {
03007   fsp_vv_accu<dotprecision,rvector_slice,srvector,sparse_dot>(dot, v1, srvector(v2));
03008 }
03009 
03011 
03016 inline void accumulate_approx(dotprecision& dot, const srmatrix_subv& v1, const srmatrix_subv& v2) {
03017   spsp_vv_accuapprox<dotprecision,srvector,srvector,sparse_dot>(dot, srvector(v1), srvector(v2));
03018 }
03019 
03021 
03026 inline void accumulate_approx(dotprecision& dot, const srmatrix_subv& v1, const srvector& v2) {
03027   spsp_vv_accuapprox<dotprecision,srvector,srvector,sparse_dot>(dot, srvector(v1), v2);
03028 }
03029 
03031 
03036 inline void accumulate_approx(dotprecision& dot, const srmatrix_subv& v1, const srvector_slice& v2) {
03037   spsl_vv_accuapprox<dotprecision,srvector,srvector_slice,sparse_dot>(dot, srvector(v1), v2);
03038 }
03039 
03041 
03046 inline void accumulate_approx(dotprecision& dot, const srmatrix_subv& v1, const rvector& v2) {
03047   spf_vv_accuapprox<dotprecision,srvector,rvector,sparse_dot>(dot, srvector(v1), v2);
03048 }
03049 
03051 
03056 inline void accumulate_approx(dotprecision& dot, const srmatrix_subv& v1, const rvector_slice& v2) {
03057   spf_vv_accuapprox<dotprecision,srvector,rvector_slice,sparse_dot>(dot, srvector(v1), v2);
03058 }
03059 
03061 
03066 inline void accumulate_approx(dotprecision& dot, const srvector& v1, const srmatrix_subv& v2) {
03067   spsp_vv_accuapprox<dotprecision,srvector,srvector,sparse_dot>(dot, v1, srvector(v2));
03068 }
03069 
03071 
03076 inline void accumulate_approx(dotprecision& dot, const srvector_slice& v1, const srmatrix_subv& v2) {
03077   slsp_vv_accuapprox<dotprecision,srvector_slice,srvector,sparse_dot>(dot, v1, srvector(v2));
03078 }
03079 
03081 
03086 inline void accumulate_approx(dotprecision& dot, const rvector& v1, const srmatrix_subv& v2) {
03087   fsp_vv_accuapprox<dotprecision,rvector,srvector,sparse_dot>(dot, v1, srvector(v2));
03088 }
03089 
03091 
03096 inline void accumulate_approx(dotprecision& dot, const rvector_slice& v1, const srmatrix_subv& v2) {
03097   fsp_vv_accuapprox<dotprecision,rvector_slice,srvector,sparse_dot>(dot, v1, srvector(v2));
03098 }
03099 
03101 
03104 inline void accumulate(idotprecision& dot, const srmatrix_subv& v1, const srmatrix_subv& v2) {
03105   dotprecision tmp(0.0);
03106   tmp.set_k(dot.get_k());
03107   accumulate(tmp,srvector(v1),srvector(v2));
03108   dot += tmp;
03109 }
03110 
03112 
03115 inline void accumulate(idotprecision& dot, const srmatrix_subv& v1, const srvector& v2) {
03116   dotprecision tmp(0.0);
03117   tmp.set_k(dot.get_k());
03118   accumulate(tmp,srvector(v1),v2);
03119   dot += tmp;
03120 }
03121 
03123 
03126 inline void accumulate(idotprecision& dot, const srmatrix_subv& v1, const srvector_slice& v2) {
03127   dotprecision tmp(0.0);
03128   tmp.set_k(dot.get_k());
03129   accumulate(tmp,srvector(v1),v2);
03130   dot += tmp;
03131 }
03132 
03134 
03137 inline void accumulate(idotprecision& dot, const srmatrix_subv& v1, const rvector& v2) {
03138   dotprecision tmp(0.0);
03139   tmp.set_k(dot.get_k());
03140   accumulate(tmp,srvector(v1),v2);
03141   dot += tmp;
03142 }
03143 
03145 
03148 inline void accumulate(idotprecision& dot, const srmatrix_subv& v1, const rvector_slice& v2) {
03149   dotprecision tmp(0.0);
03150   tmp.set_k(dot.get_k());
03151   accumulate(tmp,srvector(v1),v2);
03152   dot += tmp;
03153 }
03154 
03156 
03159 inline void accumulate(idotprecision& dot, const srvector& v1, const srmatrix_subv& v2) {
03160   dotprecision tmp(0.0);
03161   tmp.set_k(dot.get_k());
03162   accumulate(tmp,v1,srvector(v2));
03163   dot += tmp;
03164 }
03165 
03167 
03170 inline void accumulate(idotprecision& dot, const srvector_slice& v1, const srmatrix_subv& v2) {
03171   dotprecision tmp(0.0);
03172   tmp.set_k(dot.get_k());
03173   accumulate(tmp,v1,srvector(v2));
03174   dot += tmp;
03175 }
03176 
03178 
03181 inline void accumulate(idotprecision& dot, const rvector& v1, const srmatrix_subv& v2) {
03182   dotprecision tmp(0.0);
03183   tmp.set_k(dot.get_k());
03184   accumulate(tmp,v1,srvector(v2));
03185   dot += tmp;
03186 }
03187 
03189 
03192 inline void accumulate(idotprecision& dot, const rvector_slice& v1, const srmatrix_subv& v2) {
03193   dotprecision tmp(0.0);
03194   tmp.set_k(dot.get_k());
03195   accumulate(tmp,v1,srvector(v2));
03196   dot += tmp;
03197 }
03198 
03200 
03203 inline void accumulate(cdotprecision& dot, const srmatrix_subv& v1, const srmatrix_subv& v2) {
03204   dotprecision tmp(0.0);
03205   tmp.set_k(dot.get_k());
03206   accumulate(tmp,srvector(v1),srvector(v2));
03207   dot += tmp;
03208 }
03209 
03211 
03214 inline void accumulate(cdotprecision& dot, const srmatrix_subv& v1, const srvector& v2) {
03215   dotprecision tmp(0.0);
03216   tmp.set_k(dot.get_k());
03217   accumulate(tmp,srvector(v1),v2);
03218   dot += tmp;
03219 }
03220 
03222 
03225 inline void accumulate(cdotprecision& dot, const srmatrix_subv& v1, const srvector_slice& v2) {
03226   dotprecision tmp(0.0);
03227   tmp.set_k(dot.get_k());
03228   accumulate(tmp,srvector(v1),v2);
03229   dot += tmp;
03230 }
03231 
03233 
03236 inline void accumulate(cdotprecision& dot, const srmatrix_subv& v1, const rvector& v2) {
03237   dotprecision tmp(0.0);
03238   tmp.set_k(dot.get_k());
03239   accumulate(tmp,srvector(v1),v2);
03240   dot += tmp;
03241 }
03242 
03244 
03247 inline void accumulate(cdotprecision& dot, const srmatrix_subv& v1, const rvector_slice& v2) {
03248   dotprecision tmp(0.0);
03249   tmp.set_k(dot.get_k());
03250   accumulate(tmp,srvector(v1),v2);
03251   dot += tmp;
03252 }
03253 
03255 
03258 inline void accumulate(cdotprecision& dot, const srvector& v1, const srmatrix_subv& v2) {
03259   dotprecision tmp(0.0);
03260   tmp.set_k(dot.get_k());
03261   accumulate(tmp,v1,srvector(v2));
03262   dot += tmp;
03263 }
03264 
03266 
03269 inline void accumulate(cdotprecision& dot, const srvector_slice& v1, const srmatrix_subv& v2) {
03270   dotprecision tmp(0.0);
03271   tmp.set_k(dot.get_k());
03272   accumulate(tmp,v1,srvector(v2));
03273   dot += tmp;
03274 }
03275 
03277 
03280 inline void accumulate(cdotprecision& dot, const rvector& v1, const srmatrix_subv& v2) {
03281   dotprecision tmp(0.0);
03282   tmp.set_k(dot.get_k());
03283   accumulate(tmp,v1,srvector(v2));
03284   dot += tmp;
03285 }
03286 
03288 
03293 inline void accumulate_approx(cdotprecision& dot, const srmatrix_subv& v1, const srmatrix_subv& v2) {
03294   dotprecision tmp(0.0);
03295   tmp.set_k(dot.get_k());
03296   accumulate_approx(tmp,srvector(v1),srvector(v2));
03297   dot += tmp;
03298 }
03299 
03301 
03306 inline void accumulate_approx(cdotprecision& dot, const srmatrix_subv& v1, const srvector& v2) {
03307   dotprecision tmp(0.0);
03308   tmp.set_k(dot.get_k());
03309   accumulate_approx(tmp,srvector(v1),v2);
03310   dot += tmp;
03311 }
03312 
03314 
03319 inline void accumulate_approx(cdotprecision& dot, const srmatrix_subv& v1, const srvector_slice& v2) {
03320   dotprecision tmp(0.0);
03321   tmp.set_k(dot.get_k());
03322   accumulate_approx(tmp,srvector(v1),v2);
03323   dot += tmp;
03324 }
03325 
03327 
03332 inline void accumulate_approx(cdotprecision& dot, const srmatrix_subv& v1, const rvector& v2) {
03333   dotprecision tmp(0.0);
03334   tmp.set_k(dot.get_k());
03335   accumulate_approx(tmp,srvector(v1),v2);
03336   dot += tmp;
03337 }
03338 
03340 
03345 inline void accumulate_approx(cdotprecision& dot, const srmatrix_subv& v1, const rvector_slice& v2) {
03346   dotprecision tmp(0.0);
03347   tmp.set_k(dot.get_k());
03348   accumulate_approx(tmp,srvector(v1),v2);
03349   dot += tmp;
03350 }
03351 
03353 
03358 inline void accumulate_approx(cdotprecision& dot, const srvector& v1, const srmatrix_subv& v2) {
03359   dotprecision tmp(0.0);
03360   tmp.set_k(dot.get_k());
03361   accumulate_approx(tmp,v1,srvector(v2));
03362   dot += tmp;
03363 }
03364 
03366 
03371 inline void accumulate_approx(cdotprecision& dot, const srvector_slice& v1, const srmatrix_subv& v2) {
03372   dotprecision tmp(0.0);
03373   tmp.set_k(dot.get_k());
03374   accumulate_approx(tmp,v1,srvector(v2));
03375   dot += tmp;
03376 }
03377 
03379 
03384 inline void accumulate_approx(cdotprecision& dot, const rvector& v1, const srmatrix_subv& v2) {
03385   dotprecision tmp(0.0);
03386   tmp.set_k(dot.get_k());
03387   accumulate_approx(tmp,v1,srvector(v2));
03388   dot += tmp;
03389 }
03390 
03392 
03395 inline void accumulate(cdotprecision& dot, const rvector_slice& v1, const srmatrix_subv& v2) {
03396   dotprecision tmp(0.0);
03397   tmp.set_k(dot.get_k());
03398   accumulate(tmp,v1,srvector(v2));
03399   dot += tmp;
03400 }
03401 
03403 
03406 inline void accumulate(cidotprecision& dot, const srmatrix_subv& v1, const srmatrix_subv& v2) {
03407   dotprecision tmp(0.0);
03408   tmp.set_k(dot.get_k());
03409   accumulate(tmp,srvector(v1),srvector(v2));
03410   dot += tmp;
03411 }
03412 
03414 
03417 inline void accumulate(cidotprecision& dot, const srmatrix_subv& v1, const srvector& v2) {
03418   dotprecision tmp(0.0);
03419   tmp.set_k(dot.get_k());
03420   accumulate(tmp,srvector(v1),v2);
03421   dot += tmp;
03422 }
03423 
03425 
03428 inline void accumulate(cidotprecision& dot, const srmatrix_subv& v1, const srvector_slice& v2) {
03429   dotprecision tmp(0.0);
03430   tmp.set_k(dot.get_k());
03431   accumulate(tmp,srvector(v1),v2);
03432   dot += tmp;
03433 }
03434 
03436 
03439 inline void accumulate(cidotprecision& dot, const srmatrix_subv& v1, const rvector& v2) {
03440   dotprecision tmp(0.0);
03441   tmp.set_k(dot.get_k());
03442   accumulate(tmp,srvector(v1),v2);
03443   dot += tmp;
03444 }
03445 
03447 
03450 inline void accumulate(cidotprecision& dot, const srmatrix_subv& v1, const rvector_slice& v2) {
03451   dotprecision tmp(0.0);
03452   tmp.set_k(dot.get_k());
03453   accumulate(tmp,srvector(v1),v2);
03454   dot += tmp;
03455 }
03456 
03458 
03461 inline void accumulate(cidotprecision& dot, const srvector& v1, const srmatrix_subv& v2) {
03462   dotprecision tmp(0.0);
03463   tmp.set_k(dot.get_k());
03464   accumulate(tmp,v1,srvector(v2));
03465   dot += tmp;
03466 }
03467 
03469 
03472 inline void accumulate(cidotprecision& dot, const srvector_slice& v1, const srmatrix_subv& v2) {
03473   dotprecision tmp(0.0);
03474   tmp.set_k(dot.get_k());
03475   accumulate(tmp,v1,srvector(v2));
03476   dot += tmp;
03477 }
03478 
03480 
03483 inline void accumulate(cidotprecision& dot, const rvector& v1, const srmatrix_subv& v2) {
03484   dotprecision tmp(0.0);
03485   tmp.set_k(dot.get_k());
03486   accumulate(tmp,v1,srvector(v2));
03487   dot += tmp;
03488 }
03489 
03491 
03494 inline void accumulate(cidotprecision& dot, const rvector_slice& v1, const srmatrix_subv& v2) {
03495   dotprecision tmp(0.0);
03496   tmp.set_k(dot.get_k());
03497   accumulate(tmp,v1,srvector(v2));
03498   dot += tmp;
03499 }
03500 
03501 }  //namespace cxsc;
03502 
03503 #include "sparsematrix.inl"
03504 
03505 #endif