C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
srmatrix.hpp
1 /*
2 ** CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4)
3 **
4 ** Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik,
5 ** Universitaet Karlsruhe, Germany
6 ** (C) 2000-2014 Wiss. Rechnen/Softwaretechnologie
7 ** Universitaet Wuppertal, Germany
8 **
9 ** This library is free software; you can redistribute it and/or
10 ** modify it under the terms of the GNU Library General Public
11 ** License as published by the Free Software Foundation; either
12 ** version 2 of the License, or (at your option) any later version.
13 **
14 ** This library is distributed in the hope that it will be useful,
15 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
16 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 ** Library General Public License for more details.
18 **
19 ** You should have received a copy of the GNU Library General Public
20 ** License along with this library; if not, write to the Free
21 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 */
23 
24 /* CVS $Id: srmatrix.hpp,v 1.21 2014/01/30 17:23:49 cxsc Exp $ */
25 
26 #ifndef _CXSC_SRMATRIX_HPP_INCLUDED
27 #define _CXSC_SRMATRIX_HPP_INCLUDED
28 
29 #include <real.hpp>
30 #include <rmatrix.hpp>
31 #include <srvector.hpp>
32 #include <cidot.hpp>
33 #include <vector>
34 #include <algorithm>
35 #include <iostream>
36 #include <sparsedot.hpp>
37 #include <sparsematrix.hpp>
38 
39 namespace cxsc {
40 
42 enum STORAGE_TYPE{triplet,compressed_row,compressed_column};
43 
44 class srmatrix_slice;
45 class srmatrix_subv;
46 class scmatrix;
47 class scmatrix_slice;
48 class scmatrix_subv;
49 class simatrix;
50 class simatrix_slice;
51 class simatrix_subv;
52 class sivector_slice;
53 class scimatrix;
54 class scimatrix_slice;
55 class scimatrix_subv;
56 class scivector_slice;
57 
58 
59 inline bool comp_pair(std::pair<int,real> p1, std::pair<int,real> p2) {
60  return p1.first < p2.first;
61 }
62 
63 
65 
77 class srmatrix {
78 
79  private:
80  std::vector<int> p;
81  std::vector<int> ind;
82  std::vector<real> x;
83  int m;
84  int n;
85  int lb1,ub1,lb2,ub2;
86 
87  public:
88 
90  std::vector<int>& column_pointers() {
91  return p;
92  }
93 
95  std::vector<int>& row_indices() {
96  return ind;
97  }
98 
100  std::vector<real>& values() {
101  return x;
102  }
103 
105  const std::vector<int>& column_pointers() const {
106  return p;
107  }
108 
110  const std::vector<int>& row_indices() const {
111  return ind;
112  }
113 
115  const std::vector<real>& values() const {
116  return x;
117  }
118 
121  p.push_back(0);
122  m = n = 0;
123  lb1 = lb2 = ub1 = ub2 = 0;
124  }
125 
127  srmatrix(const int r, const int c) : m(r),n(c),lb1(1),ub1(r),lb2(1),ub2(c) {
128  p = std::vector<int>((n>0) ? n+1 : 1, 0);
129  ind.reserve(2*(m+n));
130  x.reserve(2*(m+n));
131 
132  p[0] = 0;
133  }
134 
136  srmatrix(const int r, const int c, const int e) : m(r),n(c),lb1(1),ub1(r),lb2(1),ub2(c) {
137  p = std::vector<int>((n>0) ? n+1 : 1, 0);
138  ind.reserve(e);
139  x.reserve(e);
140 
141  p[0] = 0;
142  }
143 
145 
151  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) {
152  if(t == triplet) {
153  this->m = m;
154  this->n = n;
155  p = std::vector<int>(n+1,0);
156  ind.reserve(nnz);
157  x.reserve(nnz);
158  lb1 = lb2 = 1;
159  ub1 = m; ub2 = n;
160 
161  std::vector<triplet_store<real> > work;
162  work.reserve(nnz);
163 
164  for(int k=0 ; k<nnz ; k++) {
165  work.push_back(triplet_store<real>(rows[Lb(rows)+k],cols[Lb(cols)+k],values[Lb(values)+k]));
166  }
167 
168  sort(work.begin(), work.end());
169 
170  int i=0;
171 
172  for(int j=0 ; j<n ; j++) {
173 
174  while((unsigned int)i < work.size() && work[i].col == j ) {
175  ind.push_back(work[i].row);
176  x.push_back(work[i].val);
177  i++;
178  }
179 
180  p[j+1] = i;
181  }
182 
183  } else if(t == compressed_row) {
184 
185  this->m = m;
186  this->n = n;
187  p = std::vector<int>(n+1,0);
188  ind.reserve(nnz);
189  x.reserve(nnz);
190  lb1 = lb2 = 1;
191  ub1 = m; ub2 = n;
192 
193  for(int i=0 ; i<n+1 ; i++)
194  p[i] = rows[Lb(rows)+i];
195 
196  std::vector<triplet_store<real> > work;
197  work.reserve(nnz);
198 
199  for(int j=0 ; j<n ; j++) {
200  for(int k=p[j] ; k<p[j+1] ; k++) {
201  work.push_back(triplet_store<real>(j,cols[Lb(cols)+k],values[Lb(values)+k]));
202  }
203  }
204 
205  sort(work.begin(), work.end());
206 
207  int i=0;
208 
209  for(int j=0 ; j<n ; j++) {
210 
211  while((unsigned int)i < work.size() && work[i].col == j ) {
212  ind.push_back(work[i].row);
213  x.push_back(work[i].val);
214  i++;
215  }
216 
217  p[j+1] = i;
218  }
219 
220  } else if(t == compressed_column) {
221  this->m = m;
222  this->n = n;
223  p = std::vector<int>(n+1,0);
224  ind.reserve(nnz);
225  x.reserve(nnz);
226  lb1 = lb2 = 1;
227  ub1 = m; ub2 = n;
228 
229  for(int i=0 ; i<n+1 ; i++)
230  p[i] = rows[Lb(rows)+i];
231 
232  std::vector<std::pair<int,real> > work;
233  work.reserve(n);
234 
235  for(int j=0 ; j<n ; j++) {
236  work.clear();
237 
238  for(int k=p[j] ; k<p[j+1] ; k++) {
239  work.push_back(std::make_pair(cols[Lb(cols)+k],values[Lb(values)+k]));
240  }
241 
242  std::sort(work.begin(),work.end(),comp_pair);
243 
244  for(unsigned int i=0 ; i<work.size() ; i++) {
245  ind.push_back(work[i].first);
246  x.push_back(work[i].second);
247  }
248  }
249 
250  }
251 
252  }
253 
254 
256 
263  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) {
264  if(t == triplet) {
265  this->m = m;
266  this->n = n;
267  p = std::vector<int>(n+1,0);
268  ind.reserve(nnz);
269  x.reserve(nnz);
270  lb1 = lb2 = 1;
271  ub1 = m; ub2 = n;
272 
273  std::vector<triplet_store<real> > work;
274  work.reserve(nnz);
275 
276  for(int k=0 ; k<nnz ; k++) {
277  work.push_back(triplet_store<real>(rows[k],cols[k],values[k]));
278  }
279 
280  sort(work.begin(), work.end());
281 
282  int i=0;
283 
284  for(int j=0 ; j<n ; j++) {
285 
286  while((unsigned int)i < work.size() && work[i].col == j ) {
287  ind.push_back(work[i].row);
288  x.push_back(work[i].val);
289  i++;
290  }
291 
292  p[j+1] = i;
293  }
294 
295  } else if(t == compressed_row) {
296 
297  this->m = m;
298  this->n = n;
299  p = std::vector<int>(n+1,0);
300  ind.reserve(nnz);
301  x.reserve(nnz);
302  lb1 = lb2 = 1;
303  ub1 = m; ub2 = n;
304 
305  for(int i=0 ; i<n+1 ; i++)
306  p[i] = rows[i];
307 
308  std::vector<triplet_store<real> > work;
309  work.reserve(nnz);
310 
311  for(int j=0 ; j<n ; j++) {
312  for(int k=p[j] ; k<p[j+1] ; k++) {
313  work.push_back(triplet_store<real>(j,cols[k],values[k]));
314  }
315  }
316 
317  sort(work.begin(), work.end());
318 
319  int i=0;
320 
321  for(int j=0 ; j<n ; j++) {
322 
323  while((unsigned int)i < work.size() && work[i].col == j ) {
324  ind.push_back(work[i].row);
325  x.push_back(work[i].val);
326  i++;
327  }
328 
329  p[j+1] = i;
330  }
331 
332  } else if(t == compressed_column) {
333  this->m = m;
334  this->n = n;
335  p = std::vector<int>(n+1,0);
336  ind.reserve(nnz);
337  x.reserve(nnz);
338  lb1 = lb2 = 1;
339  ub1 = m; ub2 = n;
340 
341  for(int i=0 ; i<n+1 ; i++)
342  p[i] = rows[i];
343 
344  std::vector<std::pair<int,real> > work;
345  work.reserve(n);
346 
347  for(int j=0 ; j<n ; j++) {
348  work.clear();
349 
350  for(int k=p[j] ; k<p[j+1] ; k++) {
351  work.push_back(std::make_pair(cols[k],values[k]));
352  }
353 
354  std::sort(work.begin(),work.end(),comp_pair);
355 
356  for(unsigned int i=0 ; i<work.size() ; i++) {
357  ind.push_back(work[i].first);
358  x.push_back(work[i].second);
359  }
360  }
361 
362  }
363 
364  }
365 
367  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)) {
368  p = std::vector<int>((n>0) ? n+1 : 1, 0);
369  ind.reserve((m*n*0.1 < 2*m) ? (int)(m*n*0.1) : 2*m);
370  x.reserve((m*n*0.1 < 2*m) ? (int)(m*n*0.1) : 2*m);
371 
372  p[0] = 0;
373  int nnz = 0;
374 
375  for(int j=0 ; j<n ; j++) {
376  for(int i=0 ; i<m ; i++) {
377  if(A[i+lb1][j+lb2] != 0.0) {
378  ind.push_back(i);
379  x.push_back(A[i+lb1][j+lb2]);
380  nnz++;
381  }
382  }
383 
384  p[j+1] = nnz;
385  }
386 
387  }
388 
389 
391 
394  srmatrix(const int ms, const int ns, const rmatrix& A) : m(ms), n(ns), lb1(1), ub1(ms), lb2(1), ub2(ns) {
395  //Banded matrix constructor
396  int nnz = RowLen(A)*ColLen(A);
397  p = std::vector<int>((n>0) ? n+1 : 1, 0);
398  ind.reserve(nnz);
399  x.reserve(nnz);
400 
401  std::vector<triplet_store<real> > work;
402  work.reserve(nnz);
403 
404 
405  for(int i=0 ; i<ColLen(A) ; i++) {
406  for(int j=Lb(A,2) ; j<=Ub(A,2) ; j++) {
407  if(i+j >=0 && i+j < n) {
408  work.push_back(triplet_store<real>(i,i+j,A[i+Lb(A,1)][j]));
409  }
410  }
411  }
412 
413  sort(work.begin(), work.end());
414 
415  int i=0;
416 
417  for(int j=0 ; j<n ; j++) {
418 
419  while((unsigned int)i < work.size() && work[i].col == j ) {
420  ind.push_back(work[i].row);
421  x.push_back(work[i].val);
422  i++;
423  }
424 
425  p[j+1] = i;
426  }
427 
428  }
429 
431  srmatrix(const srmatrix_slice&);
432 
434  void full(rmatrix& A) const {
435  A = rmatrix(lb1,ub1,lb2,ub2);
436  A = 0.0;
437  for(int j=0 ; j<n ; j++) {
438  for(int k=p[j] ; k<p[j+1] ; k++) {
439  A[ind[k]+lb1][j+lb2] = x[k];
440  }
441  }
442  }
443 
445 
449  void dropzeros() {
450  std::vector<int> pnew(n+1,0);
451  std::vector<int> indnew;
452  std::vector<real> xnew;
453  int nnznew = 0;
454 
455  for(int j=0 ; j<n ; j++) {
456  for(int k=p[j] ; k<p[j+1] ; k++) {
457  if(x[k] != 0.0) {
458  xnew.push_back(x[k]);
459  indnew.push_back(ind[k]);
460  nnznew++;
461  }
462  }
463  pnew[j+1] = nnznew;
464  }
465 
466  p = pnew;
467  ind = indnew;
468  x = xnew;
469  }
470 
471 
473  srmatrix& operator=(const real& A) {
474  return sp_ms_assign<srmatrix,real,real>(*this,A);
475  }
476 
479  return spf_mm_assign<srmatrix,rmatrix,real>(*this,A);
480  }
481 
484  return spf_mm_assign<srmatrix,rmatrix,real>(*this,A);
485  }
486 
487 /* srmatrix& operator=(const srmatrix& A) {
488  p = A.p;
489  ind = A.ind;
490  x = A.x;
491  return *this;
492  }*/
493 
496 
498 
504  const real operator()(int i, int j) const {
505 #if(CXSC_INDEX_CHECK)
506  if(i<lb1 || i>ub1 || j<lb2 || j>ub2)
507  cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix::operator()(int, int)"));
508 #endif
509  real r = 0.0;
510  for(int k=p[j-lb2] ; k<p[j-lb2+1] && ind[k]<=i-lb1 ; k++) {
511  if(ind[k] == i-lb1) r = x[k];
512  }
513  return r;
514  }
515 
517 
525  real& element(int i, int j) {
526 #if(CXSC_INDEX_CHECK)
527  if(i<lb1 || i>ub1 || j<lb2 || j>ub2)
528  cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix::element(int, int)"));
529 #endif
530  int k;
531  for(k=p[j-lb2] ; k<p[j-lb2+1] && ind[k]<=i-lb1 ; k++) {
532  if(ind[k] == i-lb1) return x[k];
533  }
534 
535  //Nicht gefunden, Element muss angelegt werden, da Schreibzugriff moeglich
536  std::vector<int>::iterator ind_it = ind.begin() + k;
537  std::vector<real>::iterator x_it = x.begin() + k;
538  ind.insert(ind_it, i-lb1);
539  x_it = x.insert(x_it, 0.0);
540  for(k=j-lb2+1 ; k<(int)p.size() ; k++)
541  p[k]++;
542 
543  return *x_it;
544  }
545 
547  srmatrix_subv operator[](const cxscmatrix_column&);
549  srmatrix_subv operator[](const int);
551  const srmatrix_subv operator[](const cxscmatrix_column&) const;
553  const srmatrix_subv operator[](const int) const;
554 
556  srmatrix_slice operator()(const int, const int , const int, const int);
558  const srmatrix_slice operator()(const int, const int , const int, const int) const;
559 
561  srmatrix operator()(const intvector& pervec, const intvector& q) {
562  srmatrix A(m,n,get_nnz());
563  intvector per = perminv(pervec);
564 
565  int nnz=0;
566  for(int k=0 ; k<n ; k++) {
567  A.p[k] = nnz;
568 
569  std::map<int,real> work;
570  for(int j=p[q[Lb(q)+k]] ; j<p[q[Lb(q)+k]+1] ; j++)
571  work.insert(std::make_pair(per[Lb(per)+ind[j]], x[j]));
572 
573  for(std::map<int,real>::iterator it = work.begin() ; it != work.end() ; it++) {
574  A.ind.push_back(it->first);
575  A.x.push_back(it->second);
576  }
577 
578  nnz += work.size();
579 
580  }
581 
582  A.p[n] = nnz;
583 
584  return A;
585  }
586 
588  srmatrix operator()(const intvector& pervec) {
589  srmatrix A(m,n,get_nnz());
590  intvector per = perminv(pervec);
591 
592  for(int k=0 ; k<n ; k++) {
593  A.p[k] = p[k];
594 
595  std::map<int,real> work;
596  for(int j=p[k] ; j<p[k+1] ; j++)
597  work.insert(std::make_pair(per[Lb(per)+ind[j]], x[j]));
598 
599  for(std::map<int,real>::iterator it = work.begin() ; it != work.end() ; it++) {
600  A.ind.push_back(it->first);
601  A.x.push_back(it->second);
602  }
603 
604  }
605 
606  A.p[n] = p[n];
607 
608  return A;
609  }
610 
612  srmatrix operator()(const intmatrix& P, const intmatrix& Q) {
613  intvector p = permvec(P);
614  intvector q = perminv(permvec(Q));
615  return (*this)(p,q);
616  }
617 
620  intvector p = permvec(P);
621  return (*this)(p);
622  }
623 
625  real density() const {
626  return p[n]/((double)m*n);
627  }
628 
630  int get_nnz() const {
631  return p[n];
632  }
633 
636  return spf_mm_addassign<srmatrix,rmatrix,rmatrix>(*this,B);
637  }
638 
641  return spf_mm_addassign<srmatrix,rmatrix_slice,rmatrix>(*this,B);
642  }
643 
646  return spsp_mm_addassign<srmatrix,srmatrix,real>(*this,B);
647  }
648 
651  return spf_mm_subassign<srmatrix,rmatrix,rmatrix>(*this,B);
652  }
653 
656  return spf_mm_subassign<srmatrix,rmatrix_slice,rmatrix>(*this,B);
657  }
658 
661  return spsp_mm_subassign<srmatrix,srmatrix,real>(*this,B);
662  }
663 
666  return spf_mm_multassign<srmatrix,rmatrix,sparse_dot,rmatrix>(*this,B);
667  }
668 
671  return spf_mm_multassign<srmatrix,rmatrix,sparse_dot,rmatrix>(*this,B);
672  }
673 
676  return spsp_mm_multassign<srmatrix,srmatrix,sparse_dot,real>(*this,B);
677  }
678 
680  srmatrix& operator*=(const real& r) {
681  return sp_ms_multassign(*this,r);
682  }
683 
685  srmatrix& operator/=(const real& r) {
686  return sp_ms_divassign(*this,r);
687  }
688 
689  friend int Lb(const srmatrix&, int);
690  friend int Ub(const srmatrix&, int);
691  friend void SetLb(srmatrix&, const int, const int);
692  friend void SetUb(srmatrix&, const int, const int);
693  friend int RowLen(const srmatrix&);
694  friend int ColLen(const srmatrix&);
695  friend srmatrix Re(const scmatrix&);
696  friend srmatrix Im(const scmatrix&);
697  friend srmatrix Inf(const simatrix&);
698  friend srmatrix Sup(const simatrix&);
699  friend srmatrix InfRe(const scimatrix&);
700  friend srmatrix SupRe(const scimatrix&);
701  friend srmatrix InfIm(const scimatrix&);
702  friend srmatrix SupIm(const scimatrix&);
703  friend srmatrix mid(const simatrix&);
704  friend srmatrix diam(const simatrix&);
705  friend srmatrix absmin(const simatrix&);
706  friend srmatrix absmax(const simatrix&);
707  friend srmatrix abs(const srmatrix&);
708  friend srmatrix abs(const scmatrix&);
709 
710  friend srmatrix CompMat(const simatrix&);
711  friend srmatrix CompMat(const srmatrix&);
712  friend srmatrix CompMat(const scmatrix&);
713  friend srmatrix CompMat(const scimatrix&);
714  friend srmatrix transp(const srmatrix&);
715  friend srmatrix Id(const srmatrix&);
716 
717  friend std::istream& operator>>(std::istream&, srmatrix_slice&);
718  friend std::istream& operator>>(std::istream&, srmatrix_subv&);
719 
720  friend class srmatrix_slice;
721  friend class srmatrix_subv;
722  friend class srvector;
723  friend class scmatrix;
724  friend class scmatrix_slice;
725  friend class scmatrix_subv;
726  friend class simatrix;
727  friend class simatrix_slice;
728  friend class simatrix_subv;
729  friend class scivector;
730  friend class scimatrix;
731  friend class scimatrix_slice;
732  friend class scimatrix_subv;
733  friend class rmatrix;
734  friend class imatrix;
735  friend class cmatrix;
736  friend class cimatrix;
737 
738 #include "matrix_friend_declarations.inl"
739 };
740 
741 inline rmatrix::rmatrix(const srmatrix& A) {
742  dat = new real[A.m*A.n];
743  lb1 = A.lb1; lb2 = A.lb2; ub1 = A.ub1; ub2 = A.ub2;
744  xsize = A.n;
745  ysize = A.m;
746  *this = 0.0;
747 
748  for(int j=0 ; j<A.n ; j++) {
749  for(int k=A.p[j] ; k<A.p[j+1] ; k++) {
750  dat[A.ind[k]*A.n+j] = A.x[k];
751  }
752  }
753 }
754 
756 inline srmatrix Id(const srmatrix& A) {
757  srmatrix I(A.m, A.n, (A.m>A.n) ? A.m : A.n);
758  I.lb1 = A.lb1; I.lb2 = A.lb2;
759  I.ub1 = A.ub1; I.ub2 = A.ub2;
760 
761  if(A.m < A.n) {
762  for(int i=0 ; i<A.m ; i++) {
763  I.p[i+1] = I.p[i] + 1;
764  I.ind.push_back(i);
765  I.x.push_back(1.0);
766  }
767  } else {
768  for(int i=0 ; i<A.n ; i++) {
769  I.p[i+1] = I.p[i] + 1;
770  I.ind.push_back(i);
771  I.x.push_back(1.0);
772  }
773  }
774 
775  return I;
776 }
777 
779 inline srmatrix transp(const srmatrix& A) {
780  srmatrix B(A.n, A.m, A.get_nnz());
781 
782  //NIchtnullen pro Zeile bestimmen
783  std::vector<int> w(A.m,0);
784  for(unsigned int i=0 ; i<A.ind.size() ; i++)
785  w[A.ind[i]]++;
786 
787  //Spalten"pointer" setzen
788  B.p.resize(A.m+1);
789  B.p[0] = 0;
790  for(unsigned int i=1 ; i<B.p.size() ; i++)
791  B.p[i] = w[i-1] + B.p[i-1];
792 
793  //w vorbereiten
794  w.insert(w.begin(), 0);
795  for(unsigned int i=1 ; i<w.size() ; i++) {
796  w[i] += w[i-1];
797  }
798 
799  //neuer zeilenindex und wert wird gesetzt
800  int q;
801  B.ind.resize(A.get_nnz());
802  B.x.resize(A.get_nnz());
803  for(int j=0 ; j<A.n ; j++) {
804  for(int k=A.p[j] ; k<A.p[j+1] ; k++) {
805  q = w[A.ind[k]]++;
806  B.ind[q] = j;
807  B.x[q] = A.x[k];
808  }
809  }
810 
811  return B;
812 }
813 
815 inline srmatrix abs(const srmatrix& A) {
816  srmatrix ret(A);
817  for(unsigned int i=0 ; i<ret.x.size() ; i++)
818  ret.x[i] = abs(ret.x[i]);
819  return ret;
820 }
821 
823 inline srmatrix CompMat(const srmatrix& A) {
824  srmatrix res(A.m,A.n,A.get_nnz());
825  res.lb1 = A.lb1;
826  res.lb2 = A.lb2;
827  res.ub1 = A.ub1;
828  res.ub2 = A.ub2;
829  res.p = A.p;
830  res.ind = A.ind;
831  res.p[A.n] = A.p[A.n];
832 
833  for(int j=0 ; j<res.n ; j++) {
834  for(int k=A.p[j] ; k<A.p[j+1] ; k++) {
835  if(A.ind[k] == j)
836  res.x.push_back(abs(A.x[k]));
837  else
838  res.x.push_back(-abs(A.x[k]));
839  }
840  }
841 
842  res.dropzeros();
843 
844  return res;
845 }
846 
848 
852 inline void SetLb(srmatrix& A, const int i, const int j) {
853  if(i==1) {
854  A.lb1 = j;
855  A.ub1 = j + A.m - 1;
856  } else if(i==2) {
857  A.lb2 = j;
858  A.ub2 = j + A.n - 1;
859  }
860 }
861 
863 
867 inline void SetUb(srmatrix& A, const int i, const int j) {
868  if(i==1) {
869  A.ub1 = j;
870  A.lb1 = j - A.m + 1;
871  } else if(i==2) {
872  A.ub2 = j;
873  A.lb2 = j - A.n + 1;
874  }
875 }
876 
878 
881 inline int Lb(const srmatrix& A, int i) {
882  if(i==1)
883  return A.lb1;
884  else if(i==2)
885  return A.lb2;
886  else
887  return 1;
888 }
889 
891 
894 inline int Ub(const srmatrix& A, int i) {
895  if(i==1)
896  return A.ub1;
897  else if(i==2)
898  return A.ub2;
899  else
900  return 1;
901 }
902 
904 inline int RowLen(const srmatrix& A) {
905  return A.n;
906 }
907 
909 inline int ColLen(const srmatrix& A) {
910  return A.m;
911 }
912 
914 inline void Resize(srmatrix& A) {
915  sp_m_resize(A);
916 }
917 
919 inline void Resize(srmatrix& A, const int m, const int n) {
920  sp_m_resize(A,m,n);
921 }
922 
924 inline void Resize(srmatrix& A, const int l1, const int u1, const int l2, const int u2) {
925  sp_m_resize(A,l1,u1,l2,u2);
926 }
927 
929 
935 inline rmatrix operator*(const rmatrix& A, const srmatrix& B) {
936  return fsp_mm_mult<rmatrix,srmatrix,rmatrix,sparse_dot>(A,B);
937 }
938 
940 
946 inline rmatrix operator*(const srmatrix& A, const rmatrix& B) {
947  return spf_mm_mult<srmatrix,rmatrix,rmatrix,sparse_dot>(A,B);
948 }
949 
951 
957 inline rmatrix operator*(const rmatrix_slice& A, const srmatrix& B) {
958  return fsp_mm_mult<rmatrix_slice,srmatrix,rmatrix,sparse_dot>(A,B);
959 }
960 
962 
968 inline rmatrix operator*(const srmatrix& A, const rmatrix_slice& B) {
969  return spf_mm_mult<srmatrix,rmatrix_slice,rmatrix,sparse_dot>(A,B);
970 }
971 
973 
979 inline srmatrix operator*(const srmatrix& A, const srmatrix& B) {
980  return spsp_mm_mult<srmatrix,srmatrix,srmatrix,sparse_dot,real>(A,B);
981 }
982 
984 inline srmatrix operator*(const srmatrix& A, const real& r) {
985  return sp_ms_mult<srmatrix,real,srmatrix>(A,r);
986 }
987 
989 inline srmatrix operator/(const srmatrix& A, const real& r) {
990  return sp_ms_div<srmatrix,real,srmatrix>(A,r);
991 }
992 
994 inline srmatrix operator*(const real& r, const srmatrix& A) {
995  return sp_sm_mult<real,srmatrix,srmatrix>(r,A);
996 }
997 
999 
1005 inline rvector operator*(const srmatrix& A, const rvector& v) {
1006  return spf_mv_mult<srmatrix,rvector,rvector,sparse_dot>(A,v);
1007 }
1008 
1010 
1016 inline rvector operator*(const srmatrix& A, const rvector_slice& v) {
1017  return spf_mv_mult<srmatrix,rvector_slice,rvector,sparse_dot>(A,v);
1018 }
1019 
1021 
1027 inline srvector operator*(const srmatrix& A, const srvector& v) {
1028  return spsp_mv_mult<srmatrix,srvector,srvector,sparse_dot,real>(A,v);
1029 }
1030 
1032 
1038 inline srvector operator*(const srmatrix& A, const srvector_slice& v) {
1039  return spsl_mv_mult<srmatrix,srvector_slice,srvector,sparse_dot,real>(A,v);
1040 }
1041 
1043 
1049 inline rvector operator*(const rmatrix& A, const srvector& v) {
1050  return fsp_mv_mult<rmatrix,srvector,rvector,sparse_dot>(A,v);
1051 }
1052 
1054 
1060 inline rvector operator*(const rmatrix_slice& A, const srvector& v) {
1061  return fsp_mv_mult<rmatrix_slice,srvector,rvector,sparse_dot>(A,v);
1062 }
1063 
1065 
1071 inline rvector operator*(const rmatrix& A, const srvector_slice& v) {
1072  return fsl_mv_mult<rmatrix,srvector_slice,rvector,sparse_dot>(A,v);
1073 }
1074 
1076 
1082 inline rvector operator*(const rmatrix_slice& A, const srvector_slice& v) {
1083  return fsl_mv_mult<rmatrix_slice,srvector_slice,rvector,sparse_dot>(A,v);
1084 }
1085 
1087 inline rmatrix operator+(const rmatrix& A, const srmatrix& B) {
1088  return fsp_mm_add<rmatrix,srmatrix,rmatrix>(A,B);
1089 }
1090 
1092 inline rmatrix operator+(const srmatrix& A, const rmatrix& B) {
1093  return spf_mm_add<srmatrix,rmatrix,rmatrix>(A,B);
1094 }
1095 
1097 inline rmatrix operator+(const rmatrix_slice& A, const srmatrix& B) {
1098  return fsp_mm_add<rmatrix_slice,srmatrix,rmatrix>(A,B);
1099 }
1100 
1102 inline rmatrix operator+(const srmatrix& A, const rmatrix_slice& B) {
1103  return spf_mm_add<srmatrix,rmatrix_slice,rmatrix>(A,B);
1104 }
1105 
1107 inline srmatrix operator+(const srmatrix& A, const srmatrix& B) {
1108  return spsp_mm_add<srmatrix,srmatrix,srmatrix,real>(A,B);
1109 }
1110 
1112 inline rmatrix operator-(const rmatrix& A, const srmatrix& B) {
1113  return fsp_mm_sub<rmatrix,srmatrix,rmatrix>(A,B);
1114 }
1115 
1117 inline rmatrix operator-(const srmatrix& A, const rmatrix& B) {
1118  return spf_mm_sub<srmatrix,rmatrix,rmatrix>(A,B);
1119 }
1120 
1122 inline rmatrix operator-(const rmatrix_slice& A, const srmatrix& B) {
1123  return fsp_mm_sub<rmatrix_slice,srmatrix,rmatrix>(A,B);
1124 }
1125 
1127 inline rmatrix operator-(const srmatrix& A, const rmatrix_slice& B) {
1128  return spf_mm_sub<srmatrix,rmatrix_slice,rmatrix>(A,B);
1129 }
1130 
1132 inline srmatrix operator-(const srmatrix& A, const srmatrix& B) {
1133  return spsp_mm_sub<srmatrix,srmatrix,srmatrix,real>(A,B);
1134 }
1135 
1137 inline srmatrix operator-(const srmatrix& M) {
1138  return sp_m_negative<srmatrix,srmatrix>(M);
1139 }
1140 
1142 inline srmatrix& operator+(srmatrix& A) {
1143  return A;
1144 }
1145 
1147  *this = rmatrix(B);
1148  return *this;
1149 }
1150 
1152  *this = rmatrix(B);
1153  return *this;
1154 }
1155 
1157  return fsp_mm_addassign(*this,B);
1158 }
1159 
1161  return fsp_mm_addassign(*this,B);
1162 }
1163 
1165  return fsp_mm_subassign(*this,B);
1166 }
1167 
1169  return fsp_mm_subassign(*this,B);
1170 }
1171 
1173  return fsp_mm_multassign<rmatrix,srmatrix,sparse_dot,rmatrix>(*this,B);
1174 }
1175 
1177  return fsp_mm_multassign<rmatrix_slice,srmatrix,sparse_dot,rmatrix>(*this,B);
1178 }
1179 
1181 inline bool operator==(const srmatrix& A, const srmatrix& B) {
1182  return spsp_mm_comp(A,B);
1183 }
1184 
1186 inline bool operator==(const srmatrix& A, const rmatrix& B) {
1187  return spf_mm_comp(A,B);
1188 }
1189 
1191 inline bool operator==(const rmatrix& A, const srmatrix& B) {
1192  return fsp_mm_comp(A,B);
1193 }
1194 
1196 inline bool operator==(const rmatrix_slice& A, const srmatrix& B) {
1197  return fsp_mm_comp(A,B);
1198 }
1199 
1201 inline bool operator==(const srmatrix& A, const rmatrix_slice& B) {
1202  return spf_mm_comp(A,B);
1203 }
1204 
1206 inline bool operator!=(const srmatrix& A, const srmatrix& B) {
1207  return !spsp_mm_comp(A,B);
1208 }
1209 
1211 inline bool operator!=(const srmatrix& A, const rmatrix& B) {
1212  return !spf_mm_comp(A,B);
1213 }
1214 
1216 inline bool operator!=(const rmatrix& A, const srmatrix& B) {
1217  return !fsp_mm_comp(A,B);
1218 }
1219 
1221 inline bool operator!=(const rmatrix_slice& A, const srmatrix& B) {
1222  return !fsp_mm_comp(A,B);
1223 }
1224 
1226 inline bool operator!=(const srmatrix& A, const rmatrix_slice& B) {
1227  return !spf_mm_comp(A,B);
1228 }
1229 
1231 inline bool operator<(const srmatrix& A, const srmatrix& B) {
1232  return spsp_mm_less<srmatrix,srmatrix,real>(A,B);
1233 }
1234 
1236 inline bool operator<(const srmatrix& A, const rmatrix& B) {
1237  return spf_mm_less<srmatrix,rmatrix,real>(A,B);
1238 }
1239 
1241 inline bool operator<(const rmatrix& A, const srmatrix& B) {
1242  return fsp_mm_less<rmatrix,srmatrix,real>(A,B);
1243 }
1244 
1246 inline bool operator<(const rmatrix_slice& A, const srmatrix& B) {
1247  return fsp_mm_less<rmatrix_slice,srmatrix,real>(A,B);
1248 }
1249 
1251 inline bool operator<(const srmatrix& A, const rmatrix_slice& B) {
1252  return spf_mm_less<srmatrix,rmatrix_slice,real>(A,B);
1253 }
1254 
1256 inline bool operator<=(const srmatrix& A, const srmatrix& B) {
1257  return spsp_mm_leq<srmatrix,srmatrix,real>(A,B);
1258 }
1259 
1261 inline bool operator<=(const srmatrix& A, const rmatrix& B) {
1262  return spf_mm_leq<srmatrix,rmatrix,real>(A,B);
1263 }
1264 
1266 inline bool operator<=(const rmatrix& A, const srmatrix& B) {
1267  return fsp_mm_leq<rmatrix,srmatrix,real>(A,B);
1268 }
1269 
1271 inline bool operator<=(const rmatrix_slice& A, const srmatrix& B) {
1272  return fsp_mm_leq<rmatrix_slice,srmatrix,real>(A,B);
1273 }
1274 
1276 inline bool operator<=(const srmatrix& A, const rmatrix_slice& B) {
1277  return spf_mm_leq<srmatrix,rmatrix_slice,real>(A,B);
1278 }
1279 
1281 inline bool operator>(const srmatrix& A, const srmatrix& B) {
1282  return spsp_mm_greater<srmatrix,srmatrix,real>(A,B);
1283 }
1284 
1286 inline bool operator>(const srmatrix& A, const rmatrix& B) {
1287  return spf_mm_greater<srmatrix,rmatrix,real>(A,B);
1288 }
1289 
1291 inline bool operator>(const rmatrix& A, const srmatrix& B) {
1292  return fsp_mm_greater<rmatrix,srmatrix,real>(A,B);
1293 }
1294 
1296 inline bool operator>(const rmatrix_slice& A, const srmatrix& B) {
1297  return fsp_mm_greater<rmatrix_slice,srmatrix,real>(A,B);
1298 }
1299 
1301 inline bool operator>(const srmatrix& A, const rmatrix_slice& B) {
1302  return spf_mm_greater<srmatrix,rmatrix_slice,real>(A,B);
1303 }
1304 
1306 inline bool operator>=(const srmatrix& A, const srmatrix& B) {
1307  return spsp_mm_geq<srmatrix,srmatrix,real>(A,B);
1308 }
1309 
1311 inline bool operator>=(const srmatrix& A, const rmatrix& B) {
1312  return spf_mm_geq<srmatrix,rmatrix,real>(A,B);
1313 }
1314 
1316 inline bool operator>=(const rmatrix& A, const srmatrix& B) {
1317  return fsp_mm_geq<rmatrix,srmatrix,real>(A,B);
1318 }
1319 
1321 inline bool operator>=(const rmatrix_slice& A, const srmatrix& B) {
1322  return fsp_mm_geq<rmatrix_slice,srmatrix,real>(A,B);
1323 }
1324 
1326 inline bool operator>=(const srmatrix& A, const rmatrix_slice& B) {
1327  return spf_mm_geq<srmatrix,rmatrix_slice,real>(A,B);
1328 }
1329 
1331 inline bool operator!(const srmatrix& A) {
1332  return sp_m_not(A);
1333 }
1334 
1336 
1341 inline std::ostream& operator<<(std::ostream& os, const srmatrix& A) {
1342  return sp_m_output<srmatrix,real>(os,A);
1343 }
1344 
1346 
1351 inline std::istream& operator>>(std::istream& is, srmatrix& A) {
1352  return sp_m_input<srmatrix,real>(is,A);
1353 }
1354 
1356 
1361  public:
1362  srmatrix A;
1363  srmatrix* M; //Originalmatrix
1364 
1365  private:
1366  srmatrix_slice(srmatrix& Mat, int sl1l, int sl1u, int sl2l, int sl2u) {
1367  A.lb1 = sl1l;
1368  A.lb2 = sl2l;
1369  A.ub1 = sl1u;
1370  A.ub2 = sl2u;
1371  A.m = sl1u-sl1l+1;
1372  A.n = sl2u-sl2l+1;
1373 
1374  //Kopieren der Werte aus A
1375  A.p = std::vector<int>(A.n+1, 0);
1376  A.ind.reserve(A.m + A.n);
1377  A.x.reserve(A.m + A.n);
1378 
1379  for(int i=0 ; i<A.n ; i++) {
1380  A.p[i+1] = A.p[i];
1381  for(int j=Mat.p[sl2l-Mat.lb2+i] ; j<Mat.p[sl2l-Mat.lb2+i+1] ; j++) {
1382  if(Mat.ind[j] >= sl1l-Mat.lb1 && Mat.ind[j] <= sl1u-Mat.lb1) {
1383  A.ind.push_back(Mat.ind[j]-(sl1l-Mat.lb1));
1384  A.x.push_back(Mat.x[j]);
1385  A.p[i+1]++;
1386  }
1387  }
1388  }
1389 
1390  //Zeiger auf A fuer Datenmanipulationen
1391  M = &Mat;
1392  }
1393 
1394  srmatrix_slice(const srmatrix& Mat, int sl1l, int sl1u, int sl2l, int sl2u) {
1395  A.lb1 = sl1l;
1396  A.lb2 = sl2l;
1397  A.ub1 = sl1u;
1398  A.ub2 = sl2u;
1399  A.m = sl1u-sl1l+1;
1400  A.n = sl2u-sl2l+1;
1401 
1402  //Kopieren der Werte aus A
1403  A.p = std::vector<int>(A.n+1, 0);
1404  A.ind.reserve(A.m + A.n);
1405  A.x.reserve(A.m + A.n);
1406 
1407  for(int i=0 ; i<A.n ; i++) {
1408  A.p[i+1] = A.p[i];
1409  for(int j=Mat.p[sl2l-Mat.lb2+i] ; j<Mat.p[sl2l-Mat.lb2+i+1] ; j++) {
1410  if(Mat.ind[j] >= sl1l-Mat.lb1 && Mat.ind[j] <= sl1u-Mat.lb1) {
1411  A.ind.push_back(Mat.ind[j]-(sl1l-Mat.lb1));
1412  A.x.push_back(Mat.x[j]);
1413  A.p[i+1]++;
1414  }
1415  }
1416  }
1417 
1418  //Zeiger auf A fuer Datenmanipulationen
1419  M = const_cast<srmatrix*>(&Mat);
1420  }
1421 
1422  public:
1425  return sl_ms_assign<srmatrix_slice, real, std::vector<real>::iterator, real>(*this,C);
1426  }
1427 
1430  return slsp_mm_assign<srmatrix_slice, srmatrix, std::vector<real>::iterator>(*this,C);
1431  }
1432 
1435  return slf_mm_assign<srmatrix_slice, rmatrix, std::vector<real>::iterator, real>(*this,C);
1436  }
1437 
1440  return slf_mm_assign<srmatrix_slice, rmatrix_slice, std::vector<real>::iterator, real>(*this,C);
1441  }
1442 
1445  *this = C.A;
1446  return *this;
1447  }
1448 
1451  *this = A*M.A;
1452  return *this;
1453  }
1454 
1457  *this = A*M;
1458  return *this;
1459  }
1460 
1463  *this = A*M;
1464  return *this;
1465  }
1466 
1469  *this = A*M;
1470  return *this;
1471  }
1472 
1475  *this = A*r;
1476  return *this;
1477  }
1478 
1481  *this = A/r;
1482  return *this;
1483  }
1484 
1487  *this = A+M.A;
1488  return *this;
1489  }
1490 
1493  *this = A+M;
1494  return *this;
1495  }
1496 
1499  *this = A+M;
1500  return *this;
1501  }
1502 
1505  *this = A+M;
1506  return *this;
1507  }
1508 
1511  *this = A-M.A;
1512  return *this;
1513  }
1514 
1517  *this = A-M;
1518  return *this;
1519  }
1520 
1523  *this = A-M;
1524  return *this;
1525  }
1526 
1529  *this = A-M;
1530  return *this;
1531  }
1532 
1534 
1538  const real operator()(const int i, const int j) const {
1539 #if(CXSC_INDEX_CHECK)
1540  if(i<A.lb1 || i>A.ub1 || j<A.lb2 || j>A.ub2)
1541  cxscthrow(ELEMENT_NOT_IN_VEC("srmatrix_slice::operator()(int, int)"));
1542 #endif
1543  real r = A(i,j);
1544  return r;
1545  }
1546 
1548 
1552  real& element(const int i, const int j) {
1553 #if(CXSC_INDEX_CHECK)
1554  if(i<A.lb1 || i>A.ub1 || j<A.lb2 || j>A.ub2)
1555  cxscthrow(ELEMENT_NOT_IN_VEC("srmatrix::element(int, int)"));
1556 #endif
1557  return M->element(i,j);
1558  }
1559 
1561  srmatrix_subv operator[](const int);
1563  srmatrix_subv operator[](const cxscmatrix_column&);
1565  const srmatrix_subv operator[](const int) const;
1567  const srmatrix_subv operator[](const cxscmatrix_column&) const;
1568 
1569  friend int Lb(const srmatrix_slice&, const int);
1570  friend int Ub(const srmatrix_slice&, const int);
1571  friend int RowLen(const srmatrix_slice&);
1572  friend int ColLen(const srmatrix_slice&);
1573 
1574  friend class srmatrix;
1575  friend class srmatrix_subv;
1576  friend class srvector;
1577  friend class scmatrix;
1578  friend class scmatrix_slice;
1579  friend class scmatrix_subv;
1580  friend class scvector;
1581  friend class simatrix;
1582  friend class simatrix_slice;
1583  friend class simatrix_subv;
1584  friend class sivector;
1585  friend class scimatrix;
1586  friend class scimatrix_slice;
1587  friend class scimatrix_subv;
1588  friend class scivector;
1589  friend class rmatrix;
1590  friend class imatrix;
1591  friend class cmatrix;
1592  friend class cimatrix;
1593 
1594 #include "matrix_friend_declarations.inl"
1595 };
1596 
1598  dat = new real[A.A.m*A.A.n];
1599  lb1 = A.A.lb1; lb2 = A.A.lb2; ub1 = A.A.ub1; ub2 = A.A.ub2;
1600  xsize = A.A.n;
1601  ysize = A.A.m;
1602  *this = 0.0;
1603  for(int j=0 ; j<A.A.n ; j++) {
1604  for(int k=A.A.p[j] ; k<A.A.p[j+1] ; k++) {
1605  dat[A.A.ind[k]*A.A.n+j] = A.A.x[k];
1606  }
1607  }
1608 }
1609 
1611 inline int Lb(const srmatrix_slice& S, const int i) {
1612  return Lb(S.A, i);
1613 }
1614 
1616 inline int Ub(const srmatrix_slice& S, const int i) {
1617  return Ub(S.A, i);
1618 }
1619 
1621 inline int RowLen(const srmatrix_slice& S) {
1622  return RowLen(S.A);
1623 }
1624 
1626 inline int ColLen(const srmatrix_slice& S) {
1627  return ColLen(S.A);
1628 }
1629 
1630 inline srmatrix_slice srmatrix::operator()(const int i, const int j, const int k, const int l) {
1631 #if(CXSC_INDEX_CHECK)
1632  if(i<lb1 || j>ub1 || k<lb2 || l>ub2)
1633  cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix::operator()(int, int, int, int)"));
1634 #endif
1635  return srmatrix_slice(*this, i, j, k, l);
1636 }
1637 
1638 inline const srmatrix_slice srmatrix::operator()(const int i, const int j, const int k, const int l) const {
1639 #if(CXSC_INDEX_CHECK)
1640  if(i<lb1 || j>ub1 || k<lb2 || l>ub2)
1641  cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix::operator()(int, int, int, int) const"));
1642 #endif
1643  return srmatrix_slice(*this, i, j, k, l);
1644 }
1645 
1647  *this = S.A;
1648  return *this;
1649 }
1650 
1652  m = S.A.m;
1653  n = S.A.n;
1654  lb1 = S.A.lb1;
1655  ub1 = S.A.ub1;
1656  lb2 = S.A.lb2;
1657  ub2 = S.A.ub2;
1658  *this = S.A;
1659 }
1660 
1662 inline srmatrix operator-(const srmatrix_slice& M) {
1663  return sp_m_negative<srmatrix,srmatrix>(M.A);
1664 }
1665 
1667 inline srmatrix operator+(const srmatrix_slice& M) {
1668  return M.A;
1669 }
1670 
1672 
1678 inline srmatrix operator*(const srmatrix_slice& M1, const srmatrix_slice& M2) {
1679  return spsp_mm_mult<srmatrix,srmatrix,srmatrix,sparse_dot,real>(M1.A,M2.A);
1680 }
1681 
1683 
1689 inline srmatrix operator*(const srmatrix_slice& M1, const srmatrix& M2) {
1690  return spsp_mm_mult<srmatrix,srmatrix,srmatrix,sparse_dot,real>(M1.A,M2);
1691 }
1692 
1694 
1700 inline srmatrix operator*(const srmatrix& M1, const srmatrix_slice& M2) {
1701  return spsp_mm_mult<srmatrix,srmatrix,srmatrix,sparse_dot,real>(M1,M2.A);
1702 }
1703 
1705 
1711 inline rmatrix operator*(const srmatrix_slice& M1, const rmatrix& M2) {
1712  return spf_mm_mult<srmatrix,rmatrix,rmatrix,sparse_dot>(M1.A,M2);
1713 }
1714 
1716 
1722 inline rmatrix operator*(const rmatrix& M1, const srmatrix_slice& M2) {
1723  return fsp_mm_mult<rmatrix,srmatrix,rmatrix,sparse_dot>(M1,M2.A);
1724 }
1725 
1727 
1733 inline rmatrix operator*(const srmatrix_slice& M1, const rmatrix_slice& M2) {
1734  return spf_mm_mult<srmatrix,rmatrix_slice,rmatrix,sparse_dot>(M1.A,M2);
1735 }
1736 
1738 
1744 inline rmatrix operator*(const rmatrix_slice& M1, const srmatrix_slice& M2) {
1745  return fsp_mm_mult<rmatrix,srmatrix,rmatrix,sparse_dot>(M1,M2.A);
1746 }
1747 
1749 
1755 inline srvector operator*(const srmatrix_slice& M, const srvector& v) {
1756  return spsp_mv_mult<srmatrix,srvector,srvector,sparse_dot,real>(M.A,v);
1757 }
1758 
1760 
1766 inline srvector operator*(const srmatrix_slice& M, const srvector_slice& v) {
1767  return spsl_mv_mult<srmatrix,srvector_slice,srvector,sparse_dot,real>(M.A,v);
1768 }
1769 
1771 
1777 inline rvector operator*(const srmatrix_slice& M, const rvector& v) {
1778  return spf_mv_mult<srmatrix,rvector,rvector,sparse_dot>(M.A,v);
1779 }
1780 
1782 
1788 inline rvector operator*(const srmatrix_slice& M, const rvector_slice& v) {
1789  return spf_mv_mult<srmatrix,rvector_slice,rvector,sparse_dot>(M.A,v);
1790 }
1791 
1793 inline srmatrix operator*(const srmatrix_slice& M, const real& r) {
1794  return sp_ms_mult<srmatrix,real,srmatrix>(M.A,r);
1795 }
1796 
1798 inline srmatrix operator/(const srmatrix_slice& M, const real& r) {
1799  return sp_ms_div<srmatrix,real,srmatrix>(M.A,r);
1800 }
1801 
1803 inline srmatrix operator*(const real& r, const srmatrix_slice& M) {
1804  return sp_sm_mult<real,srmatrix,srmatrix>(r,M.A);
1805 }
1806 
1808 inline srmatrix operator+(const srmatrix_slice& M1, const srmatrix_slice& M2) {
1809  return spsp_mm_add<srmatrix,srmatrix,srmatrix,real>(M1.A,M2.A);
1810 }
1811 
1813 inline srmatrix operator+(const srmatrix_slice& M1, const srmatrix& M2) {
1814  return spsp_mm_add<srmatrix,srmatrix,srmatrix,real>(M1.A,M2);
1815 }
1816 
1818 inline srmatrix operator+(const srmatrix& M1, const srmatrix_slice& M2) {
1819  return spsp_mm_add<srmatrix,srmatrix,srmatrix,real>(M1,M2.A);
1820 }
1821 
1823 inline rmatrix operator+(const srmatrix_slice& M1, const rmatrix& M2) {
1824  return spf_mm_add<srmatrix,rmatrix,rmatrix>(M1.A,M2);
1825 }
1826 
1828 inline rmatrix operator+(const rmatrix& M1, const srmatrix_slice& M2) {
1829  return fsp_mm_add<rmatrix,srmatrix,rmatrix>(M1,M2.A);
1830 }
1831 
1833 inline rmatrix operator+(const srmatrix_slice& M1, const rmatrix_slice& M2) {
1834  return spf_mm_add<srmatrix,rmatrix_slice,rmatrix>(M1.A,M2);
1835 }
1836 
1838 inline rmatrix operator+(const rmatrix_slice& M1, const srmatrix_slice& M2) {
1839  return fsp_mm_add<rmatrix_slice,srmatrix,rmatrix>(M1,M2.A);
1840 }
1841 
1843 inline srmatrix operator-(const srmatrix_slice& M1, const srmatrix_slice& M2) {
1844  return spsp_mm_sub<srmatrix,srmatrix,srmatrix,real>(M1.A,M2.A);
1845 }
1846 
1848 inline srmatrix operator-(const srmatrix_slice& M1, const srmatrix& M2) {
1849  return spsp_mm_sub<srmatrix,srmatrix,srmatrix,real>(M1.A,M2);
1850 }
1851 
1853 inline srmatrix operator-(const srmatrix& M1, const srmatrix_slice& M2) {
1854  return spsp_mm_sub<srmatrix,srmatrix,srmatrix,real>(M1,M2.A);
1855 }
1856 
1858 inline rmatrix operator-(const srmatrix_slice& M1, const rmatrix& M2) {
1859  return spf_mm_sub<srmatrix,rmatrix,rmatrix>(M1.A,M2);
1860 }
1861 
1863 inline rmatrix operator-(const rmatrix& M1, const srmatrix_slice& M2) {
1864  return fsp_mm_sub<rmatrix,srmatrix,rmatrix>(M1,M2.A);
1865 }
1866 
1868 inline rmatrix operator-(const srmatrix_slice& M1, const rmatrix_slice& M2) {
1869  return spf_mm_sub<srmatrix,rmatrix_slice,rmatrix>(M1.A,M2);
1870 }
1871 
1873 inline rmatrix operator-(const rmatrix_slice& M1, const srmatrix_slice& M2) {
1874  return fsp_mm_sub<rmatrix_slice,srmatrix,rmatrix>(M1,M2.A);
1875 }
1876 
1878  *this = rmatrix(M);
1879  return *this;
1880 }
1881 
1883  *this += M.A;
1884  return *this;
1885 }
1886 
1888  *this += M.A;
1889  return *this;
1890 }
1891 
1893  *this -= M.A;
1894  return *this;
1895 }
1896 
1898  *this -= M.A;
1899  return *this;
1900 }
1901 
1903  *this *= M.A;
1904  return *this;
1905 }
1906 
1908  *this *= M.A;
1909  return *this;
1910 }
1911 
1913 inline bool operator==(const srmatrix_slice& M1, const srmatrix_slice& M2) {
1914  return spsp_mm_comp(M1.A,M2.A);
1915 }
1916 
1918 inline bool operator==(const srmatrix_slice& M1, const srmatrix& M2) {
1919  return spsp_mm_comp(M1.A,M2);
1920 }
1921 
1923 inline bool operator==(const srmatrix& M1, const srmatrix_slice& M2) {
1924  return spsp_mm_comp(M1,M2.A);
1925 }
1926 
1928 inline bool operator==(const srmatrix_slice& M1, const rmatrix& M2) {
1929  return spf_mm_comp(M1.A,M2);
1930 }
1931 
1933 inline bool operator==(const rmatrix& M1, const srmatrix_slice& M2) {
1934  return fsp_mm_comp(M1,M2.A);
1935 }
1936 
1938 inline bool operator==(const rmatrix_slice& M1, const srmatrix_slice& M2) {
1939  return fsp_mm_comp(M1,M2.A);
1940 }
1941 
1943 inline bool operator==(const srmatrix_slice& M1, const rmatrix_slice& M2) {
1944  return spf_mm_comp(M1.A,M2);
1945 }
1946 
1948 inline bool operator!=(const srmatrix_slice& M1, const srmatrix_slice& M2) {
1949  return !spsp_mm_comp(M1.A,M2.A);
1950 }
1951 
1953 inline bool operator!=(const srmatrix_slice& M1, const srmatrix& M2) {
1954  return !spsp_mm_comp(M1.A,M2);
1955 }
1956 
1958 inline bool operator!=(const srmatrix& M1, const srmatrix_slice& M2) {
1959  return !spsp_mm_comp(M1,M2.A);
1960 }
1961 
1963 inline bool operator!=(const srmatrix_slice& M1, const rmatrix& M2) {
1964  return !spf_mm_comp(M1.A,M2);
1965 }
1966 
1968 inline bool operator!=(const rmatrix& M1, const srmatrix_slice& M2) {
1969  return !fsp_mm_comp(M1,M2.A);
1970 }
1971 
1973 inline bool operator!=(const rmatrix_slice& M1, const srmatrix_slice& M2) {
1974  return !fsp_mm_comp(M1,M2.A);
1975 }
1976 
1978 inline bool operator!=(const srmatrix_slice& M1, const rmatrix_slice& M2) {
1979  return !spf_mm_comp(M1.A,M2);
1980 }
1981 
1983 inline bool operator<(const srmatrix_slice& M1, const srmatrix_slice& M2) {
1984  return spsp_mm_less<srmatrix,srmatrix,real>(M1.A,M2.A);
1985 }
1986 
1988 inline bool operator<(const srmatrix_slice& M1, const srmatrix& M2) {
1989  return spsp_mm_less<srmatrix,srmatrix,real>(M1.A,M2);
1990 }
1991 
1993 inline bool operator<(const srmatrix& M1, const srmatrix_slice& M2) {
1994  return spsp_mm_less<srmatrix,srmatrix,real>(M1,M2.A);
1995 }
1996 
1998 inline bool operator<(const srmatrix_slice& M1, const rmatrix& M2) {
1999  return spf_mm_less<srmatrix,rmatrix,real>(M1.A,M2);
2000 }
2001 
2003 inline bool operator<(const rmatrix& M1, const srmatrix_slice& M2) {
2004  return fsp_mm_less<rmatrix,srmatrix,real>(M1,M2.A);
2005 }
2006 
2008 inline bool operator<(const rmatrix_slice& M1, const srmatrix_slice& M2) {
2009  return fsp_mm_less<rmatrix_slice,srmatrix,real>(M1,M2.A);
2010 }
2011 
2013 inline bool operator<(const srmatrix_slice& M1, const rmatrix_slice& M2) {
2014  return spf_mm_less<srmatrix,rmatrix_slice,real>(M1.A,M2);
2015 }
2016 
2018 inline bool operator<=(const srmatrix_slice& M1, const srmatrix_slice& M2) {
2019  return spsp_mm_leq<srmatrix,srmatrix,real>(M1.A,M2.A);
2020 }
2021 
2023 inline bool operator<=(const srmatrix_slice& M1, const srmatrix& M2) {
2024  return spsp_mm_leq<srmatrix,srmatrix,real>(M1.A,M2);
2025 }
2026 
2028 inline bool operator<=(const srmatrix& M1, const srmatrix_slice& M2) {
2029  return spsp_mm_leq<srmatrix,srmatrix,real>(M1,M2.A);
2030 }
2031 
2033 inline bool operator<=(const srmatrix_slice& M1, const rmatrix& M2) {
2034  return spf_mm_leq<srmatrix,rmatrix,real>(M1.A,M2);
2035 }
2036 
2038 inline bool operator<=(const rmatrix& M1, const srmatrix_slice& M2) {
2039  return fsp_mm_leq<rmatrix,srmatrix,real>(M1,M2.A);
2040 }
2041 
2043 inline bool operator<=(const rmatrix_slice& M1, const srmatrix_slice& M2) {
2044  return fsp_mm_leq<rmatrix_slice,srmatrix,real>(M1,M2.A);
2045 }
2046 
2048 inline bool operator<=(const srmatrix_slice& M1, const rmatrix_slice& M2) {
2049  return spf_mm_leq<srmatrix,rmatrix_slice,real>(M1.A,M2);
2050 }
2051 
2053 inline bool operator>(const srmatrix_slice& M1, const srmatrix_slice& M2) {
2054  return spsp_mm_greater<srmatrix,srmatrix,real>(M1.A,M2.A);
2055 }
2056 
2058 inline bool operator>(const srmatrix_slice& M1, const srmatrix& M2) {
2059  return spsp_mm_greater<srmatrix,srmatrix,real>(M1.A,M2);
2060 }
2061 
2063 inline bool operator>(const srmatrix& M1, const srmatrix_slice& M2) {
2064  return spsp_mm_greater<srmatrix,srmatrix,real>(M1,M2.A);
2065 }
2066 
2068 inline bool operator>(const srmatrix_slice& M1, const rmatrix& M2) {
2069  return spf_mm_greater<srmatrix,rmatrix,real>(M1.A,M2);
2070 }
2071 
2073 inline bool operator>(const rmatrix& M1, const srmatrix_slice& M2) {
2074  return fsp_mm_greater<rmatrix,srmatrix,real>(M1,M2.A);
2075 }
2076 
2078 inline bool operator>(const rmatrix_slice& M1, const srmatrix_slice& M2) {
2079  return fsp_mm_greater<rmatrix_slice,srmatrix,real>(M1,M2.A);
2080 }
2081 
2083 inline bool operator>(const srmatrix_slice& M1, const rmatrix_slice& M2) {
2084  return spf_mm_greater<srmatrix,rmatrix_slice,real>(M1.A,M2);
2085 }
2086 
2088 inline bool operator>=(const srmatrix_slice& M1, const srmatrix_slice& M2) {
2089  return spsp_mm_geq<srmatrix,srmatrix,real>(M1.A,M2.A);
2090 }
2091 
2093 inline bool operator>=(const srmatrix_slice& M1, const srmatrix& M2) {
2094  return spsp_mm_geq<srmatrix,srmatrix,real>(M1.A,M2);
2095 }
2096 
2098 inline bool operator>=(const srmatrix& M1, const srmatrix_slice& M2) {
2099  return spsp_mm_geq<srmatrix,srmatrix,real>(M1,M2.A);
2100 }
2101 
2103 inline bool operator>=(const srmatrix_slice& M1, const rmatrix& M2) {
2104  return spf_mm_geq<srmatrix,rmatrix,real>(M1.A,M2);
2105 }
2106 
2108 inline bool operator>=(const rmatrix& M1, const srmatrix_slice& M2) {
2109  return fsp_mm_geq<rmatrix,srmatrix,real>(M1,M2.A);
2110 }
2111 
2113 inline bool operator>=(const rmatrix_slice& M1, const srmatrix_slice& M2) {
2114  return fsp_mm_geq<rmatrix_slice,srmatrix,real>(M1,M2.A);
2115 }
2116 
2118 inline bool operator>=(const srmatrix_slice& M1, const rmatrix_slice& M2) {
2119  return spf_mm_geq<srmatrix,rmatrix_slice,real>(M1.A,M2);
2120 }
2121 
2123 inline bool operator!(const srmatrix_slice& M) {
2124  return sp_m_not(M.A);
2125 }
2126 
2128 
2133 inline std::ostream& operator<<(std::ostream& os, const srmatrix_slice& M) {
2134  return sp_m_output<srmatrix,real>(os, M.A);
2135 }
2136 
2138 
2143 inline std::istream& operator>>(std::istream& is, srmatrix_slice& M) {
2144  srmatrix tmp(M.A.m,M.A.n);
2145  sp_m_input<srmatrix,real>(is, tmp);
2146  M = tmp;
2147  return is;
2148 }
2149 
2150 
2152 
2158  private:
2159  srmatrix_slice dat;
2160  bool row;
2161  int index;
2162 
2163 
2164  srmatrix_subv(srmatrix& A, bool r, int i, int j, int k, int l) : dat(A,i,j,k,l), row(r) {
2165  if(row) index=i; else index=k;
2166  }
2167 
2168  srmatrix_subv(const srmatrix& A, bool r, int i, int j, int k, int l) : dat(A,i,j,k,l), row(r){
2169  if(row) index=i; else index=k;
2170  }
2171 
2172  public:
2173 
2175 
2179  real& operator[](const int i) {
2180  if(row) {
2181 #if(CXSC_INDEX_CHECK)
2182  if(i<dat.A.lb2 || i>dat.A.ub2)
2183  cxscthrow(ELEMENT_NOT_IN_VEC("srmatrix_subv::operator[](int)"));
2184 #endif
2185  return dat.element(index,i);
2186  } else {
2187 #if(CXSC_INDEX_CHECK)
2188  if(i<dat.A.lb1 || i>dat.A.ub1)
2189  cxscthrow(ELEMENT_NOT_IN_VEC("srmatrix_subv::operator[](int)"));
2190 #endif
2191  return dat.element(i,index);
2192  }
2193  }
2194 
2196 
2199  const real operator[](const int i) const {
2200  if(row) {
2201 #if(CXSC_INDEX_CHECK)
2202  if(i<dat.A.lb2 || i>dat.A.ub2)
2203  cxscthrow(ELEMENT_NOT_IN_VEC("srmatrix_subv::operator[](int)"));
2204 #endif
2205  return dat(index,i);
2206  } else {
2207 #if(CXSC_INDEX_CHECK)
2208  if(i<dat.A.lb1 || i>dat.A.ub1)
2209  cxscthrow(ELEMENT_NOT_IN_VEC("srmatrix_subv::operator[](int)"));
2210 #endif
2211  return dat(i,index);
2212  }
2213  }
2214 
2217  return svsp_vv_assign(*this,srvector(v));
2218  }
2219 
2222  return sv_vs_assign(*this,v);
2223  }
2224 
2227  return svsp_vv_assign(*this,v);
2228  }
2229 
2232  return svsl_vv_assign(*this,v);
2233  }
2234 
2237  return svf_vv_assign(*this,v);
2238  }
2239 
2242  return svf_vv_assign(*this,v);
2243  }
2244 
2246  srmatrix_subv& operator*=(const real&);
2248  srmatrix_subv& operator/=(const real&);
2254  srmatrix_subv& operator+=(const rvector&);
2262  srmatrix_subv& operator-=(const rvector&);
2265 
2266  friend srvector operator-(const srmatrix_subv&);
2267 
2268  friend std::istream& operator>>(std::istream&, srmatrix_subv&);
2269 
2270  friend int Lb(const srmatrix_subv&);
2271  friend int Ub(const srmatrix_subv&);
2272  friend int VecLen(const srmatrix_subv&);
2273 
2274  friend class srvector;
2275  friend class srmatrix;
2276  friend class srmatrix_slice;
2277  friend class scvector;
2278  friend class scmatrix;
2279  friend class scmatrix_slice;
2280  friend class sivector;
2281  friend class simatrix;
2282  friend class simatrix_slice;
2283  friend class scivector;
2284  friend class scimatrix;
2285  friend class scimatrix_slice;
2286 
2287 #include "vector_friend_declarations.inl"
2288 };
2289 
2291 inline int Lb(const srmatrix_subv& S) {
2292  if(S.row)
2293  return Lb(S.dat, 2);
2294  else
2295  return Lb(S.dat, 1);
2296 }
2297 
2299 inline int Ub(const srmatrix_subv& S) {
2300  if(S.row)
2301  return Ub(S.dat, 2);
2302  else
2303  return Ub(S.dat, 1);
2304 }
2305 
2307 inline int VecLen(const srmatrix_subv& S) {
2308  return Ub(S)-Lb(S)+1;
2309 }
2310 
2312 inline std::ostream& operator<<(std::ostream& os, const srmatrix_subv& v) {
2313  os << srvector(v);
2314  return os;
2315 }
2316 
2318 inline std::istream& operator>>(std::istream& is, srmatrix_subv& v) {
2319  int n = 0;
2320  if(v.row) n=v.dat.A.n; else n=v.dat.A.m;
2321  srvector tmp(n);
2322  is >> tmp;
2323  v = tmp;
2324  return is;
2325 }
2326 
2327 inline srmatrix_subv srmatrix::operator[](const cxscmatrix_column& c) {
2328 #if(CXSC_INDEX_CHECK)
2329  if(c.col()<lb2 || c.col()>ub2)
2330  cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix::operator[](const cxscmatrix_column&)"));
2331 #endif
2332  return srmatrix_subv(*this, false, lb1, ub1, c.col(), c.col());
2333 }
2334 
2336 #if(CXSC_INDEX_CHECK)
2337  if(i<lb1 || i>ub1)
2338  cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix::operator[](const int)"));
2339 #endif
2340  return srmatrix_subv(*this, true, i, i, lb2, ub2);
2341 }
2342 
2343 inline const srmatrix_subv srmatrix::operator[](const cxscmatrix_column& c) const {
2344 #if(CXSC_INDEX_CHECK)
2345  if(c.col()<lb2 || c.col()>ub2)
2346  cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix::operator[](const cxscmatrix_column&)"));
2347 #endif
2348  return srmatrix_subv(*this, false, lb1, ub1, c.col(), c.col());
2349 }
2350 
2351 inline const srmatrix_subv srmatrix::operator[](const int i) const {
2352 #if(CXSC_INDEX_CHECK)
2353  if(i<lb1 || i>ub1)
2354  cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix::operator[](const int)"));
2355 #endif
2356  return srmatrix_subv(*this, true, i, i, lb2, ub2);
2357 }
2358 
2360 #if(CXSC_INDEX_CHECK)
2361  if(i<A.lb1 || i>A.ub1)
2362  cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix_slice::operator[](const int)"));
2363 #endif
2364  return srmatrix_subv(*M, true, i, i, A.lb2, A.ub2);
2365 }
2366 
2367 inline srmatrix_subv srmatrix_slice::operator[](const cxscmatrix_column& c) {
2368 #if(CXSC_INDEX_CHECK)
2369  if(c.col()<A.lb2 || c.col()>A.ub2)
2370  cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix_slice::operator[](const cxscmatrix_column&)"));
2371 #endif
2372  return srmatrix_subv(*M, false, A.lb1, A.ub1, c.col(), c.col());
2373 }
2374 
2375 inline const srmatrix_subv srmatrix_slice::operator[](const int i) const {
2376 #if(CXSC_INDEX_CHECK)
2377  if(i<A.lb1 || i>A.ub1)
2378  cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix_slice::operator[](const int)"));
2379 #endif
2380  return srmatrix_subv(*M, true, i, i, A.lb2, A.ub2);
2381 }
2382 
2383 inline const srmatrix_subv srmatrix_slice::operator[](const cxscmatrix_column& c) const {
2384 #if(CXSC_INDEX_CHECK)
2385  if(c.col()<A.lb2 || c.col()>A.ub2)
2386  cxscthrow(ROW_OR_COL_NOT_IN_MAT("srmatrix_slice::operator[](const cxscmatrix_column&)"));
2387 #endif
2388  return srmatrix_subv(*M, false, A.lb1, A.ub1, c.col(), c.col());
2389 }
2390 
2392  p.reserve(A.dat.A.get_nnz());
2393  x.reserve(A.dat.A.get_nnz());
2394 
2395  if(A.row) {
2396  lb = A.dat.A.lb2;
2397  ub = A.dat.A.ub2;
2398  n = ub-lb+1;
2399 
2400  for(int j=0 ; j<n ; j++) {
2401  for(int k=A.dat.A.p[j] ; k<A.dat.A.p[j+1] ; k++) {
2402  p.push_back(j);
2403  x.push_back(A.dat.A.x[k]);
2404  }
2405  }
2406 
2407  } else {
2408  lb = A.dat.A.lb1;
2409  ub = A.dat.A.ub1;
2410  n = ub-lb+1;
2411 
2412  for(unsigned int k=0 ; k<A.dat.A.ind.size() ; k++) {
2413  p.push_back(A.dat.A.ind[k]);
2414  x.push_back(A.dat.A.x[k]);
2415  }
2416  }
2417 }
2418 
2420 inline srvector operator-(const srmatrix_subv& v) {
2421  srvector s(v);
2422  return -s;
2423 }
2424 
2426 inline srvector operator*(const srmatrix_subv& v1, const real& v2) {
2427  return srvector(v1) * v2;
2428 }
2429 
2431 inline srvector operator/(const srmatrix_subv& v1, const real& v2) {
2432  return srvector(v1) / v2;
2433 }
2434 
2436 inline srvector operator*(const real& v1, const srmatrix_subv& v2) {
2437  return v1 * srvector(v2);
2438 }
2439 
2441 
2447 inline real operator*(const srmatrix_subv& v1, const srvector& v2) {
2448  return srvector(v1) * v2;
2449 }
2450 
2452 
2458 inline real operator*(const srmatrix_subv& v1, const srvector_slice& v2) {
2459  return srvector(v1) * v2;
2460 }
2461 
2463 
2469 inline real operator*(const srmatrix_subv& v1, const rvector& v2) {
2470  return srvector(v1) * v2;
2471 }
2472 
2474 
2480 inline real operator*(const srmatrix_subv& v1, const rvector_slice& v2) {
2481  return srvector(v1) * v2;
2482 }
2483 
2485 
2491 inline real operator*(const srvector& v1, const srmatrix_subv& v2) {
2492  return v1 * srvector(v2);
2493 }
2494 
2496 
2502 inline real operator*(const srvector_slice& v1, const srmatrix_subv& v2) {
2503  return v1 * srvector(v2);
2504 }
2505 
2507 
2513 inline real operator*(const rvector& v1, const srmatrix_subv& v2) {
2514  return v1 * srvector(v2);
2515 }
2516 
2518 
2524 inline real operator*(const rvector_slice& v1, const srmatrix_subv& v2) {
2525  return v1 * srvector(v2);
2526 }
2527 
2529 inline srvector operator+(const srmatrix_subv& v1, const srvector& v2) {
2530  return srvector(v1) + v2;
2531 }
2532 
2534 inline srvector operator+(const srmatrix_subv& v1, const srvector_slice& v2) {
2535  return srvector(v1) + v2;
2536 }
2537 
2539 inline rvector operator+(const srmatrix_subv& v1, const rvector& v2) {
2540  return srvector(v1) + v2;
2541 }
2542 
2544 inline rvector operator+(const srmatrix_subv& v1, const rvector_slice& v2) {
2545  return srvector(v1) + v2;
2546 }
2547 
2549 inline srvector operator+(const srvector& v1, const srmatrix_subv& v2) {
2550  return v1 + srvector(v2);
2551 }
2552 
2554 inline srvector operator+(const srvector_slice& v1, const srmatrix_subv& v2) {
2555  return v1 + srvector(v2);
2556 }
2557 
2559 inline rvector operator+(const rvector& v1, const srmatrix_subv& v2) {
2560  return v1 + srvector(v2);
2561 }
2562 
2564 inline rvector operator+(const rvector_slice& v1, const srmatrix_subv& v2) {
2565  return v1 + srvector(v2);
2566 }
2567 
2569 inline srvector operator-(const srmatrix_subv& v1, const srvector& v2) {
2570  return srvector(v1) - v2;
2571 }
2572 
2574 inline srvector operator-(const srmatrix_subv& v1, const srvector_slice& v2) {
2575  return srvector(v1) - v2;
2576 }
2577 
2579 inline rvector operator-(const srmatrix_subv& v1, const rvector& v2) {
2580  return srvector(v1) - v2;
2581 }
2582 
2584 inline rvector operator-(const srmatrix_subv& v1, const rvector_slice& v2) {
2585  return srvector(v1) - v2;
2586 }
2587 
2589 inline srvector operator-(const srvector& v1, const srmatrix_subv& v2) {
2590  return v1 - srvector(v2);
2591 }
2592 
2594 inline srvector operator-(const srvector_slice& v1, const srmatrix_subv& v2) {
2595  return v1 - srvector(v2);
2596 }
2597 
2599 inline rvector operator-(const rvector& v1, const srmatrix_subv& v2) {
2600  return v1 - srvector(v2);
2601 }
2602 
2604 inline rvector operator-(const rvector_slice& v1, const srmatrix_subv& v2) {
2605  return v1 - srvector(v2);
2606 }
2607 
2609  *this = *this * v;
2610  return *this;
2611 }
2612 
2614  *this = *this * v;
2615  return *this;
2616 }
2617 
2619  *this = *this + v;
2620  return *this;
2621 }
2622 
2624  *this = *this + v;
2625  return *this;
2626 }
2627 
2629  *this = *this + v;
2630  return *this;
2631 }
2632 
2634  *this = *this + v;
2635  return *this;
2636 }
2637 
2639  *this = *this - v;
2640  return *this;
2641 }
2642 
2644  *this = *this - v;
2645  return *this;
2646 }
2647 
2649  *this = *this - v;
2650  return *this;
2651 }
2652 
2654  *this = *this - v;
2655  return *this;
2656 }
2657 
2659  *this += rvector(v);
2660  return *this;
2661 }
2662 
2664  *this += rvector(v);
2665  return *this;
2666 }
2667 
2669  *this += rvector(v);
2670  return *this;
2671 }
2672 
2674  *this -= rvector(v);
2675  return *this;
2676 }
2677 
2679  *this = rvector(v);
2680  return *this;
2681 }
2682 
2684  *this = rvector(v);
2685  return *this;
2686 }
2687 
2689  *this = rvector(v);
2690  return *this;
2691 }
2692 
2694 inline bool operator==(const srmatrix_subv& v1, const srvector& v2) {
2695  return srvector(v1) == v2;
2696 }
2697 
2699 inline bool operator==(const srmatrix_subv& v1, const srvector_slice& v2) {
2700  return srvector(v1) == v2;
2701 }
2702 
2704 inline bool operator==(const srmatrix_subv& v1, const rvector& v2) {
2705  return srvector(v1) == v2;
2706 }
2707 
2709 inline bool operator==(const srmatrix_subv& v1, const rvector_slice& v2) {
2710  return srvector(v1) == v2;
2711 }
2712 
2714 inline bool operator==(const srvector& v1, const srmatrix_subv& v2) {
2715  return v1 == srvector(v2);
2716 }
2717 
2719 inline bool operator==(const srvector_slice& v1, const srmatrix_subv& v2) {
2720  return v1 == srvector(v2);
2721 }
2722 
2724 inline bool operator==(const rvector& v1, const srmatrix_subv& v2) {
2725  return v1 == srvector(v2);
2726 }
2727 
2729 inline bool operator==(const rvector_slice& v1, const srmatrix_subv& v2) {
2730  return v1 == srvector(v2);
2731 }
2732 
2734 inline bool operator!=(const srmatrix_subv& v1, const srvector& v2) {
2735  return srvector(v1) != v2;
2736 }
2737 
2739 inline bool operator!=(const srmatrix_subv& v1, const srvector_slice& v2) {
2740  return srvector(v1) != v2;
2741 }
2742 
2744 inline bool operator!=(const srmatrix_subv& v1, const rvector& v2) {
2745  return srvector(v1) != v2;
2746 }
2747 
2749 inline bool operator!=(const srmatrix_subv& v1, const rvector_slice& v2) {
2750  return srvector(v1) != v2;
2751 }
2752 
2754 inline bool operator!=(const srvector& v1, const srmatrix_subv& v2) {
2755  return v1 != srvector(v2);
2756 }
2757 
2759 inline bool operator!=(const srvector_slice& v1, const srmatrix_subv& v2) {
2760  return v1 != srvector(v2);
2761 }
2762 
2764 inline bool operator!=(const rvector& v1, const srmatrix_subv& v2) {
2765  return v1 != srvector(v2);
2766 }
2767 
2769 inline bool operator!=(const rvector_slice& v1, const srmatrix_subv& v2) {
2770  return v1 != srvector(v2);
2771 }
2772 
2774 inline bool operator<(const srmatrix_subv& v1, const srvector& v2) {
2775  return srvector(v1) < v2;
2776 }
2777 
2779 inline bool operator<(const srmatrix_subv& v1, const srvector_slice& v2) {
2780  return srvector(v1) < v2;
2781 }
2782 
2784 inline bool operator<(const srmatrix_subv& v1, const rvector& v2) {
2785  return srvector(v1) < v2;
2786 }
2787 
2789 inline bool operator<(const srmatrix_subv& v1, const rvector_slice& v2) {
2790  return srvector(v1) < v2;
2791 }
2792 
2794 inline bool operator<(const srvector& v1, const srmatrix_subv& v2) {
2795  return v1 < srvector(v2);
2796 }
2797 
2799 inline bool operator<(const srvector_slice& v1, const srmatrix_subv& v2) {
2800  return v1 < srvector(v2);
2801 }
2802 
2804 inline bool operator<(const rvector& v1, const srmatrix_subv& v2) {
2805  return v1 < srvector(v2);
2806 }
2807 
2809 inline bool operator<(const rvector_slice& v1, const srmatrix_subv& v2) {
2810  return v1 < srvector(v2);
2811 }
2812 
2814 inline bool operator<=(const srmatrix_subv& v1, const srvector& v2) {
2815  return srvector(v1) <= v2;
2816 }
2817 
2819 inline bool operator<=(const srmatrix_subv& v1, const srvector_slice& v2) {
2820  return srvector(v1) <= v2;
2821 }
2822 
2824 inline bool operator<=(const srmatrix_subv& v1, const rvector& v2) {
2825  return srvector(v1) <= v2;
2826 }
2827 
2829 inline bool operator<=(const srmatrix_subv& v1, const rvector_slice& v2) {
2830  return srvector(v1) <= v2;
2831 }
2832 
2834 inline bool operator<=(const srvector& v1, const srmatrix_subv& v2) {
2835  return v1 <= srvector(v2);
2836 }
2837 
2839 inline bool operator<=(const srvector_slice& v1, const srmatrix_subv& v2) {
2840  return v1 <= srvector(v2);
2841 }
2842 
2844 inline bool operator<=(const rvector& v1, const srmatrix_subv& v2) {
2845  return v1 <= srvector(v2);
2846 }
2847 
2849 inline bool operator<=(const rvector_slice& v1, const srmatrix_subv& v2) {
2850  return v1 <= srvector(v2);
2851 }
2852 
2854 inline bool operator>(const srmatrix_subv& v1, const srvector& v2) {
2855  return srvector(v1) > v2;
2856 }
2857 
2859 inline bool operator>(const srmatrix_subv& v1, const srvector_slice& v2) {
2860  return srvector(v1) > v2;
2861 }
2862 
2864 inline bool operator>(const srmatrix_subv& v1, const rvector& v2) {
2865  return srvector(v1) > v2;
2866 }
2867 
2869 inline bool operator>(const srmatrix_subv& v1, const rvector_slice& v2) {
2870  return srvector(v1) > v2;
2871 }
2872 
2874 inline bool operator>(const srvector& v1, const srmatrix_subv& v2) {
2875  return v1 > srvector(v2);
2876 }
2877 
2879 inline bool operator>(const srvector_slice& v1, const srmatrix_subv& v2) {
2880  return v1 > srvector(v2);
2881 }
2882 
2884 inline bool operator>(const rvector& v1, const srmatrix_subv& v2) {
2885  return v1 > srvector(v2);
2886 }
2887 
2889 inline bool operator>(const rvector_slice& v1, const srmatrix_subv& v2) {
2890  return v1 > srvector(v2);
2891 }
2892 
2894 inline bool operator>=(const srmatrix_subv& v1, const srvector& v2) {
2895  return srvector(v1) >= v2;
2896 }
2897 
2899 inline bool operator>=(const srmatrix_subv& v1, const srvector_slice& v2) {
2900  return srvector(v1) >= v2;
2901 }
2902 
2904 inline bool operator>=(const srmatrix_subv& v1, const rvector& v2) {
2905  return srvector(v1) >= v2;
2906 }
2907 
2909 inline bool operator>=(const srmatrix_subv& v1, const rvector_slice& v2) {
2910  return srvector(v1) >= v2;
2911 }
2912 
2914 inline bool operator>=(const srvector& v1, const srmatrix_subv& v2) {
2915  return v1 >= srvector(v2);
2916 }
2917 
2919 inline bool operator>=(const srvector_slice& v1, const srmatrix_subv& v2) {
2920  return v1 >= srvector(v2);
2921 }
2922 
2924 inline bool operator>=(const rvector& v1, const srmatrix_subv& v2) {
2925  return v1 >= srvector(v2);
2926 }
2927 
2929 inline bool operator>=(const rvector_slice& v1, const srmatrix_subv& v2) {
2930  return v1 >= srvector(v2);
2931 }
2932 
2934 inline bool operator!(const srmatrix_subv& v) {
2935  return sv_v_not(v);
2936 }
2937 
2939 
2942 inline void accumulate(dotprecision& dot, const srmatrix_subv& v1, const srmatrix_subv& v2) {
2943  spsp_vv_accu<dotprecision,srvector,srvector,sparse_dot>(dot, srvector(v1), srvector(v2));
2944 }
2945 
2947 
2950 inline void accumulate(dotprecision& dot, const srmatrix_subv& v1, const srvector& v2) {
2951  spsp_vv_accu<dotprecision,srvector,srvector,sparse_dot>(dot, srvector(v1), v2);
2952 }
2953 
2955 
2958 inline void accumulate(dotprecision& dot, const srmatrix_subv& v1, const srvector_slice& v2) {
2959  spsl_vv_accu<dotprecision,srvector,srvector_slice,sparse_dot>(dot, srvector(v1), v2);
2960 }
2961 
2963 
2966 inline void accumulate(dotprecision& dot, const srmatrix_subv& v1, const rvector& v2) {
2967  spf_vv_accu<dotprecision,srvector,rvector,sparse_dot>(dot, srvector(v1), v2);
2968 }
2969 
2971 
2974 inline void accumulate(dotprecision& dot, const srmatrix_subv& v1, const rvector_slice& v2) {
2975  spf_vv_accu<dotprecision,srvector,rvector_slice,sparse_dot>(dot, srvector(v1), v2);
2976 }
2977 
2979 
2982 inline void accumulate(dotprecision& dot, const srvector& v1, const srmatrix_subv& v2) {
2983  spsp_vv_accu<dotprecision,srvector,srvector,sparse_dot>(dot, v1, srvector(v2));
2984 }
2985 
2987 
2990 inline void accumulate(dotprecision& dot, const srvector_slice& v1, const srmatrix_subv& v2) {
2991  slsp_vv_accu<dotprecision,srvector_slice,srvector,sparse_dot>(dot, v1, srvector(v2));
2992 }
2993 
2995 
2998 inline void accumulate(dotprecision& dot, const rvector& v1, const srmatrix_subv& v2) {
2999  fsp_vv_accu<dotprecision,rvector,srvector,sparse_dot>(dot, v1, srvector(v2));
3000 }
3001 
3003 
3006 inline void accumulate(dotprecision& dot, const rvector_slice& v1, const srmatrix_subv& v2) {
3007  fsp_vv_accu<dotprecision,rvector_slice,srvector,sparse_dot>(dot, v1, srvector(v2));
3008 }
3009 
3011 
3016 inline void accumulate_approx(dotprecision& dot, const srmatrix_subv& v1, const srmatrix_subv& v2) {
3017  spsp_vv_accuapprox<dotprecision,srvector,srvector,sparse_dot>(dot, srvector(v1), srvector(v2));
3018 }
3019 
3021 
3026 inline void accumulate_approx(dotprecision& dot, const srmatrix_subv& v1, const srvector& v2) {
3027  spsp_vv_accuapprox<dotprecision,srvector,srvector,sparse_dot>(dot, srvector(v1), v2);
3028 }
3029 
3031 
3036 inline void accumulate_approx(dotprecision& dot, const srmatrix_subv& v1, const srvector_slice& v2) {
3037  spsl_vv_accuapprox<dotprecision,srvector,srvector_slice,sparse_dot>(dot, srvector(v1), v2);
3038 }
3039 
3041 
3046 inline void accumulate_approx(dotprecision& dot, const srmatrix_subv& v1, const rvector& v2) {
3047  spf_vv_accuapprox<dotprecision,srvector,rvector,sparse_dot>(dot, srvector(v1), v2);
3048 }
3049 
3051 
3056 inline void accumulate_approx(dotprecision& dot, const srmatrix_subv& v1, const rvector_slice& v2) {
3057  spf_vv_accuapprox<dotprecision,srvector,rvector_slice,sparse_dot>(dot, srvector(v1), v2);
3058 }
3059 
3061 
3066 inline void accumulate_approx(dotprecision& dot, const srvector& v1, const srmatrix_subv& v2) {
3067  spsp_vv_accuapprox<dotprecision,srvector,srvector,sparse_dot>(dot, v1, srvector(v2));
3068 }
3069 
3071 
3076 inline void accumulate_approx(dotprecision& dot, const srvector_slice& v1, const srmatrix_subv& v2) {
3077  slsp_vv_accuapprox<dotprecision,srvector_slice,srvector,sparse_dot>(dot, v1, srvector(v2));
3078 }
3079 
3081 
3086 inline void accumulate_approx(dotprecision& dot, const rvector& v1, const srmatrix_subv& v2) {
3087  fsp_vv_accuapprox<dotprecision,rvector,srvector,sparse_dot>(dot, v1, srvector(v2));
3088 }
3089 
3091 
3096 inline void accumulate_approx(dotprecision& dot, const rvector_slice& v1, const srmatrix_subv& v2) {
3097  fsp_vv_accuapprox<dotprecision,rvector_slice,srvector,sparse_dot>(dot, v1, srvector(v2));
3098 }
3099 
3101 
3104 inline void accumulate(idotprecision& dot, const srmatrix_subv& v1, const srmatrix_subv& v2) {
3105  dotprecision tmp(0.0);
3106  tmp.set_k(dot.get_k());
3107  accumulate(tmp,srvector(v1),srvector(v2));
3108  dot += tmp;
3109 }
3110 
3112 
3115 inline void accumulate(idotprecision& dot, const srmatrix_subv& v1, const srvector& v2) {
3116  dotprecision tmp(0.0);
3117  tmp.set_k(dot.get_k());
3118  accumulate(tmp,srvector(v1),v2);
3119  dot += tmp;
3120 }
3121 
3123 
3126 inline void accumulate(idotprecision& dot, const srmatrix_subv& v1, const srvector_slice& v2) {
3127  dotprecision tmp(0.0);
3128  tmp.set_k(dot.get_k());
3129  accumulate(tmp,srvector(v1),v2);
3130  dot += tmp;
3131 }
3132 
3134 
3137 inline void accumulate(idotprecision& dot, const srmatrix_subv& v1, const rvector& v2) {
3138  dotprecision tmp(0.0);
3139  tmp.set_k(dot.get_k());
3140  accumulate(tmp,srvector(v1),v2);
3141  dot += tmp;
3142 }
3143 
3145 
3148 inline void accumulate(idotprecision& dot, const srmatrix_subv& v1, const rvector_slice& v2) {
3149  dotprecision tmp(0.0);
3150  tmp.set_k(dot.get_k());
3151  accumulate(tmp,srvector(v1),v2);
3152  dot += tmp;
3153 }
3154 
3156 
3159 inline void accumulate(idotprecision& dot, const srvector& v1, const srmatrix_subv& v2) {
3160  dotprecision tmp(0.0);
3161  tmp.set_k(dot.get_k());
3162  accumulate(tmp,v1,srvector(v2));
3163  dot += tmp;
3164 }
3165 
3167 
3170 inline void accumulate(idotprecision& dot, const srvector_slice& v1, const srmatrix_subv& v2) {
3171  dotprecision tmp(0.0);
3172  tmp.set_k(dot.get_k());
3173  accumulate(tmp,v1,srvector(v2));
3174  dot += tmp;
3175 }
3176 
3178 
3181 inline void accumulate(idotprecision& dot, const rvector& v1, const srmatrix_subv& v2) {
3182  dotprecision tmp(0.0);
3183  tmp.set_k(dot.get_k());
3184  accumulate(tmp,v1,srvector(v2));
3185  dot += tmp;
3186 }
3187 
3189 
3192 inline void accumulate(idotprecision& dot, const rvector_slice& v1, const srmatrix_subv& v2) {
3193  dotprecision tmp(0.0);
3194  tmp.set_k(dot.get_k());
3195  accumulate(tmp,v1,srvector(v2));
3196  dot += tmp;
3197 }
3198 
3200 
3203 inline void accumulate(cdotprecision& dot, const srmatrix_subv& v1, const srmatrix_subv& v2) {
3204  dotprecision tmp(0.0);
3205  tmp.set_k(dot.get_k());
3206  accumulate(tmp,srvector(v1),srvector(v2));
3207  dot += tmp;
3208 }
3209 
3211 
3214 inline void accumulate(cdotprecision& dot, const srmatrix_subv& v1, const srvector& v2) {
3215  dotprecision tmp(0.0);
3216  tmp.set_k(dot.get_k());
3217  accumulate(tmp,srvector(v1),v2);
3218  dot += tmp;
3219 }
3220 
3222 
3225 inline void accumulate(cdotprecision& dot, const srmatrix_subv& v1, const srvector_slice& v2) {
3226  dotprecision tmp(0.0);
3227  tmp.set_k(dot.get_k());
3228  accumulate(tmp,srvector(v1),v2);
3229  dot += tmp;
3230 }
3231 
3233 
3236 inline void accumulate(cdotprecision& dot, const srmatrix_subv& v1, const rvector& v2) {
3237  dotprecision tmp(0.0);
3238  tmp.set_k(dot.get_k());
3239  accumulate(tmp,srvector(v1),v2);
3240  dot += tmp;
3241 }
3242 
3244 
3247 inline void accumulate(cdotprecision& dot, const srmatrix_subv& v1, const rvector_slice& v2) {
3248  dotprecision tmp(0.0);
3249  tmp.set_k(dot.get_k());
3250  accumulate(tmp,srvector(v1),v2);
3251  dot += tmp;
3252 }
3253 
3255 
3258 inline void accumulate(cdotprecision& dot, const srvector& v1, const srmatrix_subv& v2) {
3259  dotprecision tmp(0.0);
3260  tmp.set_k(dot.get_k());
3261  accumulate(tmp,v1,srvector(v2));
3262  dot += tmp;
3263 }
3264 
3266 
3269 inline void accumulate(cdotprecision& dot, const srvector_slice& v1, const srmatrix_subv& v2) {
3270  dotprecision tmp(0.0);
3271  tmp.set_k(dot.get_k());
3272  accumulate(tmp,v1,srvector(v2));
3273  dot += tmp;
3274 }
3275 
3277 
3280 inline void accumulate(cdotprecision& dot, const rvector& v1, const srmatrix_subv& v2) {
3281  dotprecision tmp(0.0);
3282  tmp.set_k(dot.get_k());
3283  accumulate(tmp,v1,srvector(v2));
3284  dot += tmp;
3285 }
3286 
3288 
3293 inline void accumulate_approx(cdotprecision& dot, const srmatrix_subv& v1, const srmatrix_subv& v2) {
3294  dotprecision tmp(0.0);
3295  tmp.set_k(dot.get_k());
3296  accumulate_approx(tmp,srvector(v1),srvector(v2));
3297  dot += tmp;
3298 }
3299 
3301 
3306 inline void accumulate_approx(cdotprecision& dot, const srmatrix_subv& v1, const srvector& v2) {
3307  dotprecision tmp(0.0);
3308  tmp.set_k(dot.get_k());
3309  accumulate_approx(tmp,srvector(v1),v2);
3310  dot += tmp;
3311 }
3312 
3314 
3319 inline void accumulate_approx(cdotprecision& dot, const srmatrix_subv& v1, const srvector_slice& v2) {
3320  dotprecision tmp(0.0);
3321  tmp.set_k(dot.get_k());
3322  accumulate_approx(tmp,srvector(v1),v2);
3323  dot += tmp;
3324 }
3325 
3327 
3332 inline void accumulate_approx(cdotprecision& dot, const srmatrix_subv& v1, const rvector& v2) {
3333  dotprecision tmp(0.0);
3334  tmp.set_k(dot.get_k());
3335  accumulate_approx(tmp,srvector(v1),v2);
3336  dot += tmp;
3337 }
3338 
3340 
3345 inline void accumulate_approx(cdotprecision& dot, const srmatrix_subv& v1, const rvector_slice& v2) {
3346  dotprecision tmp(0.0);
3347  tmp.set_k(dot.get_k());
3348  accumulate_approx(tmp,srvector(v1),v2);
3349  dot += tmp;
3350 }
3351 
3353 
3358 inline void accumulate_approx(cdotprecision& dot, const srvector& v1, const srmatrix_subv& v2) {
3359  dotprecision tmp(0.0);
3360  tmp.set_k(dot.get_k());
3361  accumulate_approx(tmp,v1,srvector(v2));
3362  dot += tmp;
3363 }
3364 
3366 
3371 inline void accumulate_approx(cdotprecision& dot, const srvector_slice& v1, const srmatrix_subv& v2) {
3372  dotprecision tmp(0.0);
3373  tmp.set_k(dot.get_k());
3374  accumulate_approx(tmp,v1,srvector(v2));
3375  dot += tmp;
3376 }
3377 
3379 
3384 inline void accumulate_approx(cdotprecision& dot, const rvector& v1, const srmatrix_subv& v2) {
3385  dotprecision tmp(0.0);
3386  tmp.set_k(dot.get_k());
3387  accumulate_approx(tmp,v1,srvector(v2));
3388  dot += tmp;
3389 }
3390 
3392 
3395 inline void accumulate(cdotprecision& dot, const rvector_slice& v1, const srmatrix_subv& v2) {
3396  dotprecision tmp(0.0);
3397  tmp.set_k(dot.get_k());
3398  accumulate(tmp,v1,srvector(v2));
3399  dot += tmp;
3400 }
3401 
3403 
3406 inline void accumulate(cidotprecision& dot, const srmatrix_subv& v1, const srmatrix_subv& v2) {
3407  dotprecision tmp(0.0);
3408  tmp.set_k(dot.get_k());
3409  accumulate(tmp,srvector(v1),srvector(v2));
3410  dot += tmp;
3411 }
3412 
3414 
3417 inline void accumulate(cidotprecision& dot, const srmatrix_subv& v1, const srvector& v2) {
3418  dotprecision tmp(0.0);
3419  tmp.set_k(dot.get_k());
3420  accumulate(tmp,srvector(v1),v2);
3421  dot += tmp;
3422 }
3423 
3425 
3428 inline void accumulate(cidotprecision& dot, const srmatrix_subv& v1, const srvector_slice& v2) {
3429  dotprecision tmp(0.0);
3430  tmp.set_k(dot.get_k());
3431  accumulate(tmp,srvector(v1),v2);
3432  dot += tmp;
3433 }
3434 
3436 
3439 inline void accumulate(cidotprecision& dot, const srmatrix_subv& v1, const rvector& v2) {
3440  dotprecision tmp(0.0);
3441  tmp.set_k(dot.get_k());
3442  accumulate(tmp,srvector(v1),v2);
3443  dot += tmp;
3444 }
3445 
3447 
3450 inline void accumulate(cidotprecision& dot, const srmatrix_subv& v1, const rvector_slice& v2) {
3451  dotprecision tmp(0.0);
3452  tmp.set_k(dot.get_k());
3453  accumulate(tmp,srvector(v1),v2);
3454  dot += tmp;
3455 }
3456 
3458 
3461 inline void accumulate(cidotprecision& dot, const srvector& v1, const srmatrix_subv& v2) {
3462  dotprecision tmp(0.0);
3463  tmp.set_k(dot.get_k());
3464  accumulate(tmp,v1,srvector(v2));
3465  dot += tmp;
3466 }
3467 
3469 
3472 inline void accumulate(cidotprecision& dot, const srvector_slice& v1, const srmatrix_subv& v2) {
3473  dotprecision tmp(0.0);
3474  tmp.set_k(dot.get_k());
3475  accumulate(tmp,v1,srvector(v2));
3476  dot += tmp;
3477 }
3478 
3480 
3483 inline void accumulate(cidotprecision& dot, const rvector& v1, const srmatrix_subv& v2) {
3484  dotprecision tmp(0.0);
3485  tmp.set_k(dot.get_k());
3486  accumulate(tmp,v1,srvector(v2));
3487  dot += tmp;
3488 }
3489 
3491 
3494 inline void accumulate(cidotprecision& dot, const rvector_slice& v1, const srmatrix_subv& v2) {
3495  dotprecision tmp(0.0);
3496  tmp.set_k(dot.get_k());
3497  accumulate(tmp,v1,srvector(v2));
3498  dot += tmp;
3499 }
3500 
3501 } //namespace cxsc;
3502 
3503 #include "sparsematrix.inl"
3504 
3505 #endif
rmatrix_subv & operator=(const rmatrix_subv &rv)
Implementation of standard assigning operator.
Definition: rmatrix.inl:322
friend srmatrix Re(const scmatrix &)
Returns the real part of the sparse matrix A.
Definition: scmatrix.hpp:1004
rmatrix & operator-=(const srmatrix &m)
Implementation of substraction and allocation operation.
Definition: srmatrix.hpp:1164
void set_k(unsigned int i)
Set precision for computation of dot products.
Definition: dot.hpp:131
srmatrix & operator*=(const rmatrix &B)
Multiply the sparse matrix by B and assign the result to it.
Definition: srmatrix.hpp:665
srmatrix & operator=(const rmatrix_slice &A)
Assigns a dense matrix slice to the sparse matrix. Only the non zero entries of the dense matrix slic...
Definition: srmatrix.hpp:483
The Data Type rmatrix_slice.
Definition: rmatrix.hpp:1442
srmatrix_slice & operator=(const rmatrix &C)
Assing C to the slice.
Definition: srmatrix.hpp:1434
srmatrix_subv & operator*=(const real &)
Assign the componentwise product of the subvector with a scalar to the subvector. ...
Definition: srmatrix.hpp:2608
srmatrix_slice & operator-=(const srmatrix &M)
Assigns the element wise difference of the sparse slice and M to the slice.
Definition: srmatrix.hpp:1516
friend void SetLb(srmatrix &, const int, const int)
Sets the lower index bound of the rows (i==ROW) or columns (i==COL) to j.
Definition: srmatrix.hpp:852
srmatrix(const int r, const int c)
Creates an empty matrix with r rows and c columns, pre-reserving space for 2*(r+c) elements...
Definition: srmatrix.hpp:127
srmatrix_slice & operator-=(const rmatrix_slice &M)
Assigns the element wise difference of the sparse slice and M to the slice.
Definition: srmatrix.hpp:1528
The Data Type idotprecision.
Definition: idot.hpp:47
friend void SetUb(srmatrix &, const int, const int)
Sets the upper index bound of the rows (i==ROW) or columns (i==COL) to j.
Definition: srmatrix.hpp:867
real density() const
Returns the density (the number of non-zeros divided by the number of elements) of the matrix...
Definition: srmatrix.hpp:625
The Data Type dotprecision.
Definition: dot.hpp:111
srvector()
Default constructor, creates an empty vector of size 0.
Definition: srvector.hpp:68
friend srmatrix Im(const scmatrix &)
Returns the imaginary part of the sparse matrix A.
Definition: scmatrix.hpp:1022
const std::vector< int > & column_pointers() const
Returns a constant reference to the vector containing the column pointers (the array) ...
Definition: srmatrix.hpp:105
srmatrix_slice & operator+=(const srmatrix &M)
Assigns the element wise sum of the sparse slice and M to the slice.
Definition: srmatrix.hpp:1492
The Data Type intmatrix.
Definition: intmatrix.hpp:313
real & element(const int i, const int j)
Returns a reference to the element (i,j) of the matrix.
Definition: srmatrix.hpp:1552
srmatrix & operator*=(const rmatrix_slice &B)
Multiply the sparse matrix by B and assign the result to it.
Definition: srmatrix.hpp:670
Represents a row or column vector of a sparse matrix.
Definition: simatrix.hpp:4383
friend srmatrix Sup(const simatrix &)
Returns the Supremum of the matrix A.
Definition: simatrix.hpp:1069
friend srmatrix absmin(const simatrix &)
Returns the componentwise minimum absolute value.
Definition: simatrix.hpp:1105
srmatrix_subv & operator=(const real &v)
Assigns v to all elements of the subvector.
Definition: srmatrix.hpp:2221
srmatrix & operator+=(const rmatrix &B)
Add B to the sparse matrix and assign the result to it.
Definition: srmatrix.hpp:635
srmatrix_subv & operator=(const srmatrix_subv &v)
Assigns a vector to a subvector.
Definition: srmatrix.hpp:2216
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition: cdot.cpp:29
civector operator*(const cimatrix_subv &rv, const cinterval &s)
Implementation of multiplication operation.
Definition: cimatrix.inl:731
const real operator()(int i, int j) const
Returns a copy of the element in row i and column j.
Definition: srmatrix.hpp:504
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)
Creates a sparse matrix out of three vectors (arrays) forming a matrix stored in triplet, compressed row or compressed column storage.
Definition: srmatrix.hpp:151
std::vector< int > & row_indices()
Returns a reference to the vector containing the row indices (the array)
Definition: srmatrix.hpp:95
int get_k() const
Get currently set precision for computation of dot products.
Definition: cdot.hpp:91
srmatrix(const int r, const int c, const int e)
Creates an empty matrix with r rows and c columns, pre-reserving space for e elements.
Definition: srmatrix.hpp:136
A sparse interval vector.
Definition: sivector.hpp:59
srmatrix & operator/=(const real &r)
Divide all elements of the sparse matrix by r and assign the result to it.
Definition: srmatrix.hpp:685
real & operator[](const int i)
Returns a reference to the i-th element of the subvector.
Definition: srmatrix.hpp:2179
srmatrix_slice & operator*=(const srmatrix &M)
Assigns the product of the sparse slice and M to the slice.
Definition: srmatrix.hpp:1456
int VecLen(const scimatrix_subv &S)
Returns the length of the subvector.
Definition: scimatrix.hpp:9966
srmatrix operator()(const intvector &pervec, const intvector &q)
Performs a row and column permutation using two permutation vectors.
Definition: srmatrix.hpp:561
real & element(int i, int j)
Returns a reference to the element (i,j) of the matrix.
Definition: srmatrix.hpp:525
srmatrix_subv operator[](const cxscmatrix_column &)
Returns a column of the matrix as a sparse subvector object.
Definition: srmatrix.hpp:2327
srmatrix_slice & operator+=(const rmatrix &M)
Assigns the element wise sum of the sparse slice and M to the slice.
Definition: srmatrix.hpp:1498
A sparse complex interval matrix.
Definition: scimatrix.hpp:71
void accumulate_approx(cdotprecision &dp, const cmatrix_subv &rv1, const cmatrix_subv &rv2)
The accurate scalar product of the last two arguments added to the value of the first argument (witho...
Definition: cmatrix.cpp:99
std::vector< int > & column_pointers()
Returns a reference to the vector containing the column pointers (the array)
Definition: srmatrix.hpp:90
A sparse real matrix.
Definition: srmatrix.hpp:77
srmatrix()
Standard constructor, creates an empty matrix of dimension 0x0.
Definition: srmatrix.hpp:120
srmatrix_subv & operator=(const rvector &v)
Assigns a vector to a subvector.
Definition: srmatrix.hpp:2236
srmatrix_slice & operator-=(const rmatrix &M)
Assigns the element wise difference of the sparse slice and M to the slice.
Definition: srmatrix.hpp:1522
friend int Ub(const srmatrix &, int)
Returns the upper index bound for the rows or columns of A.
Definition: srmatrix.hpp:894
friend srmatrix SupIm(const scimatrix &)
Returns the imaginary part of the supremum of the matrix A.
Definition: scimatrix.hpp:1581
rmatrix & operator+=(const srmatrix &m)
Implementation of addition and allocation operation.
Definition: srmatrix.hpp:1156
friend int RowLen(const srmatrix &)
Returns the number of columns of the matrix.
Definition: srmatrix.hpp:904
srmatrix(const rmatrix &A)
Creates a sparse matrix out of a dense matrix A. Only the non zero elements of A are stored explicitl...
Definition: srmatrix.hpp:367
rmatrix_slice & operator*=(const srmatrix &m)
Implementation of multiplication and allocation operation.
Definition: srmatrix.hpp:1176
A sparse complex interval vector.
Definition: scivector.hpp:62
srmatrix & operator*=(const real &r)
Multiply all elements of the sparse matrix by r and assign the result to it.
Definition: srmatrix.hpp:680
srmatrix & operator-=(const rmatrix &B)
Subtract B from the sparse matrix and assign the result to it.
Definition: srmatrix.hpp:650
srmatrix_slice & operator=(const srmatrix_slice &C)
Assing C to the slice.
Definition: srmatrix.hpp:1444
srmatrix_slice & operator+=(const srmatrix_slice &M)
Assigns the element wise sum of the sparse slice and M to the slice.
Definition: srmatrix.hpp:1486
The Data Type cidotprecision.
Definition: cidot.hpp:57
A sparse real vector.
Definition: srvector.hpp:58
friend srmatrix mid(const simatrix &)
Returns the midpoint matrix for A.
Definition: simatrix.hpp:1166
The Data Type rvector_slice.
Definition: rvector.hpp:1063
friend srmatrix transp(const srmatrix &)
Returns the transpose of A.
Definition: srmatrix.hpp:779
A sparse complex vector.
Definition: scvector.hpp:58
A slice of a sparse real interval matrix.
Definition: simatrix.hpp:2431
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)
Creates a sparse matrix out of three arrays forming a matrix stored in triplet, compressed row or com...
Definition: srmatrix.hpp:263
friend srmatrix Id(const srmatrix &)
Return a sparse unity matrix of the same dimension as A.
Definition: srmatrix.hpp:756
srmatrix_slice & operator*=(const real &r)
Assigns the component wise product of the sparse slice and r to the slice.
Definition: srmatrix.hpp:1474
The Data Type rmatrix_subv.
Definition: rmatrix.hpp:53
srmatrix & operator-=(const srmatrix &B)
Subtract B from the sparse matrix and assign the result to it.
Definition: srmatrix.hpp:660
srmatrix_subv & operator=(const srvector &v)
Assigns a vector to a subvector.
Definition: srmatrix.hpp:2226
The Data Type rvector.
Definition: rvector.hpp:57
std::vector< real > & values()
Returns a reference to the vector containing the stored values (the array)
Definition: srmatrix.hpp:100
srmatrix operator()(const intmatrix &P)
Performs a row permutation using the permutation matrix P. Faster than explicitly computing the produ...
Definition: srmatrix.hpp:619
friend int Lb(const srmatrix &, int)
Returns the lower index bound for the rows or columns of A.
Definition: srmatrix.hpp:881
The Data Type cmatrix.
Definition: cmatrix.hpp:513
int get_nnz() const
Returns the number of non-zero entries (including explicitly stored zeros).
Definition: srmatrix.hpp:630
A slice of a sparse real matrix.
Definition: srmatrix.hpp:1360
rmatrix_slice & operator+=(const srmatrix &m)
Implementation of multiplication and allocation operation.
Definition: srmatrix.hpp:1160
The Data Type cdotprecision.
Definition: cdot.hpp:60
A sparse complex matrix.
Definition: scmatrix.hpp:69
srmatrix_slice & operator=(const rmatrix_slice &C)
Assing C to the slice.
Definition: srmatrix.hpp:1439
A slice of a sparse complex matrix.
Definition: scmatrix.hpp:1956
void Resize(cimatrix &A)
Resizes the matrix.
Definition: cimatrix.inl:1211
srmatrix & operator=(const real &A)
Assigns a real value to all elements of the matrix (resulting in a dense matrix!) ...
Definition: srmatrix.hpp:473
srmatrix_subv & operator+=(const srvector &)
Assign the sum of the subvector with a vector to the subvector.
Definition: srmatrix.hpp:2618
Helper class for slices of sparse vectors.
Definition: sivector.hpp:1831
rmatrix_subv & operator-=(const real &c)
Implementation of subtraction and allocation operation.
Definition: rmatrix.inl:429
const std::vector< real > & values() const
Returns a constant reference to the vector containing the stored values (the array) ...
Definition: srmatrix.hpp:115
srmatrix_subv & operator=(const rvector_slice &v)
Assigns a vector to a subvector.
Definition: srmatrix.hpp:2241
rmatrix_slice & operator-=(const srmatrix &m)
Implementation of multiplication and allocation operation.
Definition: srmatrix.hpp:1168
friend srmatrix abs(const srmatrix &)
Returns the componentwise absolute value of A.
Definition: srmatrix.hpp:815
srmatrix & operator-=(const rmatrix_slice &B)
Subtract B from the sparse matrix and assign the result to it.
Definition: srmatrix.hpp:655
srmatrix & operator+=(const rmatrix_slice &B)
Add B to the sparse matrix and assign the result to it.
Definition: srmatrix.hpp:640
Represents a row or column vector of a sparse matrix.
Definition: srmatrix.hpp:2157
The Data Type rmatrix.
Definition: rmatrix.hpp:470
const std::vector< int > & row_indices() const
Returns a constant reference to the vector containing the row indices (the array) ...
Definition: srmatrix.hpp:110
friend srmatrix absmax(const simatrix &)
Returns the componentwise maximum absolute value.
Definition: simatrix.hpp:1123
srmatrix & operator*=(const srmatrix &B)
Multiply the sparse matrix by B and assign the result to it.
Definition: srmatrix.hpp:675
friend srmatrix InfRe(const scimatrix &)
Returns the real part of the infimum of the matrix A.
Definition: scimatrix.hpp:1527
friend srmatrix CompMat(const simatrix &)
Returns Ostroswkis comparison matrix for A.
Definition: simatrix.hpp:1141
rmatrix_subv & operator+=(const real &c)
Implementation of addition and allocation operation.
Definition: rmatrix.inl:428
srmatrix_subv operator[](const int)
Returns a row of the matrix.
Definition: srmatrix.hpp:2359
rmatrix_slice & operator=(const srmatrix &m)
Implementation of standard assigning operator.
Definition: srmatrix.hpp:1151
srmatrix_slice & operator=(const real &C)
Assing C to all elements of the slice.
Definition: srmatrix.hpp:1424
srmatrix_slice & operator*=(const rmatrix_slice &M)
Assigns the product of the sparse slice and M to the slice.
Definition: srmatrix.hpp:1468
srmatrix_slice & operator-=(const srmatrix_slice &M)
Assigns the element wise difference of the sparse slice and M to the slice.
Definition: srmatrix.hpp:1510
rmatrix & operator*=(const srmatrix &m)
Implementation of multiplication and allocation operation.
Definition: srmatrix.hpp:1172
friend srmatrix SupRe(const scimatrix &)
Returns the real part of the supremum of the matrix A.
Definition: scimatrix.hpp:1563
int get_k() const
Get currently set precision for computation of dot products.
Definition: cidot.hpp:89
int get_k() const
Get currently set precision for computation of dot products.
Definition: idot.hpp:86
srmatrix_subv & operator-=(const srvector &)
Assign the difference of the subvector with a vector to the subvector.
Definition: srmatrix.hpp:2638
srmatrix_subv & operator=(const srvector_slice &v)
Assigns a vector to a subvector.
Definition: srmatrix.hpp:2231
friend std::istream & operator>>(std::istream &, srmatrix_slice &)
Standard input operator for sparse matrix slice.
Definition: srmatrix.hpp:2143
srmatrix_slice & operator*=(const srmatrix_slice &M)
Assigns the product of the sparse slice and M to the slice.
Definition: srmatrix.hpp:1450
const real operator[](const int i) const
Returns a copy of the i-th element of the subvector.
Definition: srmatrix.hpp:2199
srmatrix_subv & operator/=(const real &)
Assign the componentwise division of the subvector with a scalar to the subvector.
Definition: srmatrix.hpp:2613
friend srmatrix diam(const simatrix &)
Returns the componentwise diameter of A.
Definition: simatrix.hpp:1184
srmatrix_slice & operator*=(const rmatrix &M)
Assigns the product of the sparse slice and M to the slice.
Definition: srmatrix.hpp:1462
void dropzeros()
Drops explicitly stored zeros from the data structure.
Definition: srmatrix.hpp:449
Helper class for slices of sparse vectors.
Definition: scivector.hpp:4063
const real operator()(const int i, const int j) const
Returns a copy of the element (i,j) of the matrix.
Definition: srmatrix.hpp:1538
srmatrix_slice & operator+=(const rmatrix_slice &M)
Assigns the element wise sum of the sparse slice and M to the slice.
Definition: srmatrix.hpp:1504
The Data Type imatrix.
Definition: imatrix.hpp:659
Represents a row or column vector of a sparse matrix.
Definition: scimatrix.hpp:9628
friend srmatrix InfIm(const scimatrix &)
Returns the imaginary part of the infimum of the matrix A.
Definition: scimatrix.hpp:1545
STORAGE_TYPE
Enumeration depicting the storage type of a sparse matrix (Triplet storage, Compressed column storage...
Definition: srmatrix.hpp:42
friend srmatrix Inf(const simatrix &)
Returns the Infimum of the matrix A.
Definition: simatrix.hpp:1051
Helper class for slices of sparse vectors.
Definition: srvector.hpp:868
friend int ColLen(const srmatrix &)
Returns the number of rows of the matrix.
Definition: srmatrix.hpp:909
srmatrix_slice & operator/=(const real &r)
Assigns the component wise division of the sparse slice and M to the slice.
Definition: srmatrix.hpp:1480
Represents a row or column vector of a sparse matrix.
Definition: scmatrix.hpp:3345
rmatrix()
Constructor of class rmatrix.
Definition: rmatrix.inl:33
rmatrix & operator=(const real &r)
Implementation of standard assigning operator.
Definition: rmatrix.inl:338
civector operator/(const cimatrix_subv &rv, const cinterval &s)
Implementation of division operation.
Definition: cimatrix.inl:730
A slice of a sparse complex interval matrix.
Definition: scimatrix.hpp:4918
srmatrix_slice & operator=(const srmatrix &C)
Assing C to the slice.
Definition: srmatrix.hpp:1429
srmatrix & operator=(const rmatrix &A)
Assigns a dense matrix to the sparse matrix. Only the non zero entries of the dense matrix are used...
Definition: srmatrix.hpp:478
The Scalar Type real.
Definition: real.hpp:113
srmatrix operator()(const intmatrix &P, const intmatrix &Q)
Performs row and column permutations using the two permutation matrices P and Q. Faster than explicit...
Definition: srmatrix.hpp:612
srmatrix & operator+=(const srmatrix &B)
Add B to the sparse matrix and assign the result to it.
Definition: srmatrix.hpp:645
void full(rmatrix &A) const
Creates a full matrix out of the sparse matrix and stores it in A. This should normally be done using...
Definition: srmatrix.hpp:434
The Data Type intvector.
Definition: intvector.hpp:51
srmatrix operator()(const intvector &pervec)
Performs a row permutation using a permutation vector.
Definition: srmatrix.hpp:588
A sparse interval matrix.
Definition: simatrix.hpp:69
The Data Type cimatrix.
Definition: cimatrix.hpp:907
srmatrix(const int ms, const int ns, const rmatrix &A)
Constructor for banded matrices.
Definition: srmatrix.hpp:394