C-XSC - A C++ Class Library for Extended Scientific Computing
2.5.4
|
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