C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
scmatrix.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: scmatrix.hpp,v 1.20 2014/01/30 17:23:48 cxsc Exp $ */
25 
26 #ifndef _CXSC_SCMATRIX_HPP_INCLUDED
27 #define _CXSC_SCMATRIX_HPP_INCLUDED
28 
29 #include <complex.hpp>
30 #include <cmatrix.hpp>
31 #include <scvector.hpp>
32 #include <cidot.hpp>
33 #include <vector>
34 #include <algorithm>
35 #include <iostream>
36 #include <sparsecdot.hpp>
37 #include <sparsematrix.hpp>
38 #include <srmatrix.hpp>
39 
40 namespace cxsc {
41 
42 //definiert in srmatrix.hpp
43 //enum STORAGE_TYPE{triplet,compressed_row,compressed_column};
44 
45 class scmatrix_slice;
46 class scmatrix_subv;
47 class scimatrix;
48 class scimatrix_slice;
49 class scimatrix_subv;
50 
51 
52 inline bool comp_pair_c(std::pair<int,complex> p1, std::pair<int,complex> p2) {
53  return p1.first < p2.first;
54 }
55 
57 
69 class scmatrix {
70 
71  private:
72  std::vector<int> p;
73  std::vector<int> ind;
74  std::vector<complex> x;
75  int m;
76  int n;
77  int lb1,ub1,lb2,ub2;
78 
79  public:
80 
82  std::vector<int>& column_pointers() {
83  return p;
84  }
85 
87  std::vector<int>& row_indices() {
88  return ind;
89  }
90 
92  std::vector<complex>& values() {
93  return x;
94  }
95 
97  const std::vector<int>& column_pointers() const {
98  return p;
99  }
100 
102  const std::vector<int>& row_indices() const {
103  return ind;
104  }
105 
107  const std::vector<complex>& values() const {
108  return x;
109  }
110 
113  p.push_back(0);
114  m = n = 0;
115  lb1 = lb2 = ub1 = ub2 = 0;
116  }
117 
119  scmatrix(const int r, const int c) : m(r),n(c),lb1(1),ub1(r),lb2(1),ub2(c) {
120  p = std::vector<int>((n>0) ? n+1 : 1, 0);
121  ind.reserve(2*(m+n));
122  x.reserve(2*(m+n));
123 
124  p[0] = 0;
125  }
126 
128  scmatrix(const int r, const int c, const int e) : m(r),n(c),lb1(1),ub1(r),lb2(1),ub2(c) {
129  p = std::vector<int>((n>0) ? n+1 : 1, 0);
130  ind.reserve(e);
131  x.reserve(e);
132 
133  p[0] = 0;
134  }
135 
137 
143  scmatrix(const int m, const int n, const int nnz, const intvector& rows, const intvector& cols, const cvector& values, const enum STORAGE_TYPE t = triplet) {
144  if(t == triplet) {
145  this->m = m;
146  this->n = n;
147  p = std::vector<int>(n+1,0);
148  ind.reserve(nnz);
149  x.reserve(nnz);
150  lb1 = lb2 = 1;
151  ub1 = m; ub2 = n;
152 
153  std::vector<triplet_store<complex> > work;
154  work.reserve(nnz);
155 
156  for(int k=0 ; k<nnz ; k++) {
157  work.push_back(triplet_store<complex>(rows[Lb(rows)+k],cols[Lb(cols)+k],values[Lb(values)+k]));
158  }
159 
160  sort(work.begin(), work.end());
161 
162  int i=0;
163 
164  for(int j=0 ; j<n ; j++) {
165 
166  while((unsigned int)i < work.size() && work[i].col == j ) {
167  ind.push_back(work[i].row);
168  x.push_back(work[i].val);
169  i++;
170  }
171 
172  p[j+1] = i;
173  }
174 
175  } else if(t == compressed_row) {
176 
177  this->m = m;
178  this->n = n;
179  p = std::vector<int>(n+1,0);
180  ind.reserve(nnz);
181  x.reserve(nnz);
182  lb1 = lb2 = 1;
183  ub1 = m; ub2 = n;
184 
185  for(int i=0 ; i<n+1 ; i++)
186  p[i] = rows[Lb(rows)+i];
187 
188  std::vector<triplet_store<complex> > work;
189  work.reserve(nnz);
190 
191  for(int j=0 ; j<n ; j++) {
192  for(int k=p[j] ; k<p[j+1] ; k++) {
193  work.push_back(triplet_store<complex>(j,cols[Lb(cols)+k],values[Lb(values)+k]));
194  }
195  }
196 
197  sort(work.begin(), work.end());
198 
199  int i=0;
200 
201  for(int j=0 ; j<n ; j++) {
202 
203  while((unsigned int)i < work.size() && work[i].col == j ) {
204  ind.push_back(work[i].row);
205  x.push_back(work[i].val);
206  i++;
207  }
208 
209  p[j+1] = i;
210  }
211 
212  } else if(t == compressed_column) {
213  this->m = m;
214  this->n = n;
215  p = std::vector<int>(n+1,0);
216  ind.reserve(nnz);
217  x.reserve(nnz);
218  lb1 = lb2 = 1;
219  ub1 = m; ub2 = n;
220 
221  for(int i=0 ; i<n+1 ; i++)
222  p[i] = rows[Lb(rows)+i];
223 
224  std::vector<std::pair<int,complex> > work;
225  work.reserve(n);
226 
227  for(int j=0 ; j<n ; j++) {
228  work.clear();
229 
230  for(int k=p[j] ; k<p[j+1] ; k++) {
231  work.push_back(std::make_pair(cols[Lb(cols)+k],values[Lb(values)+k]));
232  }
233 
234  std::sort(work.begin(),work.end(),comp_pair_c);
235 
236  for(unsigned int i=0 ; i<work.size() ; i++) {
237  ind.push_back(work[i].first);
238  x.push_back(work[i].second);
239  }
240  }
241 
242  }
243 
244  }
245 
247 
254  scmatrix(const int m, const int n, const int nnz, const int* rows, const int* cols, const complex* values, const enum STORAGE_TYPE t = triplet) {
255  if(t == triplet) {
256  this->m = m;
257  this->n = n;
258  p = std::vector<int>(n+1,0);
259  ind.reserve(nnz);
260  x.reserve(nnz);
261  lb1 = lb2 = 1;
262  ub1 = m; ub2 = n;
263 
264  std::vector<triplet_store<complex> > work;
265  work.reserve(nnz);
266 
267  for(int k=0 ; k<nnz ; k++) {
268  work.push_back(triplet_store<complex>(rows[k],cols[k],values[k]));
269  }
270 
271  sort(work.begin(), work.end());
272 
273  int i=0;
274 
275  for(int j=0 ; j<n ; j++) {
276 
277  while((unsigned int)i < work.size() && work[i].col == j ) {
278  ind.push_back(work[i].row);
279  x.push_back(work[i].val);
280  i++;
281  }
282 
283  p[j+1] = i;
284  }
285 
286  } else if(t == compressed_row) {
287 
288  this->m = m;
289  this->n = n;
290  p = std::vector<int>(n+1,0);
291  ind.reserve(nnz);
292  x.reserve(nnz);
293  lb1 = lb2 = 1;
294  ub1 = m; ub2 = n;
295 
296  for(int i=0 ; i<n+1 ; i++)
297  p[i] = rows[i];
298 
299  std::vector<triplet_store<complex> > work;
300  work.reserve(nnz);
301 
302  for(int j=0 ; j<n ; j++) {
303  for(int k=p[j] ; k<p[j+1] ; k++) {
304  work.push_back(triplet_store<complex>(j,cols[k],values[k]));
305  }
306  }
307 
308  sort(work.begin(), work.end());
309 
310  int i=0;
311 
312  for(int j=0 ; j<n ; j++) {
313 
314  while((unsigned int)i < work.size() && work[i].col == j ) {
315  ind.push_back(work[i].row);
316  x.push_back(work[i].val);
317  i++;
318  }
319 
320  p[j+1] = i;
321  }
322 
323  } else if(t == compressed_column) {
324  this->m = m;
325  this->n = n;
326  p = std::vector<int>(n+1,0);
327  ind.reserve(nnz);
328  x.reserve(nnz);
329  lb1 = lb2 = 1;
330  ub1 = m; ub2 = n;
331 
332  for(int i=0 ; i<n+1 ; i++)
333  p[i] = rows[i];
334 
335  std::vector<std::pair<int,complex> > work;
336  work.reserve(n);
337 
338  for(int j=0 ; j<n ; j++) {
339  work.clear();
340 
341  for(int k=p[j] ; k<p[j+1] ; k++) {
342  work.push_back(std::make_pair(cols[k],values[k]));
343  }
344 
345  std::sort(work.begin(),work.end(),comp_pair_c);
346 
347  for(unsigned int i=0 ; i<work.size() ; i++) {
348  ind.push_back(work[i].first);
349  x.push_back(work[i].second);
350  }
351  }
352 
353  }
354 
355  }
356 
357 
359  scmatrix(const srmatrix& A) : p(A.p), ind(A.ind), m(A.m), n(A.n), lb1(A.lb1), ub1(A.ub1), lb2(A.lb2), ub2(A.ub2) {
360  x.reserve(A.get_nnz());
361  for(unsigned int i=0 ; i<A.x.size() ; i++)
362  x.push_back(complex(A.x[i]));
363  }
364 
365 
367  scmatrix(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(complex(A[i+lb1][j+lb2]));
380  nnz++;
381  }
382  }
383 
384  p[j+1] = nnz;
385  }
386 
387  }
388 
390  scmatrix(const cmatrix& A) : m(ColLen(A)),n(RowLen(A)),lb1(Lb(A,1)),ub1(Ub(A,1)),lb2(Lb(A,2)),ub2(Ub(A,2)) {
391  p = std::vector<int>((n>0) ? n+1 : 1, 0);
392  ind.reserve((m*n*0.1 < 2*m) ? (int)(m*n*0.1) : 2*m);
393  x.reserve((m*n*0.1 < 2*m) ? (int)(m*n*0.1) : 2*m);
394 
395  p[0] = 0;
396  int nnz = 0;
397 
398  for(int j=0 ; j<n ; j++) {
399  for(int i=0 ; i<m ; i++) {
400  if(A[i+lb1][j+lb2] != 0.0) {
401  ind.push_back(i);
402  x.push_back(complex(A[i+lb1][j+lb2]));
403  nnz++;
404  }
405  }
406 
407  p[j+1] = nnz;
408  }
409 
410  }
411 
413 
416  scmatrix(const int ms, const int ns, const cmatrix& A) : m(ms), n(ns), lb1(1), ub1(ms), lb2(1), ub2(ns) {
417  //Banded matrix constructor
418  int nnz = RowLen(A)*ColLen(A);
419  p = std::vector<int>((n>0) ? n+1 : 1, 0);
420  ind.reserve(nnz);
421  x.reserve(nnz);
422 
423  std::vector<triplet_store<complex> > work;
424  work.reserve(nnz);
425 
426 
427  for(int i=0 ; i<ColLen(A) ; i++) {
428  for(int j=Lb(A,2) ; j<=Ub(A,2) ; j++) {
429  if(i+j >=0 && i+j < n) {
430  work.push_back(triplet_store<complex>(i,i+j,A[i+Lb(A,1)][j]));
431  }
432  }
433  }
434 
435  sort(work.begin(), work.end());
436 
437  int i=0;
438 
439  for(int j=0 ; j<n ; j++) {
440 
441  while((unsigned int)i < work.size() && work[i].col == j ) {
442  ind.push_back(work[i].row);
443  x.push_back(work[i].val);
444  i++;
445  }
446 
447  p[j+1] = i;
448  }
449 
450  }
451 
453  scmatrix(const srmatrix_slice&);
455  scmatrix(const scmatrix_slice&);
456 
458  void full(cmatrix& A) const {
459  A = cmatrix(lb1,ub1,lb2,ub2);
460  A = 0.0;
461  for(int j=0 ; j<n ; j++) {
462  for(int k=p[j] ; k<p[j+1] ; k++) {
463  A[ind[k]+lb1][j+lb2] = x[k];
464  }
465  }
466  }
467 
469 
473  void dropzeros() {
474  std::vector<int> pnew(n+1,0);
475  std::vector<int> indnew;
476  std::vector<complex> xnew;
477  int nnznew = 0;
478 
479  for(int j=0 ; j<n ; j++) {
480  for(int k=p[j] ; k<p[j+1] ; k++) {
481  if(x[k] != 0.0) {
482  xnew.push_back(x[k]);
483  indnew.push_back(ind[k]);
484  nnznew++;
485  }
486  }
487  pnew[j+1] = nnznew;
488  }
489 
490  p = pnew;
491  ind = indnew;
492  x = xnew;
493  }
494 
495 
497  scmatrix& operator=(const real& A) {
498  return sp_ms_assign<scmatrix,real,complex>(*this,A);
499  }
500 
503  return sp_ms_assign<scmatrix,complex,complex>(*this,A);
504  }
505 
508  return spf_mm_assign<scmatrix,rmatrix,complex>(*this,A);
509  }
510 
513  return spf_mm_assign<scmatrix,cmatrix,complex>(*this,A);
514  }
515 
518  return spf_mm_assign<scmatrix,rmatrix_slice,complex>(*this,A);
519  }
520 
523  return spf_mm_assign<scmatrix,cmatrix_slice,complex>(*this,A);
524  }
525 
528  m = A.m;
529  n = A.n;
530  p = A.p;
531  ind = A.ind;
532  x.clear();
533  x.reserve(A.get_nnz());
534  for(unsigned int i=0 ; i<A.x.size() ; i++)
535  x.push_back(complex(A.x[i]));
536  return *this;
537  }
538 
539  /* scmatrix& operator=(const scmatrix& A) {
540  p = A.p;
541  ind = A.ind;
542  x = A.x;
543  return *this;
544  } */
545 
550 
552 
558  const complex operator()(int i, int j) const {
559 #if(CXSC_INDEX_CHECK)
560  if(i<lb1 || i>ub1 || j<lb2 || j>ub2)
561  cxscthrow(ROW_OR_COL_NOT_IN_MAT("scmatrix::operator()(int, int)"));
562 #endif
563  complex r(0.0);
564  for(int k=p[j-lb2] ; k<p[j-lb2+1] && ind[k]<=i-lb1 ; k++) {
565  if(ind[k] == i-lb1) r = x[k];
566  }
567  return r;
568  }
569 
571 
579  complex& element(int i, int j) {
580 #if(CXSC_INDEX_CHECK)
581  if(i<lb1 || i>ub1 || j<lb2 || j>ub2)
582  cxscthrow(ROW_OR_COL_NOT_IN_MAT("scmatrix::element()(int, int)"));
583 #endif
584  int k;
585  for(k=p[j-lb2] ; k<p[j-lb2+1] && ind[k]<=i-lb1 ; k++) {
586  if(ind[k] == i-lb1) return x[k];
587  }
588 
589  //Nicht gefunden, Element muss angelegt werden, da Schreibzugriff moeglich
590  std::vector<int>::iterator ind_it = ind.begin() + k;
591  std::vector<complex>::iterator x_it = x.begin() + k;
592  ind.insert(ind_it, i-lb1);
593  x_it = x.insert(x_it, complex(0.0));
594  for(k=j-lb2+1 ; k<(int)p.size() ; k++)
595  p[k]++;
596 
597  return *x_it;
598  }
599 
601  scmatrix_subv operator[](const cxscmatrix_column&);
603  scmatrix_subv operator[](const int);
605  const scmatrix_subv operator[](const cxscmatrix_column&) const;
607  const scmatrix_subv operator[](const int) const;
608 
610  scmatrix_slice operator()(const int, const int , const int, const int);
612  const scmatrix_slice operator()(const int, const int , const int, const int) const;
613 
615  scmatrix operator()(const intvector& pervec, const intvector& q) {
616  scmatrix A(m,n,get_nnz());
617  intvector per = perminv(pervec);
618 
619  int nnz=0;
620  for(int k=0 ; k<n ; k++) {
621  A.p[k] = nnz;
622 
623  std::map<int,complex> work;
624  for(int j=p[q[Lb(q)+k]] ; j<p[q[Lb(q)+k]+1] ; j++)
625  work.insert(std::make_pair(per[Lb(per)+ind[j]], x[j]));
626 
627  for(std::map<int,complex>::iterator it = work.begin() ; it != work.end() ; it++) {
628  A.ind.push_back(it->first);
629  A.x.push_back(it->second);
630  }
631 
632  nnz += work.size();
633 
634  }
635 
636  A.p[n] = nnz;
637 
638  return A;
639  }
640 
642  scmatrix operator()(const intvector& pervec) {
643  scmatrix A(m,n,get_nnz());
644  intvector per = perminv(pervec);
645 
646  for(int k=0 ; k<n ; k++) {
647  A.p[k] = p[k];
648 
649  std::map<int,complex> work;
650  for(int j=p[k] ; j<p[k+1] ; j++)
651  work.insert(std::make_pair(per[Lb(per)+ind[j]], x[j]));
652 
653  for(std::map<int,complex>::iterator it = work.begin() ; it != work.end() ; it++) {
654  A.ind.push_back(it->first);
655  A.x.push_back(it->second);
656  }
657 
658  }
659 
660  A.p[n] = p[n];
661 
662  return A;
663  }
664 
666  scmatrix operator()(const intmatrix& P, const intmatrix& Q) {
667  intvector p = permvec(P);
668  intvector q = perminv(permvec(Q));
669  return (*this)(p,q);
670  }
671 
674  intvector p = permvec(P);
675  return (*this)(p);
676  }
677 
679  real density() const {
680  return p[n]/((double)m*n);
681  }
682 
684  int get_nnz() const {
685  return p[n];
686  }
687 
690  return spf_mm_addassign<scmatrix,rmatrix,cmatrix>(*this,B);
691  }
692 
695  return spf_mm_addassign<scmatrix,cmatrix,cmatrix>(*this,B);
696  }
697 
700  return spf_mm_addassign<scmatrix,rmatrix_slice,cmatrix>(*this,B);
701  }
702 
705  return spf_mm_addassign<scmatrix,cmatrix_slice,cmatrix>(*this,B);
706  }
707 
710  return spsp_mm_addassign<scmatrix,srmatrix,complex>(*this,B);
711  }
712 
715  return spsp_mm_addassign<scmatrix,scmatrix,complex>(*this,B);
716  }
717 
720  return spf_mm_subassign<scmatrix,rmatrix,cmatrix>(*this,B);
721  }
722 
725  return spf_mm_subassign<scmatrix,cmatrix,cmatrix>(*this,B);
726  }
727 
730  return spf_mm_subassign<scmatrix,rmatrix_slice,cmatrix>(*this,B);
731  }
732 
735  return spf_mm_subassign<scmatrix,cmatrix_slice,cmatrix>(*this,B);
736  }
737 
740  return spsp_mm_subassign<scmatrix,srmatrix,complex>(*this,B);
741  }
742 
745  return spsp_mm_subassign<scmatrix,scmatrix,complex>(*this,B);
746  }
747 
750  return spf_mm_multassign<scmatrix,cmatrix,sparse_cdot,cmatrix>(*this,B);
751  }
752 
755  return spf_mm_multassign<scmatrix,rmatrix,sparse_cdot,cmatrix>(*this,B);
756  }
757 
760  return spf_mm_multassign<scmatrix,rmatrix_slice,sparse_cdot,cmatrix>(*this,B);
761  }
762 
765  return spf_mm_multassign<scmatrix,cmatrix_slice,sparse_cdot,cmatrix>(*this,B);
766  }
767 
770  return spsp_mm_multassign<scmatrix,srmatrix,sparse_cdot,complex>(*this,B);
771  }
772 
775  return spsp_mm_multassign<scmatrix,scmatrix,sparse_cdot,complex>(*this,B);
776  }
777 
779  scmatrix& operator*=(const real& r) {
780  return sp_ms_multassign(*this,r);
781  }
782 
785  return sp_ms_multassign(*this,r);
786  }
787 
789  scmatrix& operator/=(const real& r) {
790  return sp_ms_divassign(*this,r);
791  }
792 
795  return sp_ms_divassign(*this,r);
796  }
797 
798  friend void SetLb(scmatrix&, const int, const int);
799  friend void SetUb(scmatrix&, const int, const int);
800  friend int Lb(const scmatrix&, int);
801  friend int Ub(const scmatrix&, int);
802  friend int RowLen(const scmatrix&);
803  friend int ColLen(const scmatrix&);
804  friend srmatrix Re(const scmatrix&);
805  friend srmatrix Im(const scmatrix&);
806  friend scmatrix Inf(const scimatrix&);
807  friend scmatrix Sup(const scimatrix&);
808  friend scmatrix mid(const scimatrix&);
809  friend scmatrix diam(const scimatrix&);
810  friend srmatrix abs(const scmatrix&);
811 
812  friend srmatrix CompMat(const scmatrix&);
813  friend scmatrix transp(const scmatrix&);
814  friend scmatrix Id(const scmatrix&);
815 
816  friend std::istream& operator>>(std::istream&, scmatrix_slice&);
817  friend std::istream& operator>>(std::istream&, scmatrix_subv&);
818 
819  friend class srmatrix_slice;
820  friend class srmatrix_subv;
821  friend class srvector;
822  friend class scmatrix_slice;
823  friend class scmatrix_subv;
824  friend class scvector;
825  friend class scivector;
826  friend class scimatrix;
827  friend class scimatrix_slice;
828  friend class scimatrix_subv;
829  friend class cmatrix;
830  friend class cimatrix;
831 
832 
833 #include "matrix_friend_declarations.inl"
834 };
835 
836 inline cmatrix::cmatrix(const srmatrix& A) {
837  dat = new complex[A.m*A.n];
838  lb1 = A.lb1; lb2 = A.lb2; ub1 = A.ub1; ub2 = A.ub2;
839  xsize = A.n;
840  ysize = A.m;
841  *this = 0.0;
842  for(int j=0 ; j<A.n ; j++) {
843  for(int k=A.p[j] ; k<A.p[j+1] ; k++) {
844  dat[A.ind[k]*A.n+j] = A.x[k];
845  }
846  }
847 }
848 
849 inline cmatrix::cmatrix(const scmatrix& A) {
850  dat = new complex[A.m*A.n];
851  lb1 = A.lb1; lb2 = A.lb2; ub1 = A.ub1; ub2 = A.ub2;
852  xsize = A.n;
853  ysize = A.m;
854  *this = 0.0;
855  for(int j=0 ; j<A.n ; j++) {
856  for(int k=A.p[j] ; k<A.p[j+1] ; k++) {
857  dat[A.ind[k]*A.n+j] = A.x[k];
858  }
859  }
860 }
861 
863 inline scmatrix Id(const scmatrix& A) {
864  scmatrix I(A.m, A.n, (A.m>A.n) ? A.m : A.n);
865  I.lb1 = A.lb1; I.lb2 = A.lb2;
866  I.ub1 = A.ub1; I.ub2 = A.ub2;
867 
868  if(A.m < A.n) {
869  for(int i=0 ; i<A.m ; i++) {
870  I.p[i+1] = I.p[i] + 1;
871  I.ind.push_back(i);
872  I.x.push_back(complex(1.0));
873  }
874  } else {
875  for(int i=0 ; i<A.n ; i++) {
876  I.p[i+1] = I.p[i] + 1;
877  I.ind.push_back(i);
878  I.x.push_back(complex(1.0));
879  }
880  }
881 
882  return I;
883 }
884 
886 inline scmatrix transp(const scmatrix& A) {
887  scmatrix B(A.n, A.m, A.get_nnz());
888 
889  //Nichtnullen pro Zeile bestimmen
890  std::vector<int> w(A.m,0);
891  for(unsigned int i=0 ; i<A.ind.size() ; i++)
892  w[A.ind[i]]++;
893 
894  //Spalten"pointer" setzen
895  B.p.resize(A.m+1);
896  B.p[0] = 0;
897  for(unsigned int i=1 ; i<B.p.size() ; i++)
898  B.p[i] = w[i-1] + B.p[i-1];
899 
900  //w vorbereiten
901  w.insert(w.begin(), 0);
902  for(unsigned int i=1 ; i<w.size() ; i++) {
903  w[i] += w[i-1];
904  }
905 
906  //neuer zeilenindex und wert wird gesetzt
907  int q;
908  B.ind.resize(A.get_nnz());
909  B.x.resize(A.get_nnz());
910  for(int j=0 ; j<A.n ; j++) {
911  for(int k=A.p[j] ; k<A.p[j+1] ; k++) {
912  q = w[A.ind[k]]++;
913  B.ind[q] = j;
914  B.x[q] = A.x[k];
915  }
916  }
917 
918  return B;
919 }
920 
922 
926 inline void SetLb(scmatrix& A, const int i, const int j) {
927  if(i==1) {
928  A.lb1 = j;
929  A.ub1 = j + A.m - 1;
930  } else if(i==2) {
931  A.lb2 = j;
932  A.ub2 = j + A.n - 1;
933  }
934 }
935 
937 
941 inline void SetUb(scmatrix& A, const int i, const int j) {
942  if(i==1) {
943  A.ub1 = j;
944  A.lb1 = j - A.m + 1;
945  } else if(i==2) {
946  A.ub2 = j;
947  A.lb2 = j - A.n + 1;
948  }
949 }
950 
951 
953 
956 inline int Lb(const scmatrix& A, int i) {
957  if(i==1)
958  return A.lb1;
959  else if(i==2)
960  return A.lb2;
961  else
962  return 1;
963 }
964 
966 
969 inline int Ub(const scmatrix& A, int i) {
970  if(i==1)
971  return A.ub1;
972  else if(i==2)
973  return A.ub2;
974  else
975  return 1;
976 }
977 
979 inline int RowLen(const scmatrix& A) {
980  return A.n;
981 }
982 
984 inline int ColLen(const scmatrix& A) {
985  return A.m;
986 }
987 
989 inline void Resize(scmatrix& A) {
990  sp_m_resize(A);
991 }
992 
994 inline void Resize(scmatrix& A, const int m, const int n) {
995  sp_m_resize(A,m,n);
996 }
997 
999 inline void Resize(scmatrix& A, const int l1, const int u1, const int l2, const int u2) {
1000  sp_m_resize(A,l1,u1,l2,u2);
1001 }
1002 
1004 inline srmatrix Re(const scmatrix& A) {
1005  srmatrix res(A.m,A.n,A.get_nnz());
1006  res.lb1 = A.lb1;
1007  res.lb2 = A.lb2;
1008  res.ub1 = A.ub1;
1009  res.ub2 = A.ub2;
1010  res.p = A.p;
1011  res.ind = A.ind;
1012 
1013  for(int i=0 ; i<res.get_nnz() ; i++)
1014  res.x.push_back(Re(A.x[i]));
1015 
1016  res.dropzeros();
1017 
1018  return res;
1019 }
1020 
1022 inline srmatrix Im(const scmatrix& A) {
1023  srmatrix res(A.m,A.n,A.get_nnz());
1024  res.lb1 = A.lb1;
1025  res.lb2 = A.lb2;
1026  res.ub1 = A.ub1;
1027  res.ub2 = A.ub2;
1028  res.p = A.p;
1029  res.ind = A.ind;
1030 
1031  for(int i=0 ; i<res.get_nnz() ; i++)
1032  res.x.push_back(Im(A.x[i]));
1033 
1034  res.dropzeros();
1035 
1036  return res;
1037 }
1038 
1040 inline srmatrix abs(const scmatrix& A) {
1041  srmatrix ret;
1042  ret.ind = A.ind;
1043  ret.p = A.p;
1044  for(unsigned int i=0 ; i<ret.x.size() ; i++)
1045  ret.x[i] = abs(A.x[i]);
1046  return ret;
1047 }
1048 
1050 inline srmatrix CompMat(const scmatrix& A) {
1051  srmatrix res(A.m,A.n,A.get_nnz());
1052  res.lb1 = A.lb1;
1053  res.lb2 = A.lb2;
1054  res.ub1 = A.ub1;
1055  res.ub2 = A.ub2;
1056  res.p = A.p;
1057  res.ind = A.ind;
1058  res.p[A.n] = A.p[A.n];
1059 
1060  for(int j=0 ; j<res.n ; j++) {
1061  for(int k=A.p[j] ; k<A.p[j+1] ; k++) {
1062  if(A.ind[k] == j)
1063  res.x.push_back(abs(A.x[k]));
1064  else
1065  res.x.push_back(-abs(A.x[k]));
1066  }
1067  }
1068 
1069  res.dropzeros();
1070 
1071  return res;
1072 }
1073 
1075 
1081 inline cmatrix operator*(const cmatrix& A, const srmatrix& B) {
1082  return fsp_mm_mult<cmatrix,srmatrix,cmatrix,sparse_cdot>(A,B);
1083 }
1084 
1086 
1092 inline cmatrix operator*(const rmatrix& A, const scmatrix& B) {
1093  return fsp_mm_mult<rmatrix,scmatrix,cmatrix,sparse_cdot>(A,B);
1094 }
1095 
1097 
1103 inline cmatrix operator*(const cmatrix& A, const scmatrix& B) {
1104  return fsp_mm_mult<cmatrix,scmatrix,cmatrix,sparse_cdot>(A,B);
1105 }
1106 
1108 
1114 inline cmatrix operator*(const scmatrix& A, const rmatrix& B) {
1115  return spf_mm_mult<scmatrix,rmatrix,cmatrix,sparse_cdot>(A,B);
1116 }
1117 
1119 
1125 inline cmatrix operator*(const srmatrix& A, const cmatrix& B) {
1126  return spf_mm_mult<srmatrix,cmatrix,cmatrix,sparse_cdot>(A,B);
1127 }
1128 
1130 
1136 inline cmatrix operator*(const scmatrix& A, const cmatrix& B) {
1137  return spf_mm_mult<scmatrix,cmatrix,cmatrix,sparse_cdot>(A,B);
1138 }
1139 
1141 
1147 inline cmatrix operator*(const cmatrix_slice& A, const srmatrix& B) {
1148  return fsp_mm_mult<cmatrix_slice,srmatrix,cmatrix,sparse_cdot>(A,B);
1149 }
1150 
1152 
1158 inline cmatrix operator*(const rmatrix_slice& A, const scmatrix& B) {
1159  return fsp_mm_mult<rmatrix_slice,scmatrix,cmatrix,sparse_cdot>(A,B);
1160 }
1161 
1163 
1169 inline cmatrix operator*(const cmatrix_slice& A, const scmatrix& B) {
1170  return fsp_mm_mult<cmatrix_slice,scmatrix,cmatrix,sparse_cdot>(A,B);
1171 }
1172 
1174 
1180 inline cmatrix operator*(const scmatrix& A, const rmatrix_slice& B) {
1181  return spf_mm_mult<scmatrix,rmatrix_slice,cmatrix,sparse_cdot>(A,B);
1182 }
1183 
1185 
1191 inline cmatrix operator*(const srmatrix& A, const cmatrix_slice& B) {
1192  return spf_mm_mult<srmatrix,cmatrix_slice,cmatrix,sparse_cdot>(A,B);
1193 }
1194 
1196 
1202 inline cmatrix operator*(const scmatrix& A, const cmatrix_slice& B) {
1203  return spf_mm_mult<scmatrix,cmatrix_slice,cmatrix,sparse_cdot>(A,B);
1204 }
1205 
1207 
1213 inline scmatrix operator*(const scmatrix& A, const srmatrix& B) {
1214  return spsp_mm_mult<scmatrix,srmatrix,scmatrix,sparse_cdot,complex>(A,B);
1215 }
1216 
1218 
1224 inline scmatrix operator*(const srmatrix& A, const scmatrix& B) {
1225  return spsp_mm_mult<srmatrix,scmatrix,scmatrix,sparse_cdot,complex>(A,B);
1226 }
1227 
1229 
1235 inline scmatrix operator*(const scmatrix& A, const scmatrix& B) {
1236  return spsp_mm_mult<scmatrix,scmatrix,scmatrix,sparse_cdot,complex>(A,B);
1237 }
1238 
1240 inline scmatrix operator/(const scmatrix& A, const real& r) {
1241  return sp_ms_div<scmatrix,real,scmatrix>(A,r);
1242 }
1243 
1245 inline scmatrix operator/(const scmatrix& A, const complex& r) {
1246  return sp_ms_div<scmatrix,complex,scmatrix>(A,r);
1247 }
1248 
1250 inline scmatrix operator/(const srmatrix& A, const complex& r) {
1251  return sp_ms_div<srmatrix,complex,scmatrix>(A,r);
1252 }
1253 
1255 inline scmatrix operator*(const scmatrix& A, const real& r) {
1256  return sp_ms_mult<scmatrix,real,scmatrix>(A,r);
1257 }
1258 
1260 inline scmatrix operator*(const scmatrix& A, const complex& r) {
1261  return sp_ms_mult<scmatrix,complex,scmatrix>(A,r);
1262 }
1263 
1265 inline scmatrix operator*(const srmatrix& A, const complex& r) {
1266  return sp_ms_mult<srmatrix,complex,scmatrix>(A,r);
1267 }
1268 
1270 inline scmatrix operator*(const real& r, const scmatrix& A) {
1271  return sp_sm_mult<real,scmatrix,scmatrix>(r,A);
1272 }
1273 
1275 inline scmatrix operator*(const complex& r, const scmatrix& A) {
1276  return sp_sm_mult<complex,scmatrix,scmatrix>(r,A);
1277 }
1278 
1280 inline scmatrix operator*(const complex& r, const srmatrix& A) {
1281  return sp_sm_mult<complex,srmatrix,scmatrix>(r,A);
1282 }
1283 
1285 
1291 inline cvector operator*(const scmatrix& A, const rvector& v) {
1292  return spf_mv_mult<scmatrix,rvector,cvector,sparse_cdot>(A,v);
1293 }
1294 
1296 
1302 inline cvector operator*(const srmatrix& A, const cvector& v) {
1303  return spf_mv_mult<srmatrix,cvector,cvector,sparse_cdot>(A,v);
1304 }
1305 
1307 
1313 inline cvector operator*(const scmatrix& A, const cvector& v) {
1314  return spf_mv_mult<scmatrix,cvector,cvector,sparse_cdot>(A,v);
1315 }
1316 
1318 
1324 inline cvector operator*(const scmatrix& A, const rvector_slice& v) {
1325  return spf_mv_mult<scmatrix,rvector_slice,cvector,sparse_cdot>(A,v);
1326 }
1327 
1329 
1335 inline cvector operator*(const srmatrix& A, const cvector_slice& v) {
1336  return spf_mv_mult<srmatrix,cvector_slice,cvector,sparse_cdot>(A,v);
1337 }
1338 
1340 
1346 inline cvector operator*(const scmatrix& A, const cvector_slice& v) {
1347  return spf_mv_mult<scmatrix,cvector_slice,cvector,sparse_cdot>(A,v);
1348 }
1349 
1351 
1357 inline scvector operator*(const scmatrix& A, const srvector& v) {
1358  return spsp_mv_mult<scmatrix,srvector,scvector,sparse_cdot,complex>(A,v);
1359 }
1360 
1362 
1368 inline scvector operator*(const srmatrix& A, const scvector& v) {
1369  return spsp_mv_mult<srmatrix,scvector,scvector,sparse_cdot,complex>(A,v);
1370 }
1371 
1373 
1379 inline scvector operator*(const scmatrix& A, const scvector& v) {
1380  return spsp_mv_mult<scmatrix,scvector,scvector,sparse_cdot,complex>(A,v);
1381 }
1382 
1384 
1390 inline scvector operator*(const scmatrix& A, const srvector_slice& v) {
1391  return spsl_mv_mult<scmatrix,srvector_slice,scvector,sparse_cdot,complex>(A,v);
1392 }
1393 
1395 
1401 inline scvector operator*(const srmatrix& A, const scvector_slice& v) {
1402  return spsl_mv_mult<srmatrix,scvector_slice,scvector,sparse_cdot,complex>(A,v);
1403 }
1404 
1406 
1412 inline scvector operator*(const scmatrix& A, const scvector_slice& v) {
1413  return spsl_mv_mult<scmatrix,scvector_slice,scvector,sparse_cdot,complex>(A,v);
1414 }
1415 
1417 
1423 inline cvector operator*(const cmatrix& A, const srvector& v) {
1424  return fsp_mv_mult<cmatrix,srvector,cvector,sparse_cdot>(A,v);
1425 }
1426 
1428 
1434 inline cvector operator*(const rmatrix& A, const scvector& v) {
1435  return fsp_mv_mult<rmatrix,scvector,cvector,sparse_cdot>(A,v);
1436 }
1437 
1439 
1445 inline cvector operator*(const cmatrix& A, const scvector& v) {
1446  return fsp_mv_mult<cmatrix,scvector,cvector,sparse_cdot>(A,v);
1447 }
1448 
1450 
1456 inline cvector operator*(const cmatrix_slice& A, const srvector& v) {
1457  return fsp_mv_mult<cmatrix_slice,srvector,cvector,sparse_cdot>(A,v);
1458 }
1459 
1461 
1467 inline cvector operator*(const rmatrix_slice& A, const scvector& v) {
1468  return fsp_mv_mult<rmatrix_slice,scvector,cvector,sparse_cdot>(A,v);
1469 }
1470 
1472 
1478 inline cvector operator*(const cmatrix_slice& A, const scvector& v) {
1479  return fsp_mv_mult<cmatrix_slice,scvector,cvector,sparse_cdot>(A,v);
1480 }
1481 
1483 
1489 inline cvector operator*(const cmatrix& A, const srvector_slice& v) {
1490  return fsl_mv_mult<cmatrix,srvector_slice,cvector,sparse_cdot>(A,v);
1491 }
1492 
1494 
1500 inline cvector operator*(const rmatrix& A, const scvector_slice& v) {
1501  return fsl_mv_mult<rmatrix,scvector_slice,cvector,sparse_cdot>(A,v);
1502 }
1503 
1505 
1511 inline cvector operator*(const cmatrix& A, const scvector_slice& v) {
1512  return fsl_mv_mult<cmatrix,scvector_slice,cvector,sparse_cdot>(A,v);
1513 }
1514 
1516 
1522 inline cvector operator*(const cmatrix_slice& A, const srvector_slice& v) {
1523  return fsl_mv_mult<cmatrix_slice,srvector_slice,cvector,sparse_cdot>(A,v);
1524 }
1525 
1527 
1533 inline cvector operator*(const rmatrix_slice& A, const scvector_slice& v) {
1534  return fsl_mv_mult<rmatrix_slice,scvector_slice,cvector,sparse_cdot>(A,v);
1535 }
1536 
1538 
1544 inline cvector operator*(const cmatrix_slice& A, const scvector_slice& v) {
1545  return fsl_mv_mult<cmatrix_slice,scvector_slice,cvector,sparse_cdot>(A,v);
1546 }
1547 
1549 inline cmatrix operator+(const cmatrix& A, const srmatrix& B) {
1550  return fsp_mm_add<cmatrix,srmatrix,cmatrix>(A,B);
1551 }
1552 
1554 inline cmatrix operator+(const rmatrix& A, const scmatrix& B) {
1555  return fsp_mm_add<rmatrix,scmatrix,cmatrix>(A,B);
1556 }
1557 
1559 inline cmatrix operator+(const cmatrix& A, const scmatrix& B) {
1560  return fsp_mm_add<cmatrix,scmatrix,cmatrix>(A,B);
1561 }
1562 
1564 inline cmatrix operator+(const scmatrix& A, const rmatrix& B) {
1565  return spf_mm_add<scmatrix,rmatrix,cmatrix>(A,B);
1566 }
1567 
1569 inline cmatrix operator+(const srmatrix& A, const cmatrix& B) {
1570  return spf_mm_add<srmatrix,cmatrix,cmatrix>(A,B);
1571 }
1572 
1574 inline cmatrix operator+(const scmatrix& A, const cmatrix& B) {
1575  return spf_mm_add<scmatrix,cmatrix,cmatrix>(A,B);
1576 }
1577 
1579 inline cmatrix operator+(const cmatrix_slice& A, const srmatrix& B) {
1580  return fsp_mm_add<cmatrix_slice,srmatrix,cmatrix>(A,B);
1581 }
1582 
1584 inline cmatrix operator+(const rmatrix_slice& A, const scmatrix& B) {
1585  return fsp_mm_add<rmatrix_slice,scmatrix,cmatrix>(A,B);
1586 }
1587 
1589 inline cmatrix operator+(const cmatrix_slice& A, const scmatrix& B) {
1590  return fsp_mm_add<cmatrix_slice,scmatrix,cmatrix>(A,B);
1591 }
1592 
1594 inline cmatrix operator+(const scmatrix& A, const rmatrix_slice& B) {
1595  return spf_mm_add<scmatrix,rmatrix_slice,cmatrix>(A,B);
1596 }
1597 
1599 inline cmatrix operator+(const srmatrix& A, const cmatrix_slice& B) {
1600  return spf_mm_add<srmatrix,cmatrix_slice,cmatrix>(A,B);
1601 }
1602 
1604 inline cmatrix operator+(const scmatrix& A, const cmatrix_slice& B) {
1605  return spf_mm_add<scmatrix,cmatrix_slice,cmatrix>(A,B);
1606 }
1607 
1609 inline scmatrix operator+(const scmatrix& A, const srmatrix& B) {
1610  return spsp_mm_add<scmatrix,srmatrix,scmatrix,complex>(A,B);
1611 }
1612 
1614 inline scmatrix operator+(const srmatrix& A, const scmatrix& B) {
1615  return spsp_mm_add<srmatrix,scmatrix,scmatrix,complex>(A,B);
1616 }
1617 
1619 inline scmatrix operator+(const scmatrix& A, const scmatrix& B) {
1620  return spsp_mm_add<scmatrix,scmatrix,scmatrix,complex>(A,B);
1621 }
1622 
1624 inline cmatrix operator-(const cmatrix& A, const srmatrix& B) {
1625  return fsp_mm_sub<cmatrix,srmatrix,cmatrix>(A,B);
1626 }
1627 
1629 inline cmatrix operator-(const rmatrix& A, const scmatrix& B) {
1630  return fsp_mm_sub<rmatrix,scmatrix,cmatrix>(A,B);
1631 }
1632 
1634 inline cmatrix operator-(const cmatrix& A, const scmatrix& B) {
1635  return fsp_mm_sub<cmatrix,scmatrix,cmatrix>(A,B);
1636 }
1637 
1639 inline cmatrix operator-(const scmatrix& A, const rmatrix& B) {
1640  return spf_mm_sub<scmatrix,rmatrix,cmatrix>(A,B);
1641 }
1642 
1644 inline cmatrix operator-(const srmatrix& A, const cmatrix& B) {
1645  return spf_mm_sub<srmatrix,cmatrix,cmatrix>(A,B);
1646 }
1647 
1649 inline cmatrix operator-(const scmatrix& A, const cmatrix& B) {
1650  return spf_mm_sub<scmatrix,cmatrix,cmatrix>(A,B);
1651 }
1652 
1654 inline cmatrix operator-(const cmatrix_slice& A, const srmatrix& B) {
1655  return fsp_mm_sub<cmatrix_slice,srmatrix,cmatrix>(A,B);
1656 }
1657 
1659 inline cmatrix operator-(const rmatrix_slice& A, const scmatrix& B) {
1660  return fsp_mm_sub<rmatrix_slice,scmatrix,cmatrix>(A,B);
1661 }
1662 
1664 inline cmatrix operator-(const cmatrix_slice& A, const scmatrix& B) {
1665  return fsp_mm_sub<cmatrix_slice,scmatrix,cmatrix>(A,B);
1666 }
1667 
1669 inline cmatrix operator-(const scmatrix& A, const rmatrix_slice& B) {
1670  return spf_mm_sub<scmatrix,rmatrix_slice,cmatrix>(A,B);
1671 }
1672 
1674 inline cmatrix operator-(const srmatrix& A, const cmatrix_slice& B) {
1675  return spf_mm_sub<srmatrix,cmatrix_slice,cmatrix>(A,B);
1676 }
1677 
1679 inline cmatrix operator-(const scmatrix& A, const cmatrix_slice& B) {
1680  return spf_mm_sub<scmatrix,cmatrix_slice,cmatrix>(A,B);
1681 }
1682 
1684 inline scmatrix operator-(const scmatrix& A, const srmatrix& B) {
1685  return spsp_mm_sub<scmatrix,srmatrix,scmatrix,complex>(A,B);
1686 }
1687 
1689 inline scmatrix operator-(const srmatrix& A, const scmatrix& B) {
1690  return spsp_mm_sub<srmatrix,scmatrix,scmatrix,complex>(A,B);
1691 }
1692 
1694 inline scmatrix operator-(const scmatrix& A, const scmatrix& B) {
1695  return spsp_mm_sub<scmatrix,scmatrix,scmatrix,complex>(A,B);
1696 }
1697 
1699 inline scmatrix operator-(const scmatrix& M) {
1700  return sp_m_negative<scmatrix,scmatrix>(M);
1701 }
1702 
1704 inline scmatrix& operator+(scmatrix& A) {
1705  return A;
1706 }
1707 
1709  *this = rmatrix(B);
1710  return *this;
1711 }
1712 
1714  *this = rmatrix(B);
1715  return *this;
1716 }
1717 
1719  *this = cmatrix(B);
1720  return *this;
1721 }
1722 
1724  *this = cmatrix(B);
1725  return *this;
1726 }
1727 
1729  return fsp_mm_addassign(*this,B);
1730 }
1731 
1733  return fsp_mm_addassign(*this,B);
1734 }
1735 
1737  return fsp_mm_addassign(*this,B);
1738 }
1739 
1741  return fsp_mm_addassign(*this,B);
1742 }
1743 
1744 inline cmatrix& cmatrix::operator-=(const srmatrix& B) {
1745  return fsp_mm_subassign(*this,B);
1746 }
1747 
1748 inline cmatrix& cmatrix::operator-=(const scmatrix& B) {
1749  return fsp_mm_subassign(*this,B);
1750 }
1751 
1752 inline cmatrix_slice& cmatrix_slice::operator-=(const srmatrix& B) {
1753  return fsp_mm_subassign(*this,B);
1754 }
1755 
1756 inline cmatrix_slice& cmatrix_slice::operator-=(const scmatrix& B) {
1757  return fsp_mm_subassign(*this,B);
1758 }
1759 
1761  return fsp_mm_multassign<cmatrix,srmatrix,sparse_cdot,cmatrix>(*this,B);
1762 }
1763 
1765  return fsp_mm_multassign<cmatrix,scmatrix,sparse_cdot,cmatrix>(*this,B);
1766 }
1767 
1769  return fsp_mm_multassign<cmatrix_slice,srmatrix,sparse_cdot,cmatrix>(*this,B);
1770 }
1771 
1773  return fsp_mm_multassign<cmatrix_slice,scmatrix,sparse_cdot,cmatrix>(*this,B);
1774 }
1775 
1777 inline bool operator==(const scmatrix& A, const srmatrix& B) {
1778  return spsp_mm_comp(A,B);
1779 }
1780 
1782 inline bool operator==(const srmatrix& A, const scmatrix& B) {
1783  return spsp_mm_comp(A,B);
1784 }
1785 
1787 inline bool operator==(const scmatrix& A, const scmatrix& B) {
1788  return spsp_mm_comp(A,B);
1789 }
1790 
1792 inline bool operator==(const scmatrix& A, const rmatrix& B) {
1793  return spf_mm_comp(A,B);
1794 }
1795 
1797 inline bool operator==(const srmatrix& A, const cmatrix& B) {
1798  return spf_mm_comp(A,B);
1799 }
1800 
1802 inline bool operator==(const scmatrix& A, const cmatrix& B) {
1803  return spf_mm_comp(A,B);
1804 }
1805 
1807 inline bool operator==(const cmatrix& A, const srmatrix& B) {
1808  return fsp_mm_comp(A,B);
1809 }
1810 
1812 inline bool operator==(const rmatrix& A, const scmatrix& B) {
1813  return fsp_mm_comp(A,B);
1814 }
1815 
1817 inline bool operator==(const cmatrix& A, const scmatrix& B) {
1818  return fsp_mm_comp(A,B);
1819 }
1820 
1822 inline bool operator==(const cmatrix_slice& A, const srmatrix& B) {
1823  return fsp_mm_comp(A,B);
1824 }
1825 
1827 inline bool operator==(const rmatrix_slice& A, const scmatrix& B) {
1828  return fsp_mm_comp(A,B);
1829 }
1830 
1832 inline bool operator==(const cmatrix_slice& A, const scmatrix& B) {
1833  return fsp_mm_comp(A,B);
1834 }
1835 
1837 inline bool operator==(const scmatrix& A, const rmatrix_slice& B) {
1838  return spf_mm_comp(A,B);
1839 }
1840 
1842 inline bool operator==(const srmatrix& A, const cmatrix_slice& B) {
1843  return spf_mm_comp(A,B);
1844 }
1845 
1847 inline bool operator==(const scmatrix& A, const cmatrix_slice& B) {
1848  return spf_mm_comp(A,B);
1849 }
1850 
1852 inline bool operator!=(const scmatrix& A, const srmatrix& B) {
1853  return !spsp_mm_comp(A,B);
1854 }
1855 
1857 inline bool operator!=(const srmatrix& A, const scmatrix& B) {
1858  return !spsp_mm_comp(A,B);
1859 }
1860 
1862 inline bool operator!=(const scmatrix& A, const scmatrix& B) {
1863  return !spsp_mm_comp(A,B);
1864 }
1865 
1867 inline bool operator!=(const scmatrix& A, const rmatrix& B) {
1868  return !spf_mm_comp(A,B);
1869 }
1870 
1872 inline bool operator!=(const srmatrix& A, const cmatrix& B) {
1873  return !spf_mm_comp(A,B);
1874 }
1875 
1877 inline bool operator!=(const scmatrix& A, const cmatrix& B) {
1878  return !spf_mm_comp(A,B);
1879 }
1880 
1882 inline bool operator!=(const cmatrix& A, const srmatrix& B) {
1883  return !fsp_mm_comp(A,B);
1884 }
1885 
1887 inline bool operator!=(const rmatrix& A, const scmatrix& B) {
1888  return !fsp_mm_comp(A,B);
1889 }
1890 
1892 inline bool operator!=(const cmatrix& A, const scmatrix& B) {
1893  return !fsp_mm_comp(A,B);
1894 }
1895 
1897 inline bool operator!=(const cmatrix_slice& A, const srmatrix& B) {
1898  return !fsp_mm_comp(A,B);
1899 }
1900 
1902 inline bool operator!=(const rmatrix_slice& A, const scmatrix& B) {
1903  return !fsp_mm_comp(A,B);
1904 }
1905 
1907 inline bool operator!=(const cmatrix_slice& A, const scmatrix& B) {
1908  return !fsp_mm_comp(A,B);
1909 }
1910 
1912 inline bool operator!=(const scmatrix& A, const rmatrix_slice& B) {
1913  return !spf_mm_comp(A,B);
1914 }
1915 
1917 inline bool operator!=(const srmatrix& A, const cmatrix_slice& B) {
1918  return !spf_mm_comp(A,B);
1919 }
1920 
1922 inline bool operator!=(const scmatrix& A, const cmatrix_slice& B) {
1923  return !spf_mm_comp(A,B);
1924 }
1925 
1927 inline bool operator!(const scmatrix& A) {
1928  return sp_m_not(A);
1929 }
1930 
1932 
1937 inline std::ostream& operator<<(std::ostream& os, const scmatrix& A) {
1938  return sp_m_output<scmatrix,complex>(os,A);
1939 }
1940 
1942 
1947 inline std::istream& operator>>(std::istream& is, scmatrix& A) {
1948  return sp_m_input<scmatrix,complex>(is,A);
1949 }
1950 
1952 
1957  public:
1958  scmatrix A;
1959  scmatrix* M; //Originalmatrix
1960 
1961  private:
1962  scmatrix_slice(scmatrix& Mat, int sl1l, int sl1u, int sl2l, int sl2u) {
1963  A.lb1 = sl1l;
1964  A.lb2 = sl2l;
1965  A.ub1 = sl1u;
1966  A.ub2 = sl2u;
1967  A.m = sl1u-sl1l+1;
1968  A.n = sl2u-sl2l+1;
1969 
1970  //Kopieren der Werte aus A
1971  A.p = std::vector<int>(A.n+1, 0);
1972  A.ind.reserve(A.m + A.n);
1973  A.x.reserve(A.m + A.n);
1974 
1975  for(int i=0 ; i<A.n ; i++) {
1976  A.p[i+1] = A.p[i];
1977  for(int j=Mat.p[sl2l-Mat.lb2+i] ; j<Mat.p[sl2l-Mat.lb2+i+1] ; j++) {
1978  if(Mat.ind[j] >= sl1l-Mat.lb1 && Mat.ind[j] <= sl1u-Mat.lb1) {
1979  A.ind.push_back(Mat.ind[j]-(sl1l-Mat.lb1));
1980  A.x.push_back(Mat.x[j]);
1981  A.p[i+1]++;
1982  }
1983  }
1984  }
1985 
1986  //Zeiger auf A fuer Datenmanipulationen
1987  M = &Mat;
1988  }
1989 
1990  scmatrix_slice(const scmatrix& Mat, int sl1l, int sl1u, int sl2l, int sl2u) {
1991  A.lb1 = sl1l;
1992  A.lb2 = sl2l;
1993  A.ub1 = sl1u;
1994  A.ub2 = sl2u;
1995  A.m = sl1u-sl1l+1;
1996  A.n = sl2u-sl2l+1;
1997 
1998  //Kopieren der Werte aus A
1999  A.p = std::vector<int>(A.n+1, 0);
2000  A.ind.reserve(A.m + A.n);
2001  A.x.reserve(A.m + A.n);
2002 
2003  for(int i=0 ; i<A.n ; i++) {
2004  A.p[i+1] = A.p[i];
2005  for(int j=Mat.p[sl2l-Mat.lb2+i] ; j<Mat.p[sl2l-Mat.lb2+i+1] ; j++) {
2006  if(Mat.ind[j] >= sl1l-Mat.lb1 && Mat.ind[j] <= sl1u-Mat.lb1) {
2007  A.ind.push_back(Mat.ind[j]-(sl1l-Mat.lb1));
2008  A.x.push_back(Mat.x[j]);
2009  A.p[i+1]++;
2010  }
2011  }
2012  }
2013 
2014  //Zeiger auf A fuer Datenmanipulationen
2015  M = const_cast<scmatrix*>(&Mat);
2016  }
2017 
2018  public:
2021  return sl_ms_assign<scmatrix_slice, real, std::vector<complex>::iterator, complex>(*this,C);
2022  }
2023 
2026  return sl_ms_assign<scmatrix_slice, complex, std::vector<complex>::iterator, complex>(*this,C);
2027  }
2028 
2031  return slsp_mm_assign<scmatrix_slice, srmatrix, std::vector<complex>::iterator>(*this,C);
2032  }
2033 
2036  return slsp_mm_assign<scmatrix_slice, scmatrix, std::vector<complex>::iterator>(*this,C);
2037  }
2038 
2041  return slf_mm_assign<scmatrix_slice, rmatrix, std::vector<complex>::iterator, complex>(*this,C);
2042  }
2043 
2046  return slf_mm_assign<scmatrix_slice, cmatrix, std::vector<complex>::iterator, complex>(*this,C);
2047  }
2048 
2051  return slf_mm_assign<scmatrix_slice, rmatrix_slice, std::vector<complex>::iterator, complex>(*this,C);
2052  }
2053 
2056  return slf_mm_assign<scmatrix_slice, cmatrix_slice, std::vector<complex>::iterator, complex>(*this,C);
2057  }
2058 
2061  *this = C.A;
2062  return *this;
2063  }
2064 
2067  *this = C.A;
2068  return *this;
2069  }
2070 
2073  *this = A*M.A;
2074  return *this;
2075  }
2076 
2079  *this = A*M.A;
2080  return *this;
2081  }
2082 
2085  *this = A*M;
2086  return *this;
2087  }
2088 
2091  *this = A*M;
2092  return *this;
2093  }
2094 
2097  *this = A*M;
2098  return *this;
2099  }
2100 
2103  *this = A*M;
2104  return *this;
2105  }
2106 
2109  *this = A*M;
2110  return *this;
2111  }
2112 
2115  *this = A*M;
2116  return *this;
2117  }
2118 
2121  *this = A*r;
2122  return *this;
2123  }
2124 
2127  *this = A*r;
2128  return *this;
2129  }
2130 
2133  *this = A/r;
2134  return *this;
2135  }
2136 
2139  *this = A/r;
2140  return *this;
2141  }
2142 
2145  *this = A+M.A;
2146  return *this;
2147  }
2148 
2151  *this = A+M.A;
2152  return *this;
2153  }
2154 
2157  *this = A+M;
2158  return *this;
2159  }
2160 
2163  *this = A+M;
2164  return *this;
2165  }
2166 
2169  *this = A+M;
2170  return *this;
2171  }
2172 
2175  *this = A+M;
2176  return *this;
2177  }
2178 
2181  *this = A+M;
2182  return *this;
2183  }
2184 
2187  *this = A+M;
2188  return *this;
2189  }
2190 
2193  *this = A-M.A;
2194  return *this;
2195  }
2196 
2199  *this = A-M.A;
2200  return *this;
2201  }
2202 
2205  *this = A-M;
2206  return *this;
2207  }
2208 
2211  *this = A-M;
2212  return *this;
2213  }
2214 
2217  *this = A-M;
2218  return *this;
2219  }
2220 
2223  *this = A-M;
2224  return *this;
2225  }
2226 
2229  *this = A-M;
2230  return *this;
2231  }
2232 
2235  *this = A-M;
2236  return *this;
2237  }
2238 
2240 
2244  const complex operator()(const int i, const int j) const {
2245 #if(CXSC_INDEX_CHECK)
2246  if(i<A.lb1 || i>A.ub1 || j<A.lb2 || j>A.ub2)
2247  cxscthrow(ELEMENT_NOT_IN_VEC("scmatrix_slice::operator()(int, int)"));
2248 #endif
2249  complex r = A(i,j);
2250  return r;
2251  }
2252 
2254 
2258  complex& element(const int i, const int j) {
2259 #if(CXSC_INDEX_CHECK)
2260  if(i<A.lb1 || i>A.ub1 || j<A.lb2 || j>A.ub2)
2261  cxscthrow(ELEMENT_NOT_IN_VEC("scmatrix_slice::element(int, int)"));
2262 #endif
2263  return M->element(i,j);
2264  }
2265 
2267  scmatrix_subv operator[](const int);
2269  scmatrix_subv operator[](const cxscmatrix_column&);
2271  const scmatrix_subv operator[](const int) const;
2273  const scmatrix_subv operator[](const cxscmatrix_column&) const;
2274 
2275  friend std::ostream& operator<<(std::ostream&, const scmatrix_slice&);
2276  friend std::istream& operator>>(std::istream&, const scmatrix_subv&);
2277 
2278  friend int Lb(const scmatrix_slice&, const int);
2279  friend int Ub(const scmatrix_slice&, const int);
2280  friend srmatrix Re(const scmatrix_slice&);
2281  friend srmatrix Im(const scmatrix_slice&);
2282  friend int RowLen(const scmatrix_slice&);
2283  friend int ColLen(const scmatrix_slice&);
2284 
2285  friend class srmatrix;
2286  friend class srmatrix_subv;
2287  friend class srvector;
2288  friend class scmatrix;
2289  friend class scmatrix_subv;
2290  friend class scvector;
2291  friend class scimatrix;
2292  friend class scimatrix_subv;
2293  friend class scimatrix_slice;
2294  friend class scivector;
2295  friend class cmatrix;
2296  friend class cimatrix;
2297 
2298 #include "matrix_friend_declarations.inl"
2299 };
2300 
2302  dat = new complex[A.A.m*A.A.n];
2303  lb1 = A.A.lb1; lb2 = A.A.lb2; ub1 = A.A.ub1; ub2 = A.A.ub2;
2304  xsize = A.A.n;
2305  ysize = A.A.m;
2306  *this = 0.0;
2307  for(int j=0 ; j<A.A.n ; j++) {
2308  for(int k=A.A.p[j] ; k<A.A.p[j+1] ; k++) {
2309  dat[A.A.ind[k]*A.A.n+j] = A.A.x[k];
2310  }
2311  }
2312 }
2313 
2315  dat = new complex[A.A.m*A.A.n];
2316  lb1 = A.A.lb1; lb2 = A.A.lb2; ub1 = A.A.ub1; ub2 = A.A.ub2;
2317  xsize = A.A.n;
2318  ysize = A.A.m;
2319  *this = 0.0;
2320  for(int j=0 ; j<A.A.n ; j++) {
2321  for(int k=A.A.p[j] ; k<A.A.p[j+1] ; k++) {
2322  dat[A.A.ind[k]*A.A.n+j] = A.A.x[k];
2323  }
2324  }
2325 }
2326 
2328 inline int RowLen(const scmatrix_slice& S) {
2329  return RowLen(S.A);
2330 }
2331 
2333 inline int ColLen(const scmatrix_slice& S) {
2334  return ColLen(S.A);
2335 }
2336 
2337 inline scmatrix_slice scmatrix::operator()(const int i, const int j, const int k, const int l) {
2338 #if(CXSC_INDEX_CHECK)
2339  if(i<lb1 || j>ub1 || k<lb2 || l>ub2)
2340  cxscthrow(ROW_OR_COL_NOT_IN_MAT("scmatrix::operator()(int, int, int, int)"));
2341 #endif
2342  return scmatrix_slice(*this, i, j, k, l);
2343 }
2344 
2345 inline const scmatrix_slice scmatrix::operator()(const int i, const int j, const int k, const int l) const {
2346 #if(CXSC_INDEX_CHECK)
2347  if(i<lb1 || j>ub1 || k<lb2 || l>ub2)
2348  cxscthrow(ROW_OR_COL_NOT_IN_MAT("scmatrix::operator()(int, int, int, int) const"));
2349 #endif
2350  return scmatrix_slice(*this, i, j, k, l);
2351 }
2352 
2354  *this = S.A;
2355  return *this;
2356 }
2357 
2359  *this = S.A;
2360  return *this;
2361 }
2362 
2364 inline int Lb(const scmatrix_slice& S, const int i) {
2365  return Lb(S.A, i);
2366 }
2367 
2369 inline int Ub(const scmatrix_slice& S, const int i) {
2370  return Ub(S.A, i);
2371 }
2372 
2374 inline srmatrix Re(const scmatrix_slice& S) {
2375  return Re(S.A);
2376 }
2377 
2379 inline srmatrix Im(const scmatrix_slice& S) {
2380  return Im(S.A);
2381 }
2382 
2384  m = S.A.m;
2385  n = S.A.n;
2386  lb1 = S.A.lb1;
2387  ub1 = S.A.ub1;
2388  lb2 = S.A.lb2;
2389  ub2 = S.A.ub2;
2390  *this = S.A;
2391 }
2392 
2394  m = S.A.m;
2395  n = S.A.n;
2396  lb1 = S.A.lb1;
2397  ub1 = S.A.ub1;
2398  lb2 = S.A.lb2;
2399  ub2 = S.A.ub2;
2400  *this = S.A;
2401 }
2402 
2404 inline scmatrix operator-(const scmatrix_slice& M) {
2405  return sp_m_negative<scmatrix,scmatrix>(M.A);
2406 }
2407 
2409 inline scmatrix operator+(const scmatrix_slice& M) {
2410  return M.A;
2411 }
2412 
2414 
2420 inline scmatrix operator*(const scmatrix_slice& M1, const srmatrix_slice& M2) {
2421  return spsp_mm_mult<scmatrix,srmatrix,scmatrix,sparse_cdot,complex>(M1.A,M2.A);
2422 }
2423 
2425 
2431 inline scmatrix operator*(const srmatrix_slice& M1, const scmatrix_slice& M2) {
2432  return spsp_mm_mult<srmatrix,scmatrix,scmatrix,sparse_cdot,complex>(M1.A,M2.A);
2433 }
2434 
2436 
2442 inline scmatrix operator*(const scmatrix_slice& M1, const scmatrix_slice& M2) {
2443  return spsp_mm_mult<scmatrix,scmatrix,scmatrix,sparse_cdot,complex>(M1.A,M2.A);
2444 }
2445 
2447 
2453 inline scmatrix operator*(const scmatrix_slice& M1, const srmatrix& M2) {
2454  return spsp_mm_mult<scmatrix,srmatrix,scmatrix,sparse_cdot,complex>(M1.A,M2);
2455 }
2456 
2458 
2464 inline scmatrix operator*(const srmatrix_slice& M1, const scmatrix& M2) {
2465  return spsp_mm_mult<srmatrix,scmatrix,scmatrix,sparse_cdot,complex>(M1.A,M2);
2466 }
2467 
2469 
2475 inline scmatrix operator*(const scmatrix_slice& M1, const scmatrix& M2) {
2476  return spsp_mm_mult<scmatrix,scmatrix,scmatrix,sparse_cdot,complex>(M1.A,M2);
2477 }
2478 
2480 
2486 inline scmatrix operator*(const scmatrix& M1, const srmatrix_slice& M2) {
2487  return spsp_mm_mult<scmatrix,srmatrix,scmatrix,sparse_cdot,complex>(M1,M2.A);
2488 }
2489 
2491 
2497 inline scmatrix operator*(const srmatrix& M1, const scmatrix_slice& M2) {
2498  return spsp_mm_mult<srmatrix,scmatrix,scmatrix,sparse_cdot,complex>(M1,M2.A);
2499 }
2500 
2502 
2508 inline scmatrix operator*(const scmatrix& M1, const scmatrix_slice& M2) {
2509  return spsp_mm_mult<scmatrix,scmatrix,scmatrix,sparse_cdot,complex>(M1,M2.A);
2510 }
2511 
2513 
2519 inline cmatrix operator*(const scmatrix_slice& M1, const rmatrix& M2) {
2520  return spf_mm_mult<scmatrix,rmatrix,cmatrix,sparse_cdot>(M1.A,M2);
2521 }
2522 
2524 
2530 inline cmatrix operator*(const srmatrix_slice& M1, const cmatrix& M2) {
2531  return spf_mm_mult<srmatrix,cmatrix,cmatrix,sparse_cdot>(M1.A,M2);
2532 }
2533 
2535 
2541 inline cmatrix operator*(const scmatrix_slice& M1, const cmatrix& M2) {
2542  return spf_mm_mult<scmatrix,cmatrix,cmatrix,sparse_cdot>(M1.A,M2);
2543 }
2544 
2546 
2552 inline cmatrix operator*(const cmatrix& M1, const srmatrix_slice& M2) {
2553  return fsp_mm_mult<cmatrix,srmatrix,cmatrix,sparse_cdot>(M1,M2.A);
2554 }
2555 
2557 
2563 inline cmatrix operator*(const rmatrix& M1, const scmatrix_slice& M2) {
2564  return fsp_mm_mult<rmatrix,scmatrix,cmatrix,sparse_cdot>(M1,M2.A);
2565 }
2566 
2568 
2574 inline cmatrix operator*(const cmatrix& M1, const scmatrix_slice& M2) {
2575  return fsp_mm_mult<cmatrix,scmatrix,cmatrix,sparse_cdot>(M1,M2.A);
2576 }
2577 
2579 
2585 inline cmatrix operator*(const scmatrix_slice& M1, const rmatrix_slice& M2) {
2586  return spf_mm_mult<scmatrix,rmatrix_slice,cmatrix,sparse_cdot>(M1.A,M2);
2587 }
2588 
2590 
2596 inline cmatrix operator*(const srmatrix_slice& M1, const cmatrix_slice& M2) {
2597  return spf_mm_mult<srmatrix,cmatrix_slice,cmatrix,sparse_cdot>(M1.A,M2);
2598 }
2599 
2601 
2607 inline cmatrix operator*(const scmatrix_slice& M1, const cmatrix_slice& M2) {
2608  return spf_mm_mult<scmatrix,cmatrix_slice,cmatrix,sparse_cdot>(M1.A,M2);
2609 }
2610 
2612 
2618 inline cmatrix operator*(const cmatrix_slice& M1, const srmatrix_slice& M2) {
2619  return fsp_mm_mult<cmatrix,srmatrix,cmatrix,sparse_cdot>(M1,M2.A);
2620 }
2621 
2623 
2629 inline cmatrix operator*(const rmatrix_slice& M1, const scmatrix_slice& M2) {
2630  return fsp_mm_mult<rmatrix,scmatrix,cmatrix,sparse_cdot>(M1,M2.A);
2631 }
2632 
2634 
2640 inline cmatrix operator*(const cmatrix_slice& M1, const scmatrix_slice& M2) {
2641  return fsp_mm_mult<cmatrix,scmatrix,cmatrix,sparse_cdot>(M1,M2.A);
2642 }
2643 
2645 
2651 inline scvector operator*(const scmatrix_slice& M, const srvector& v) {
2652  return spsp_mv_mult<scmatrix,srvector,scvector,sparse_cdot,complex>(M.A,v);
2653 }
2654 
2656 
2662 inline scvector operator*(const srmatrix_slice& M, const scvector& v) {
2663  return spsp_mv_mult<srmatrix,scvector,scvector,sparse_cdot,complex>(M.A,v);
2664 }
2665 
2667 
2673 inline scvector operator*(const scmatrix_slice& M, const scvector& v) {
2674  return spsp_mv_mult<scmatrix,scvector,scvector,sparse_cdot,complex>(M.A,v);
2675 }
2676 
2678 
2684 inline scvector operator*(const scmatrix_slice& M, const srvector_slice& v) {
2685  return spsl_mv_mult<scmatrix,srvector_slice,scvector,sparse_cdot,complex>(M.A,v);
2686 }
2687 
2689 
2695 inline scvector operator*(const srmatrix_slice& M, const scvector_slice& v) {
2696  return spsl_mv_mult<srmatrix,scvector_slice,scvector,sparse_cdot,complex>(M.A,v);
2697 }
2698 
2700 
2706 inline scvector operator*(const scmatrix_slice& M, const scvector_slice& v) {
2707  return spsl_mv_mult<scmatrix,scvector_slice,scvector,sparse_cdot,complex>(M.A,v);
2708 }
2709 
2711 
2717 inline cvector operator*(const scmatrix_slice& M, const rvector& v) {
2718  return spf_mv_mult<scmatrix,rvector,cvector,sparse_cdot>(M.A,v);
2719 }
2720 
2722 
2728 inline cvector operator*(const srmatrix_slice& M, const cvector& v) {
2729  return spf_mv_mult<srmatrix,cvector,cvector,sparse_cdot>(M.A,v);
2730 }
2731 
2733 
2739 inline cvector operator*(const scmatrix_slice& M, const cvector& v) {
2740  return spf_mv_mult<scmatrix,cvector,cvector,sparse_cdot>(M.A,v);
2741 }
2742 
2744 
2750 inline cvector operator*(const scmatrix_slice& M, const rvector_slice& v) {
2751  return spf_mv_mult<scmatrix,rvector_slice,cvector,sparse_cdot>(M.A,v);
2752 }
2753 
2755 
2761 inline cvector operator*(const srmatrix_slice& M, const cvector_slice& v) {
2762  return spf_mv_mult<srmatrix,cvector_slice,cvector,sparse_cdot>(M.A,v);
2763 }
2764 
2766 
2772 inline cvector operator*(const scmatrix_slice& M, const cvector_slice& v) {
2773  return spf_mv_mult<scmatrix,cvector_slice,cvector,sparse_cdot>(M.A,v);
2774 }
2775 
2777 inline scmatrix operator*(const scmatrix_slice& M, const real& r) {
2778  return sp_ms_mult<scmatrix,real,scmatrix>(M.A,r);
2779 }
2780 
2782 inline scmatrix operator*(const scmatrix_slice& M, const complex& r) {
2783  return sp_ms_mult<scmatrix,complex,scmatrix>(M.A,r);
2784 }
2785 
2787 inline scmatrix operator*(const srmatrix_slice& M, const complex& r) {
2788  return sp_ms_mult<srmatrix,complex,scmatrix>(M.A,r);
2789 }
2790 
2792 inline scmatrix operator/(const scmatrix_slice& M, const real& r) {
2793  return sp_ms_div<scmatrix,real,scmatrix>(M.A,r);
2794 }
2795 
2797 inline scmatrix operator/(const scmatrix_slice& M, const complex& r) {
2798  return sp_ms_div<scmatrix,complex,scmatrix>(M.A,r);
2799 }
2800 
2802 inline scmatrix operator/(const srmatrix_slice& M, const complex& r) {
2803  return sp_ms_div<srmatrix,complex,scmatrix>(M.A,r);
2804 }
2805 
2807 inline scmatrix operator*(const real& r, const scmatrix_slice& M) {
2808  return sp_sm_mult<real,scmatrix,scmatrix>(r,M.A);
2809 }
2810 
2812 inline scmatrix operator*(const complex& r, const srmatrix_slice& M) {
2813  return sp_sm_mult<complex,srmatrix,scmatrix>(r,M.A);
2814 }
2815 
2817 inline scmatrix operator*(const complex& r, const scmatrix_slice& M) {
2818  return sp_sm_mult<complex,scmatrix,scmatrix>(r,M.A);
2819 }
2820 
2822 inline scmatrix operator+(const scmatrix_slice& M1, const srmatrix_slice& M2) {
2823  return spsp_mm_add<scmatrix,srmatrix,scmatrix,complex>(M1.A,M2.A);
2824 }
2825 
2827 inline scmatrix operator+(const srmatrix_slice& M1, const scmatrix_slice& M2) {
2828  return spsp_mm_add<srmatrix,scmatrix,scmatrix,complex>(M1.A,M2.A);
2829 }
2830 
2832 inline scmatrix operator+(const scmatrix_slice& M1, const scmatrix_slice& M2) {
2833  return spsp_mm_add<scmatrix,scmatrix,scmatrix,complex>(M1.A,M2.A);
2834 }
2835 
2837 inline scmatrix operator+(const scmatrix_slice& M1, const srmatrix& M2) {
2838  return spsp_mm_add<scmatrix,srmatrix,scmatrix,complex>(M1.A,M2);
2839 }
2840 
2842 inline scmatrix operator+(const srmatrix_slice& M1, const scmatrix& M2) {
2843  return spsp_mm_add<srmatrix,scmatrix,scmatrix,complex>(M1.A,M2);
2844 }
2845 
2847 inline scmatrix operator+(const scmatrix_slice& M1, const scmatrix& M2) {
2848  return spsp_mm_add<scmatrix,scmatrix,scmatrix,complex>(M1.A,M2);
2849 }
2850 
2852 inline scmatrix operator+(const scmatrix& M1, const srmatrix_slice& M2) {
2853  return spsp_mm_add<scmatrix,srmatrix,scmatrix,complex>(M1,M2.A);
2854 }
2855 
2857 inline scmatrix operator+(const srmatrix& M1, const scmatrix_slice& M2) {
2858  return spsp_mm_add<srmatrix,scmatrix,scmatrix,complex>(M1,M2.A);
2859 }
2860 
2862 inline scmatrix operator+(const scmatrix& M1, const scmatrix_slice& M2) {
2863  return spsp_mm_add<scmatrix,scmatrix,scmatrix,complex>(M1,M2.A);
2864 }
2865 
2867 inline cmatrix operator+(const scmatrix_slice& M1, const rmatrix& M2) {
2868  return spf_mm_add<scmatrix,rmatrix,cmatrix>(M1.A,M2);
2869 }
2870 
2872 inline cmatrix operator+(const srmatrix_slice& M1, const cmatrix& M2) {
2873  return spf_mm_add<srmatrix,cmatrix,cmatrix>(M1.A,M2);
2874 }
2875 
2877 inline cmatrix operator+(const scmatrix_slice& M1, const cmatrix& M2) {
2878  return spf_mm_add<scmatrix,cmatrix,cmatrix>(M1.A,M2);
2879 }
2880 
2882 inline cmatrix operator+(const cmatrix& M1, const srmatrix_slice& M2) {
2883  return fsp_mm_add<cmatrix,srmatrix,cmatrix>(M1,M2.A);
2884 }
2885 
2887 inline cmatrix operator+(const rmatrix& M1, const scmatrix_slice& M2) {
2888  return fsp_mm_add<rmatrix,scmatrix,cmatrix>(M1,M2.A);
2889 }
2890 
2892 inline cmatrix operator+(const cmatrix& M1, const scmatrix_slice& M2) {
2893  return fsp_mm_add<cmatrix,scmatrix,cmatrix>(M1,M2.A);
2894 }
2895 
2897 inline cmatrix operator+(const scmatrix_slice& M1, const rmatrix_slice& M2) {
2898  return spf_mm_add<scmatrix,rmatrix_slice,cmatrix>(M1.A,M2);
2899 }
2900 
2902 inline cmatrix operator+(const srmatrix_slice& M1, const cmatrix_slice& M2) {
2903  return spf_mm_add<srmatrix,cmatrix_slice,cmatrix>(M1.A,M2);
2904 }
2905 
2907 inline cmatrix operator+(const scmatrix_slice& M1, const cmatrix_slice& M2) {
2908  return spf_mm_add<scmatrix,cmatrix_slice,cmatrix>(M1.A,M2);
2909 }
2910 
2912 inline cmatrix operator+(const cmatrix_slice& M1, const srmatrix_slice& M2) {
2913  return fsp_mm_add<cmatrix_slice,srmatrix,cmatrix>(M1,M2.A);
2914 }
2915 
2917 inline cmatrix operator+(const rmatrix_slice& M1, const scmatrix_slice& M2) {
2918  return fsp_mm_add<rmatrix_slice,scmatrix,cmatrix>(M1,M2.A);
2919 }
2920 
2922 inline cmatrix operator+(const cmatrix_slice& M1, const scmatrix_slice& M2) {
2923  return fsp_mm_add<cmatrix_slice,scmatrix,cmatrix>(M1,M2.A);
2924 }
2925 
2927 inline scmatrix operator-(const scmatrix_slice& M1, const srmatrix_slice& M2) {
2928  return spsp_mm_sub<scmatrix,srmatrix,scmatrix,complex>(M1.A,M2.A);
2929 }
2930 
2932 inline scmatrix operator-(const srmatrix_slice& M1, const scmatrix_slice& M2) {
2933  return spsp_mm_sub<srmatrix,scmatrix,scmatrix,complex>(M1.A,M2.A);
2934 }
2935 
2937 inline scmatrix operator-(const scmatrix_slice& M1, const scmatrix_slice& M2) {
2938  return spsp_mm_sub<scmatrix,scmatrix,scmatrix,complex>(M1.A,M2.A);
2939 }
2940 
2942 inline scmatrix operator-(const scmatrix_slice& M1, const srmatrix& M2) {
2943  return spsp_mm_sub<scmatrix,srmatrix,scmatrix,complex>(M1.A,M2);
2944 }
2945 
2947 inline scmatrix operator-(const srmatrix_slice& M1, const scmatrix& M2) {
2948  return spsp_mm_sub<srmatrix,scmatrix,scmatrix,complex>(M1.A,M2);
2949 }
2950 
2952 inline scmatrix operator-(const scmatrix_slice& M1, const scmatrix& M2) {
2953  return spsp_mm_sub<scmatrix,scmatrix,scmatrix,complex>(M1.A,M2);
2954 }
2955 
2957 inline scmatrix operator-(const scmatrix& M1, const srmatrix_slice& M2) {
2958  return spsp_mm_sub<scmatrix,srmatrix,scmatrix,complex>(M1,M2.A);
2959 }
2960 
2962 inline scmatrix operator-(const srmatrix& M1, const scmatrix_slice& M2) {
2963  return spsp_mm_sub<srmatrix,scmatrix,scmatrix,complex>(M1,M2.A);
2964 }
2965 
2967 inline scmatrix operator-(const scmatrix& M1, const scmatrix_slice& M2) {
2968  return spsp_mm_sub<scmatrix,scmatrix,scmatrix,complex>(M1,M2.A);
2969 }
2970 
2972 inline cmatrix operator-(const scmatrix_slice& M1, const rmatrix& M2) {
2973  return spf_mm_sub<scmatrix,rmatrix,cmatrix>(M1.A,M2);
2974 }
2975 
2977 inline cmatrix operator-(const srmatrix_slice& M1, const cmatrix& M2) {
2978  return spf_mm_sub<srmatrix,cmatrix,cmatrix>(M1.A,M2);
2979 }
2980 
2982 inline cmatrix operator-(const scmatrix_slice& M1, const cmatrix& M2) {
2983  return spf_mm_sub<scmatrix,cmatrix,cmatrix>(M1.A,M2);
2984 }
2985 
2987 inline cmatrix operator-(const cmatrix& M1, const srmatrix_slice& M2) {
2988  return fsp_mm_sub<cmatrix,srmatrix,cmatrix>(M1,M2.A);
2989 }
2990 
2992 inline cmatrix operator-(const rmatrix& M1, const scmatrix_slice& M2) {
2993  return fsp_mm_sub<rmatrix,scmatrix,cmatrix>(M1,M2.A);
2994 }
2995 
2997 inline cmatrix operator-(const cmatrix& M1, const scmatrix_slice& M2) {
2998  return fsp_mm_sub<cmatrix,scmatrix,cmatrix>(M1,M2.A);
2999 }
3000 
3002 inline cmatrix operator-(const scmatrix_slice& M1, const rmatrix_slice& M2) {
3003  return spf_mm_sub<scmatrix,rmatrix_slice,cmatrix>(M1.A,M2);
3004 }
3005 
3007 inline cmatrix operator-(const srmatrix_slice& M1, const cmatrix_slice& M2) {
3008  return spf_mm_sub<srmatrix,cmatrix_slice,cmatrix>(M1.A,M2);
3009 }
3010 
3012 inline cmatrix operator-(const scmatrix_slice& M1, const cmatrix_slice& M2) {
3013  return spf_mm_sub<scmatrix,cmatrix_slice,cmatrix>(M1.A,M2);
3014 }
3015 
3017 inline cmatrix operator-(const cmatrix_slice& M1, const srmatrix_slice& M2) {
3018  return fsp_mm_sub<cmatrix_slice,srmatrix,cmatrix>(M1,M2.A);
3019 }
3020 
3022 inline cmatrix operator-(const rmatrix_slice& M1, const scmatrix_slice& M2) {
3023  return fsp_mm_sub<rmatrix_slice,scmatrix,cmatrix>(M1,M2.A);
3024 }
3025 
3027 inline cmatrix operator-(const cmatrix_slice& M1, const scmatrix_slice& M2) {
3028  return fsp_mm_sub<cmatrix_slice,scmatrix,cmatrix>(M1,M2.A);
3029 }
3030 
3032  *this = rmatrix(M);
3033  return *this;
3034 }
3035 
3037  *this = cmatrix(M);
3038  return *this;
3039 }
3040 
3042  *this += M.A;
3043  return *this;
3044 }
3045 
3047  *this += M.A;
3048  return *this;
3049 }
3050 
3052  *this += M.A;
3053  return *this;
3054 }
3055 
3057  *this += M.A;
3058  return *this;
3059 }
3060 
3061 inline cmatrix& cmatrix::operator-=(const srmatrix_slice& M) {
3062  *this -= M.A;
3063  return *this;
3064 }
3065 
3066 inline cmatrix& cmatrix::operator-=(const scmatrix_slice& M) {
3067  *this -= M.A;
3068  return *this;
3069 }
3070 
3071 inline cmatrix_slice& cmatrix_slice::operator-=(const srmatrix_slice& M) {
3072  *this -= M.A;
3073  return *this;
3074 }
3075 
3076 inline cmatrix_slice& cmatrix_slice::operator-=(const scmatrix_slice& M) {
3077  *this -= M.A;
3078  return *this;
3079 }
3080 
3082  *this *= M.A;
3083  return *this;
3084 }
3085 
3087  *this *= M.A;
3088  return *this;
3089 }
3090 
3092  *this *= M.A;
3093  return *this;
3094 }
3095 
3097  *this *= M.A;
3098  return *this;
3099 }
3100 
3102 inline bool operator==(const scmatrix_slice& M1, const srmatrix_slice& M2) {
3103  return spsp_mm_comp(M1.A,M2.A);
3104 }
3105 
3107 inline bool operator==(const srmatrix_slice& M1, const scmatrix_slice& M2) {
3108  return spsp_mm_comp(M1.A,M2.A);
3109 }
3110 
3112 inline bool operator==(const scmatrix_slice& M1, const scmatrix_slice& M2) {
3113  return spsp_mm_comp(M1.A,M2.A);
3114 }
3115 
3117 inline bool operator==(const scmatrix_slice& M1, const srmatrix& M2) {
3118  return spsp_mm_comp(M1.A,M2);
3119 }
3120 
3122 inline bool operator==(const srmatrix_slice& M1, const scmatrix& M2) {
3123  return spsp_mm_comp(M1.A,M2);
3124 }
3125 
3127 inline bool operator==(const scmatrix_slice& M1, const scmatrix& M2) {
3128  return spsp_mm_comp(M1.A,M2);
3129 }
3130 
3132 inline bool operator==(const scmatrix& M1, const srmatrix_slice& M2) {
3133  return spsp_mm_comp(M1,M2.A);
3134 }
3135 
3137 inline bool operator==(const srmatrix& M1, const scmatrix_slice& M2) {
3138  return spsp_mm_comp(M1,M2.A);
3139 }
3140 
3142 inline bool operator==(const scmatrix& M1, const scmatrix_slice& M2) {
3143  return spsp_mm_comp(M1,M2.A);
3144 }
3145 
3147 inline bool operator==(const scmatrix_slice& M1, const rmatrix& M2) {
3148  return spf_mm_comp(M1.A,M2);
3149 }
3150 
3152 inline bool operator==(const srmatrix_slice& M1, const cmatrix& M2) {
3153  return spf_mm_comp(M1.A,M2);
3154 }
3155 
3157 inline bool operator==(const scmatrix_slice& M1, const cmatrix& M2) {
3158  return spf_mm_comp(M1.A,M2);
3159 }
3160 
3162 inline bool operator==(const cmatrix& M1, const srmatrix_slice& M2) {
3163  return fsp_mm_comp(M1,M2.A);
3164 }
3165 
3167 inline bool operator==(const rmatrix& M1, const scmatrix_slice& M2) {
3168  return fsp_mm_comp(M1,M2.A);
3169 }
3170 
3172 inline bool operator==(const cmatrix& M1, const scmatrix_slice& M2) {
3173  return fsp_mm_comp(M1,M2.A);
3174 }
3175 
3177 inline bool operator==(const cmatrix_slice& M1, const srmatrix_slice& M2) {
3178  return fsp_mm_comp(M1,M2.A);
3179 }
3180 
3182 inline bool operator==(const rmatrix_slice& M1, const scmatrix_slice& M2) {
3183  return fsp_mm_comp(M1,M2.A);
3184 }
3185 
3187 inline bool operator==(const cmatrix_slice& M1, const scmatrix_slice& M2) {
3188  return fsp_mm_comp(M1,M2.A);
3189 }
3190 
3192 inline bool operator==(const scmatrix_slice& M1, const rmatrix_slice& M2) {
3193  return spf_mm_comp(M1.A,M2);
3194 }
3195 
3197 inline bool operator==(const srmatrix_slice& M1, const cmatrix_slice& M2) {
3198  return spf_mm_comp(M1.A,M2);
3199 }
3200 
3202 inline bool operator==(const scmatrix_slice& M1, const cmatrix_slice& M2) {
3203  return spf_mm_comp(M1.A,M2);
3204 }
3205 
3207 inline bool operator!=(const scmatrix_slice& M1, const srmatrix_slice& M2) {
3208  return !spsp_mm_comp(M1.A,M2.A);
3209 }
3210 
3212 inline bool operator!=(const srmatrix_slice& M1, const scmatrix_slice& M2) {
3213  return !spsp_mm_comp(M1.A,M2.A);
3214 }
3215 
3217 inline bool operator!=(const scmatrix_slice& M1, const scmatrix_slice& M2) {
3218  return !spsp_mm_comp(M1.A,M2.A);
3219 }
3220 
3222 inline bool operator!=(const scmatrix_slice& M1, const srmatrix& M2) {
3223  return !spsp_mm_comp(M1.A,M2);
3224 }
3225 
3227 inline bool operator!=(const srmatrix_slice& M1, const scmatrix& M2) {
3228  return !spsp_mm_comp(M1.A,M2);
3229 }
3230 
3232 inline bool operator!=(const scmatrix_slice& M1, const scmatrix& M2) {
3233  return !spsp_mm_comp(M1.A,M2);
3234 }
3235 
3237 inline bool operator!=(const scmatrix& M1, const srmatrix_slice& M2) {
3238  return !spsp_mm_comp(M1,M2.A);
3239 }
3240 
3242 inline bool operator!=(const srmatrix& M1, const scmatrix_slice& M2) {
3243  return !spsp_mm_comp(M1,M2.A);
3244 }
3245 
3247 inline bool operator!=(const scmatrix& M1, const scmatrix_slice& M2) {
3248  return !spsp_mm_comp(M1,M2.A);
3249 }
3250 
3252 inline bool operator!=(const scmatrix_slice& M1, const rmatrix& M2) {
3253  return !spf_mm_comp(M1.A,M2);
3254 }
3255 
3257 inline bool operator!=(const srmatrix_slice& M1, const cmatrix& M2) {
3258  return !spf_mm_comp(M1.A,M2);
3259 }
3260 
3262 inline bool operator!=(const scmatrix_slice& M1, const cmatrix& M2) {
3263  return !spf_mm_comp(M1.A,M2);
3264 }
3265 
3267 inline bool operator!=(const cmatrix& M1, const srmatrix_slice& M2) {
3268  return !fsp_mm_comp(M1,M2.A);
3269 }
3270 
3272 inline bool operator!=(const rmatrix& M1, const scmatrix_slice& M2) {
3273  return !fsp_mm_comp(M1,M2.A);
3274 }
3275 
3277 inline bool operator!=(const cmatrix& M1, const scmatrix_slice& M2) {
3278  return !fsp_mm_comp(M1,M2.A);
3279 }
3280 
3282 inline bool operator!=(const cmatrix_slice& M1, const srmatrix_slice& M2) {
3283  return !fsp_mm_comp(M1,M2.A);
3284 }
3285 
3287 inline bool operator!=(const rmatrix_slice& M1, const scmatrix_slice& M2) {
3288  return !fsp_mm_comp(M1,M2.A);
3289 }
3290 
3292 inline bool operator!=(const cmatrix_slice& M1, const scmatrix_slice& M2) {
3293  return !fsp_mm_comp(M1,M2.A);
3294 }
3295 
3297 inline bool operator!=(const scmatrix_slice& M1, const rmatrix_slice& M2) {
3298  return !spf_mm_comp(M1.A,M2);
3299 }
3300 
3302 inline bool operator!=(const srmatrix_slice& M1, const cmatrix_slice& M2) {
3303  return !spf_mm_comp(M1.A,M2);
3304 }
3305 
3307 inline bool operator!=(const scmatrix_slice& M1, const cmatrix_slice& M2) {
3308  return !spf_mm_comp(M1.A,M2);
3309 }
3310 
3312 inline bool operator!(const scmatrix_slice& M) {
3313  return sp_m_not(M.A);
3314 }
3315 
3317 
3322 inline std::ostream& operator<<(std::ostream& os, const scmatrix_slice& M) {
3323  return sp_m_output<scmatrix,complex>(os, M.A);
3324 }
3325 
3327 
3332 inline std::istream& operator>>(std::istream& is, scmatrix_slice& M) {
3333  scmatrix tmp(M.A.m, M.A.n);
3334  is >> tmp;
3335  M = tmp;
3336  return is;
3337 }
3338 
3340 
3346  private:
3347  scmatrix_slice dat;
3348  bool row;
3349  int index;
3350 
3351  scmatrix_subv(scmatrix& A, bool r, int i, int j, int k, int l) : dat(A,i,j,k,l), row(r) {
3352  if(row) index=i; else index=k;
3353  }
3354 
3355  scmatrix_subv(const scmatrix& A, bool r, int i, int j, int k, int l) : dat(A,i,j,k,l), row(r) {
3356  if(row) index=i; else index=k;
3357  }
3358 
3359 
3360  public:
3362 
3366  complex& operator[](const int i) {
3367  if(row) {
3368 #if(CXSC_INDEX_CHECK)
3369  if(i<dat.A.lb2 || i>dat.A.ub2)
3370  cxscthrow(ELEMENT_NOT_IN_VEC("scmatrix_subv::operator[](int)"));
3371 #endif
3372  return dat.element(index,i);
3373  } else {
3374 #if(CXSC_INDEX_CHECK)
3375  if(i<dat.A.lb1 || i>dat.A.ub1)
3376  cxscthrow(ELEMENT_NOT_IN_VEC("scmatrix_subv::operator[](int)"));
3377 #endif
3378  return dat.element(i,index);
3379  }
3380  }
3381 
3383 
3386  const complex operator[](const int i) const {
3387  if(row) {
3388 #if(CXSC_INDEX_CHECK)
3389  if(i<dat.A.lb2 || i>dat.A.ub2)
3390  cxscthrow(ELEMENT_NOT_IN_VEC("scmatrix_subv::operator[](int)"));
3391 #endif
3392  return dat(index,i);
3393  } else {
3394 #if(CXSC_INDEX_CHECK)
3395  if(i<dat.A.lb1 || i>dat.A.ub1)
3396  cxscthrow(ELEMENT_NOT_IN_VEC("scmatrix_subv::operator[](int)"));
3397 #endif
3398  return dat(i,index);
3399  }
3400  }
3401 
3404  return sv_vs_assign(*this,v);
3405  }
3406 
3409  return sv_vs_assign(*this,v);
3410  }
3411 
3414  return svsp_vv_assign(*this,v);
3415  }
3416 
3419  return svsp_vv_assign(*this,v);
3420  }
3421 
3424  return svsl_vv_assign(*this,v);
3425  }
3426 
3429  return svsl_vv_assign(*this,v);
3430  }
3431 
3434  return svf_vv_assign(*this,v);
3435  }
3436 
3439  return svf_vv_assign(*this,v);
3440  }
3441 
3444  return svf_vv_assign(*this,v);
3445  }
3446 
3449  return svf_vv_assign(*this,v);
3450  }
3451 
3454  return svsp_vv_assign(*this,srvector(v));
3455  }
3456 
3459  return svsp_vv_assign(*this,scvector(v));
3460  }
3461 
3462 
3464  scmatrix_subv& operator*=(const real&);
3466  scmatrix_subv& operator*=(const complex&);
3468  scmatrix_subv& operator/=(const real&);
3470  scmatrix_subv& operator/=(const complex&);
3476  scmatrix_subv& operator+=(const rvector&);
3484  scmatrix_subv& operator-=(const rvector&);
3492  scmatrix_subv& operator+=(const cvector&);
3500  scmatrix_subv& operator-=(const cvector&);
3503 
3504  friend scvector operator-(const scmatrix_subv&);
3505  friend std::istream& operator>>(std::istream&, scmatrix_subv&);
3506 
3507  friend int Lb(const scmatrix_subv&);
3508  friend int Ub(const scmatrix_subv&);
3509  friend int VecLen(const scmatrix_subv&);
3510  friend srvector Re(const scmatrix_subv&);
3511  friend srvector Im(const scmatrix_subv&);
3512 
3513  friend class srvector;
3514  friend class srmatrix;
3515  friend class srmatrix_slice;
3516  friend class scvector;
3517  friend class scmatrix;
3518  friend class scmatrix_slice;
3519  friend class scivector;
3520  friend class scimatrix;
3521  friend class scimatrix_slice;
3522 
3523 #include "vector_friend_declarations.inl"
3524 };
3525 
3527 inline int Lb(const scmatrix_subv& S) {
3528  if(S.row)
3529  return Lb(S.dat, 2);
3530  else
3531  return Lb(S.dat, 1);
3532 }
3533 
3535 inline int Ub(const scmatrix_subv& S) {
3536  if(S.row)
3537  return Ub(S.dat, 2);
3538  else
3539  return Ub(S.dat, 1);
3540 }
3541 
3543 inline int VecLen(const scmatrix_subv& S) {
3544  return Ub(S)-Lb(S)+1;
3545 }
3546 
3548 inline srvector Re(const scmatrix_subv& S) {
3549  return Re(scvector(S));
3550 }
3551 
3553 inline srvector Im(const scmatrix_subv& S) {
3554  return Im(scvector(S));
3555 }
3556 
3558 inline std::ostream& operator<<(std::ostream& os, const scmatrix_subv& v) {
3559  os << scvector(v);
3560  return os;
3561 }
3562 
3564 inline std::istream& operator>>(std::istream& is, scmatrix_subv& v) {
3565  int n=0;
3566  if(v.row) n=v.dat.A.n; else n=v.dat.A.m;
3567  scvector tmp(n);
3568  is >> tmp;
3569  v = tmp;
3570  return is;
3571 }
3572 
3573 inline scmatrix_subv scmatrix::operator[](const cxscmatrix_column& c) {
3574 #if(CXSC_INDEX_CHECK)
3575  if(c.col()<lb2 || c.col()>ub2)
3576  cxscthrow(ROW_OR_COL_NOT_IN_MAT("scmatrix::operator[](const cxscmatrix_column&)"));
3577 #endif
3578  return scmatrix_subv(*this, false, lb1, ub1, c.col(), c.col());
3579 }
3580 
3582 #if(CXSC_INDEX_CHECK)
3583  if(i<lb1 || i>ub1)
3584  cxscthrow(ROW_OR_COL_NOT_IN_MAT("scmatrix::operator[](const int)"));
3585 #endif
3586  return scmatrix_subv(*this, true, i, i, lb2, ub2);
3587 }
3588 
3589 inline const scmatrix_subv scmatrix::operator[](const cxscmatrix_column& c) const {
3590 #if(CXSC_INDEX_CHECK)
3591  if(c.col()<lb2 || c.col()>ub2)
3592  cxscthrow(ROW_OR_COL_NOT_IN_MAT("scmatrix::operator[](const cxscmatrix_column&)"));
3593 #endif
3594  return scmatrix_subv(*this, false, lb1, ub1, c.col(), c.col());
3595 }
3596 
3597 inline const scmatrix_subv scmatrix::operator[](const int i) const {
3598 #if(CXSC_INDEX_CHECK)
3599  if(i<lb1 || i>ub1)
3600  cxscthrow(ROW_OR_COL_NOT_IN_MAT("scmatrix::operator[](const int)"));
3601 #endif
3602  return scmatrix_subv(*this, true, i, i, lb2, ub2);
3603 }
3604 
3606 #if(CXSC_INDEX_CHECK)
3607  if(i<A.lb1 || i>A.ub1)
3608  cxscthrow(ROW_OR_COL_NOT_IN_MAT("scmatrix_slice::operator[](const int)"));
3609 #endif
3610  return scmatrix_subv(*M, true, i, i, A.lb2, A.ub2);
3611 }
3612 
3613 inline scmatrix_subv scmatrix_slice::operator[](const cxscmatrix_column& c) {
3614 #if(CXSC_INDEX_CHECK)
3615  if(c.col()<A.lb2 || c.col()>A.ub2)
3616  cxscthrow(ROW_OR_COL_NOT_IN_MAT("scmatrix_slice::operator[](const cxscmatrix_column&)"));
3617 #endif
3618  return scmatrix_subv(*M, false, A.lb1, A.ub1, c.col(), c.col());
3619 }
3620 
3621 inline const scmatrix_subv scmatrix_slice::operator[](const int i) const {
3622 #if(CXSC_INDEX_CHECK)
3623  if(i<A.lb1 || i>A.ub1)
3624  cxscthrow(ROW_OR_COL_NOT_IN_MAT("scmatrix_slice::operator[](const int)"));
3625 #endif
3626  return scmatrix_subv(*M, true, i, i, A.lb2, A.ub2);
3627 }
3628 
3629 inline const scmatrix_subv scmatrix_slice::operator[](const cxscmatrix_column& c) const {
3630 #if(CXSC_INDEX_CHECK)
3631  if(c.col()<A.lb2 || c.col()>A.ub2)
3632  cxscthrow(ROW_OR_COL_NOT_IN_MAT("scmatrix_slice::operator[](const cxscmatrix_column&)"));
3633 #endif
3634  return scmatrix_subv(*M, false, A.lb1, A.ub1, c.col(), c.col());
3635 }
3636 
3638  int nnz = A.dat.A.get_nnz();
3639  p.reserve(nnz);
3640  x.reserve(nnz);
3641 
3642  if(A.row) {
3643  lb = A.dat.A.lb2;
3644  ub = A.dat.A.ub2;
3645  n = ub-lb+1;
3646 
3647  for(int j=0 ; j<n ; j++) {
3648  for(int k=A.dat.A.p[j] ; k<A.dat.A.p[j+1] ; k++) {
3649  p.push_back(j);
3650  x.push_back(A.dat.A.x[k]);
3651  }
3652  }
3653 
3654  } else {
3655  lb = A.dat.A.lb1;
3656  ub = A.dat.A.ub1;
3657  n = ub-lb+1;
3658 
3659  for(unsigned int k=0 ; k<A.dat.A.ind.size() ; k++) {
3660  p.push_back(A.dat.A.ind[k]);
3661  x.push_back(A.dat.A.x[k]);
3662  }
3663  }
3664 }
3665 
3667 inline scvector operator-(const scmatrix_subv& v) {
3668  scvector s(v);
3669  return -s;
3670 }
3671 
3673 inline scvector operator/(const scmatrix_subv& v1, const real& v2) {
3674  return scvector(v1) / v2;
3675 }
3676 
3678 inline scvector operator/(const scmatrix_subv& v1, const complex& v2) {
3679  return scvector(v1) / v2;
3680 }
3681 
3683 inline scvector operator/(const srmatrix_subv& v1, const complex& v2) {
3684  return srvector(v1) / v2;
3685 }
3686 
3688 inline scvector operator*(const scmatrix_subv& v1, const real& v2) {
3689  return scvector(v1) * v2;
3690 }
3691 
3693 inline scvector operator*(const scmatrix_subv& v1, const complex& v2) {
3694  return scvector(v1) * v2;
3695 }
3696 
3698 inline scvector operator*(const srmatrix_subv& v1, const complex& v2) {
3699  return srvector(v1) * v2;
3700 }
3701 
3703 inline scvector operator*(const real& v1, const scmatrix_subv& v2) {
3704  return v1 * scvector(v2);
3705 }
3706 
3708 inline scvector operator*(const complex& v1, const scmatrix_subv& v2) {
3709  return v1 * scvector(v2);
3710 }
3711 
3713 inline scvector operator*(const complex& v1, const srmatrix_subv& v2) {
3714  return v1 * srvector(v2);
3715 }
3716 
3718 
3724 inline complex operator*(const scmatrix_subv& v1, const srvector& v2) {
3725  return scvector(v1) * v2;
3726 }
3727 
3729 
3735 inline complex operator*(const srmatrix_subv& v1, const scvector& v2) {
3736  return srvector(v1) * v2;
3737 }
3738 
3740 
3746 inline complex operator*(const scmatrix_subv& v1, const scvector& v2) {
3747  return scvector(v1) * v2;
3748 }
3749 
3751 
3757 inline complex operator*(const scmatrix_subv& v1, const srvector_slice& v2) {
3758  return scvector(v1) * v2;
3759 }
3760 
3762 
3768 inline complex operator*(const srmatrix_subv& v1, const scvector_slice& v2) {
3769  return srvector(v1) * v2;
3770 }
3771 
3773 
3779 inline complex operator*(const scmatrix_subv& v1, const scvector_slice& v2) {
3780  return scvector(v1) * v2;
3781 }
3782 
3784 
3790 inline complex operator*(const scmatrix_subv& v1, const rvector& v2) {
3791  return scvector(v1) * v2;
3792 }
3793 
3795 
3801 inline complex operator*(const srmatrix_subv& v1, const cvector& v2) {
3802  return srvector(v1) * v2;
3803 }
3804 
3806 
3812 inline complex operator*(const scmatrix_subv& v1, const cvector& v2) {
3813  return scvector(v1) * v2;
3814 }
3815 
3817 
3823 inline complex operator*(const scmatrix_subv& v1, const rvector_slice& v2) {
3824  return scvector(v1) * v2;
3825 }
3826 
3828 
3834 inline complex operator*(const srmatrix_subv& v1, const cvector_slice& v2) {
3835  return srvector(v1) * v2;
3836 }
3837 
3839 
3845 inline complex operator*(const scmatrix_subv& v1, const cvector_slice& v2) {
3846  return scvector(v1) * v2;
3847 }
3848 
3850 
3856 inline complex operator*(const scvector& v1, const srmatrix_subv& v2) {
3857  return v1 * srvector(v2);
3858 }
3859 
3861 
3867 inline complex operator*(const srvector& v1, const scmatrix_subv& v2) {
3868  return v1 * scvector(v2);
3869 }
3870 
3872 
3878 inline complex operator*(const scvector& v1, const scmatrix_subv& v2) {
3879  return v1 * scvector(v2);
3880 }
3881 
3883 
3889 inline complex operator*(const scvector_slice& v1, const srmatrix_subv& v2) {
3890  return v1 * srvector(v2);
3891 }
3892 
3894 
3900 inline complex operator*(const srvector_slice& v1, const scmatrix_subv& v2) {
3901  return v1 * scvector(v2);
3902 }
3903 
3905 
3911 inline complex operator*(const scvector_slice& v1, const scmatrix_subv& v2) {
3912  return v1 * scvector(v2);
3913 }
3914 
3916 
3922 inline complex operator*(const cvector& v1, const srmatrix_subv& v2) {
3923  return v1 * srvector(v2);
3924 }
3925 
3927 
3933 inline complex operator*(const rvector& v1, const scmatrix_subv& v2) {
3934  return v1 * scvector(v2);
3935 }
3936 
3938 
3944 inline complex operator*(const cvector& v1, const scmatrix_subv& v2) {
3945  return v1 * scvector(v2);
3946 }
3947 
3949 
3955 inline complex operator*(const cvector_slice& v1, const srmatrix_subv& v2) {
3956  return v1 * srvector(v2);
3957 }
3958 
3960 
3966 inline complex operator*(const rvector_slice& v1, const scmatrix_subv& v2) {
3967  return v1 * scvector(v2);
3968 }
3969 
3971 
3977 inline complex operator*(const cvector_slice& v1, const scmatrix_subv& v2) {
3978  return v1 * scvector(v2);
3979 }
3980 
3982 inline scvector operator+(const scmatrix_subv& v1, const srvector& v2) {
3983  return scvector(v1) + v2;
3984 }
3985 
3987 inline scvector operator+(const srmatrix_subv& v1, const scvector& v2) {
3988  return srvector(v1) + v2;
3989 }
3990 
3992 inline scvector operator+(const scmatrix_subv& v1, const scvector& v2) {
3993  return scvector(v1) + v2;
3994 }
3995 
3997 inline scvector operator+(const scmatrix_subv& v1, const srvector_slice& v2) {
3998  return scvector(v1) + v2;
3999 }
4000 
4002 inline scvector operator+(const srmatrix_subv& v1, const scvector_slice& v2) {
4003  return srvector(v1) + v2;
4004 }
4005 
4007 inline scvector operator+(const scmatrix_subv& v1, const scvector_slice& v2) {
4008  return scvector(v1) + v2;
4009 }
4010 
4012 inline cvector operator+(const scmatrix_subv& v1, const rvector& v2) {
4013  return scvector(v1) + v2;
4014 }
4015 
4017 inline cvector operator+(const srmatrix_subv& v1, const cvector& v2) {
4018  return srvector(v1) + v2;
4019 }
4020 
4022 inline cvector operator+(const scmatrix_subv& v1, const cvector& v2) {
4023  return scvector(v1) + v2;
4024 }
4025 
4027 inline cvector operator+(const scmatrix_subv& v1, const rvector_slice& v2) {
4028  return scvector(v1) + v2;
4029 }
4030 
4032 inline cvector operator+(const srmatrix_subv& v1, const cvector_slice& v2) {
4033  return srvector(v1) + v2;
4034 }
4035 
4037 inline cvector operator+(const scmatrix_subv& v1, const cvector_slice& v2) {
4038  return scvector(v1) + v2;
4039 }
4040 
4042 inline scvector operator+(const scvector& v1, const srmatrix_subv& v2) {
4043  return v1 + srvector(v2);
4044 }
4045 
4047 inline scvector operator+(const srvector& v1, const scmatrix_subv& v2) {
4048  return v1 + scvector(v2);
4049 }
4050 
4052 inline scvector operator+(const scvector& v1, const scmatrix_subv& v2) {
4053  return v1 + scvector(v2);
4054 }
4055 
4057 inline scvector operator+(const scvector_slice& v1, const srmatrix_subv& v2) {
4058  return v1 + srvector(v2);
4059 }
4060 
4062 inline scvector operator+(const srvector_slice& v1, const scmatrix_subv& v2) {
4063  return v1 + scvector(v2);
4064 }
4065 
4067 inline scvector operator+(const scvector_slice& v1, const scmatrix_subv& v2) {
4068  return v1 + scvector(v2);
4069 }
4070 
4072 inline cvector operator+(const cvector& v1, const srmatrix_subv& v2) {
4073  return v1 + srvector(v2);
4074 }
4075 
4077 inline cvector operator+(const rvector& v1, const scmatrix_subv& v2) {
4078  return v1 + scvector(v2);
4079 }
4080 
4082 inline cvector operator+(const cvector& v1, const scmatrix_subv& v2) {
4083  return v1 + scvector(v2);
4084 }
4085 
4087 inline cvector operator+(const cvector_slice& v1, const srmatrix_subv& v2) {
4088  return v1 + srvector(v2);
4089 }
4090 
4092 inline cvector operator+(const rvector_slice& v1, const scmatrix_subv& v2) {
4093  return v1 + scvector(v2);
4094 }
4095 
4097 inline cvector operator+(const cvector_slice& v1, const scmatrix_subv& v2) {
4098  return v1 + scvector(v2);
4099 }
4100 
4102 inline scvector operator-(const scmatrix_subv& v1, const srvector& v2) {
4103  return scvector(v1) - v2;
4104 }
4105 
4107 inline scvector operator-(const srmatrix_subv& v1, const scvector& v2) {
4108  return srvector(v1) - v2;
4109 }
4110 
4112 inline scvector operator-(const scmatrix_subv& v1, const scvector& v2) {
4113  return scvector(v1) - v2;
4114 }
4115 
4117 inline scvector operator-(const scmatrix_subv& v1, const srvector_slice& v2) {
4118  return scvector(v1) - v2;
4119 }
4120 
4122 inline scvector operator-(const srmatrix_subv& v1, const scvector_slice& v2) {
4123  return srvector(v1) - v2;
4124 }
4125 
4127 inline scvector operator-(const scmatrix_subv& v1, const scvector_slice& v2) {
4128  return scvector(v1) - v2;
4129 }
4130 
4132 inline cvector operator-(const scmatrix_subv& v1, const rvector& v2) {
4133  return scvector(v1) - v2;
4134 }
4135 
4137 inline cvector operator-(const srmatrix_subv& v1, const cvector& v2) {
4138  return srvector(v1) - v2;
4139 }
4140 
4142 inline cvector operator-(const scmatrix_subv& v1, const cvector& v2) {
4143  return scvector(v1) - v2;
4144 }
4145 
4147 inline cvector operator-(const scmatrix_subv& v1, const rvector_slice& v2) {
4148  return scvector(v1) - v2;
4149 }
4150 
4152 inline cvector operator-(const srmatrix_subv& v1, const cvector_slice& v2) {
4153  return srvector(v1) - v2;
4154 }
4155 
4157 inline cvector operator-(const scmatrix_subv& v1, const cvector_slice& v2) {
4158  return scvector(v1) - v2;
4159 }
4160 
4162 inline scvector operator-(const scvector& v1, const srmatrix_subv& v2) {
4163  return v1 - srvector(v2);
4164 }
4165 
4167 inline scvector operator-(const srvector& v1, const scmatrix_subv& v2) {
4168  return v1 - scvector(v2);
4169 }
4170 
4172 inline scvector operator-(const scvector& v1, const scmatrix_subv& v2) {
4173  return v1 - scvector(v2);
4174 }
4175 
4177 inline scvector operator-(const scvector_slice& v1, const srmatrix_subv& v2) {
4178  return v1 - srvector(v2);
4179 }
4180 
4182 inline scvector operator-(const srvector_slice& v1, const scmatrix_subv& v2) {
4183  return v1 - scvector(v2);
4184 }
4185 
4187 inline scvector operator-(const scvector_slice& v1, const scmatrix_subv& v2) {
4188  return v1 - scvector(v2);
4189 }
4190 
4192 inline cvector operator-(const cvector& v1, const srmatrix_subv& v2) {
4193  return v1 - srvector(v2);
4194 }
4195 
4197 inline cvector operator-(const rvector& v1, const scmatrix_subv& v2) {
4198  return v1 - scvector(v2);
4199 }
4200 
4202 inline cvector operator-(const cvector& v1, const scmatrix_subv& v2) {
4203  return v1 - scvector(v2);
4204 }
4205 
4207 inline cvector operator-(const cvector_slice& v1, const srmatrix_subv& v2) {
4208  return v1 - srvector(v2);
4209 }
4210 
4212 inline cvector operator-(const rvector_slice& v1, const scmatrix_subv& v2) {
4213  return v1 - scvector(v2);
4214 }
4215 
4217 inline cvector operator-(const cvector_slice& v1, const scmatrix_subv& v2) {
4218  return v1 - scvector(v2);
4219 }
4220 
4222  *this = *this * v;
4223  return *this;
4224 }
4225 
4227  *this = *this * v;
4228  return *this;
4229 }
4230 
4232  *this = *this / v;
4233  return *this;
4234 }
4235 
4237  *this = *this / v;
4238  return *this;
4239 }
4240 
4242  *this = *this + v;
4243  return *this;
4244 }
4245 
4247  *this = *this + v;
4248  return *this;
4249 }
4250 
4252  *this = *this + v;
4253  return *this;
4254 }
4255 
4257  *this = *this + v;
4258  return *this;
4259 }
4260 
4262  *this = *this - v;
4263  return *this;
4264 }
4265 
4267  *this = *this - v;
4268  return *this;
4269 }
4270 
4272  *this = *this - v;
4273  return *this;
4274 }
4275 
4277  *this = *this - v;
4278  return *this;
4279 }
4280 
4282  *this = *this + v;
4283  return *this;
4284 }
4285 
4287  *this = *this + v;
4288  return *this;
4289 }
4290 
4292  *this = *this + v;
4293  return *this;
4294 }
4295 
4297  *this = *this + v;
4298  return *this;
4299 }
4300 
4302  *this = *this - v;
4303  return *this;
4304 }
4305 
4307  *this = *this - v;
4308  return *this;
4309 }
4310 
4312  *this = *this - v;
4313  return *this;
4314 }
4315 
4317  *this = *this - v;
4318  return *this;
4319 }
4320 
4322  *this += rvector(v);
4323  return *this;
4324 }
4325 
4327  *this += cvector(v);
4328  return *this;
4329 }
4330 
4332  *this += rvector(v);
4333  return *this;
4334 }
4335 
4337  *this += cvector(v);
4338  return *this;
4339 }
4340 
4342  *this += rvector(v);
4343  return *this;
4344 }
4345 
4347  *this += cvector(v);
4348  return *this;
4349 }
4350 
4351 inline cmatrix_subv& cmatrix_subv::operator-=(const srmatrix_subv& v) {
4352  *this -= rvector(v);
4353  return *this;
4354 }
4355 
4356 inline cmatrix_subv& cmatrix_subv::operator-=(const scmatrix_subv& v) {
4357  *this -= cvector(v);
4358  return *this;
4359 }
4360 
4362  *this = rvector(v);
4363  return *this;
4364 }
4365 
4367  *this = cvector(v);
4368  return *this;
4369 }
4370 
4372  *this = rvector(v);
4373  return *this;
4374 }
4375 
4377  *this = cvector(v);
4378  return *this;
4379 }
4380 
4382  *this = rvector(v);
4383  return *this;
4384 }
4385 
4387  *this = cvector(v);
4388  return *this;
4389 }
4390 
4392 inline bool operator==(const scmatrix_subv& v1, const srvector& v2) {
4393  return scvector(v1) == v2;
4394 }
4395 
4397 inline bool operator==(const srmatrix_subv& v1, const scvector& v2) {
4398  return srvector(v1) == v2;
4399 }
4400 
4402 inline bool operator==(const scmatrix_subv& v1, const scvector& v2) {
4403  return scvector(v1) == v2;
4404 }
4405 
4407 inline bool operator==(const scmatrix_subv& v1, const srvector_slice& v2) {
4408  return scvector(v1) == v2;
4409 }
4410 
4412 inline bool operator==(const srmatrix_subv& v1, const scvector_slice& v2) {
4413  return srvector(v1) == v2;
4414 }
4415 
4417 inline bool operator==(const scmatrix_subv& v1, const scvector_slice& v2) {
4418  return scvector(v1) == v2;
4419 }
4420 
4422 inline bool operator==(const scmatrix_subv& v1, const rvector& v2) {
4423  return scvector(v1) == v2;
4424 }
4425 
4427 inline bool operator==(const srmatrix_subv& v1, const cvector& v2) {
4428  return srvector(v1) == v2;
4429 }
4430 
4432 inline bool operator==(const scmatrix_subv& v1, const cvector& v2) {
4433  return scvector(v1) == v2;
4434 }
4435 
4437 inline bool operator==(const scmatrix_subv& v1, const rvector_slice& v2) {
4438  return scvector(v1) == v2;
4439 }
4440 
4442 inline bool operator==(const srmatrix_subv& v1, const cvector_slice& v2) {
4443  return srvector(v1) == v2;
4444 }
4445 
4447 inline bool operator==(const scmatrix_subv& v1, const cvector_slice& v2) {
4448  return scvector(v1) == v2;
4449 }
4450 
4452 inline bool operator==(const scvector& v1, const srmatrix_subv& v2) {
4453  return v1 == srvector(v2);
4454 }
4455 
4457 inline bool operator==(const srvector& v1, const scmatrix_subv& v2) {
4458  return v1 == scvector(v2);
4459 }
4460 
4462 inline bool operator==(const scvector& v1, const scmatrix_subv& v2) {
4463  return v1 == scvector(v2);
4464 }
4465 
4467 inline bool operator==(const scvector_slice& v1, const srmatrix_subv& v2) {
4468  return v1 == srvector(v2);
4469 }
4470 
4472 inline bool operator==(const srvector_slice& v1, const scmatrix_subv& v2) {
4473  return v1 == scvector(v2);
4474 }
4475 
4477 inline bool operator==(const scvector_slice& v1, const scmatrix_subv& v2) {
4478  return v1 == scvector(v2);
4479 }
4480 
4482 inline bool operator==(const cvector& v1, const srmatrix_subv& v2) {
4483  return v1 == srvector(v2);
4484 }
4485 
4487 inline bool operator==(const rvector& v1, const scmatrix_subv& v2) {
4488  return v1 == scvector(v2);
4489 }
4490 
4492 inline bool operator==(const cvector& v1, const scmatrix_subv& v2) {
4493  return v1 == scvector(v2);
4494 }
4495 
4497 inline bool operator==(const cvector_slice& v1, const srmatrix_subv& v2) {
4498  return v1 == srvector(v2);
4499 }
4500 
4502 inline bool operator==(const rvector_slice& v1, const scmatrix_subv& v2) {
4503  return v1 == scvector(v2);
4504 }
4505 
4507 inline bool operator==(const cvector_slice& v1, const scmatrix_subv& v2) {
4508  return v1 == scvector(v2);
4509 }
4510 
4512 inline bool operator!=(const scmatrix_subv& v1, const srvector& v2) {
4513  return scvector(v1) != v2;
4514 }
4515 
4517 inline bool operator!=(const srmatrix_subv& v1, const scvector& v2) {
4518  return srvector(v1) != v2;
4519 }
4520 
4522 inline bool operator!=(const scmatrix_subv& v1, const scvector& v2) {
4523  return scvector(v1) != v2;
4524 }
4525 
4527 inline bool operator!=(const scmatrix_subv& v1, const srvector_slice& v2) {
4528  return scvector(v1) != v2;
4529 }
4530 
4532 inline bool operator!=(const srmatrix_subv& v1, const scvector_slice& v2) {
4533  return srvector(v1) != v2;
4534 }
4535 
4537 inline bool operator!=(const scmatrix_subv& v1, const scvector_slice& v2) {
4538  return scvector(v1) != v2;
4539 }
4540 
4542 inline bool operator!=(const scmatrix_subv& v1, const rvector& v2) {
4543  return scvector(v1) != v2;
4544 }
4545 
4547 inline bool operator!=(const srmatrix_subv& v1, const cvector& v2) {
4548  return srvector(v1) != v2;
4549 }
4550 
4552 inline bool operator!=(const scmatrix_subv& v1, const cvector& v2) {
4553  return scvector(v1) != v2;
4554 }
4555 
4557 inline bool operator!=(const scmatrix_subv& v1, const rvector_slice& v2) {
4558  return scvector(v1) != v2;
4559 }
4561 
4562 inline bool operator!=(const srmatrix_subv& v1, const cvector_slice& v2) {
4563  return srvector(v1) != v2;
4564 }
4565 
4567 inline bool operator!=(const scmatrix_subv& v1, const cvector_slice& v2) {
4568  return scvector(v1) != v2;
4569 }
4570 
4572 inline bool operator!=(const scvector& v1, const srmatrix_subv& v2) {
4573  return v1 != srvector(v2);
4574 }
4575 
4577 inline bool operator!=(const srvector& v1, const scmatrix_subv& v2) {
4578  return v1 != scvector(v2);
4579 }
4580 
4582 inline bool operator!=(const scvector& v1, const scmatrix_subv& v2) {
4583  return v1 != scvector(v2);
4584 }
4585 
4587 inline bool operator!=(const scvector_slice& v1, const srmatrix_subv& v2) {
4588  return v1 != srvector(v2);
4589 }
4590 
4592 inline bool operator!=(const srvector_slice& v1, const scmatrix_subv& v2) {
4593  return v1 != scvector(v2);
4594 }
4595 
4597 inline bool operator!=(const scvector_slice& v1, const scmatrix_subv& v2) {
4598  return v1 != scvector(v2);
4599 }
4600 
4602 inline bool operator!=(const cvector& v1, const srmatrix_subv& v2) {
4603  return v1 != srvector(v2);
4604 }
4605 
4607 inline bool operator!=(const rvector& v1, const scmatrix_subv& v2) {
4608  return v1 != scvector(v2);
4609 }
4610 
4612 inline bool operator!=(const cvector& v1, const scmatrix_subv& v2) {
4613  return v1 != scvector(v2);
4614 }
4615 
4617 inline bool operator!=(const cvector_slice& v1, const srmatrix_subv& v2) {
4618  return v1 != srvector(v2);
4619 }
4620 
4622 inline bool operator!=(const rvector_slice& v1, const scmatrix_subv& v2) {
4623  return v1 != scvector(v2);
4624 }
4625 
4627 inline bool operator!=(const cvector_slice& v1, const scmatrix_subv& v2) {
4628  return v1 != scvector(v2);
4629 }
4630 
4632 inline bool operator!(const scmatrix_subv& x) {
4633  return sv_v_not(x);
4634 }
4635 
4637 
4640 inline void accumulate(cdotprecision& dot, const scmatrix_subv& v1, const scmatrix_subv& v2) {
4641  spsp_vv_accu<cdotprecision,scvector,scvector,sparse_cdot>(dot, scvector(v1), scvector(v2));
4642 }
4643 
4645 
4648 inline void accumulate(cdotprecision& dot, const scmatrix_subv& v1, const srmatrix_subv& v2) {
4649  spsp_vv_accu<cdotprecision,scvector,srvector,sparse_cdot>(dot, scvector(v1), srvector(v2));
4650 }
4651 
4653 
4656 inline void accumulate(cdotprecision& dot, const srmatrix_subv& v1, const scmatrix_subv& v2) {
4657  spsp_vv_accu<cdotprecision,srvector,scvector,sparse_cdot>(dot, srvector(v1), scvector(v2));
4658 }
4659 
4661 
4664 inline void accumulate(cdotprecision& dot, const scmatrix_subv& v1, const scvector& v2) {
4665  spsp_vv_accu<cdotprecision,scvector,scvector,sparse_cdot>(dot, scvector(v1), v2);
4666 }
4667 
4669 
4672 inline void accumulate(cdotprecision& dot, const scmatrix_subv& v1, const srvector& v2) {
4673  spsp_vv_accu<cdotprecision,scvector,srvector,sparse_cdot>(dot, scvector(v1), v2);
4674 }
4675 
4677 
4680 inline void accumulate(cdotprecision& dot, const srmatrix_subv& v1, const scvector& v2) {
4681  spsp_vv_accu<cdotprecision,srvector,scvector,sparse_cdot>(dot, srvector(v1), v2);
4682 }
4683 
4685 
4688 inline void accumulate(cdotprecision& dot, const scmatrix_subv& v1, const scvector_slice& v2) {
4689  spsl_vv_accu<cdotprecision,scvector,scvector_slice,sparse_cdot>(dot, scvector(v1), v2);
4690 }
4691 
4693 
4696 inline void accumulate(cdotprecision& dot, const scmatrix_subv& v1, const srvector_slice& v2) {
4697  spsl_vv_accu<cdotprecision,scvector,srvector_slice,sparse_cdot>(dot, scvector(v1), v2);
4698 }
4699 
4701 
4704 inline void accumulate(cdotprecision& dot, const srmatrix_subv& v1, const scvector_slice& v2) {
4705  spsl_vv_accu<cdotprecision,srvector,scvector_slice,sparse_cdot>(dot, srvector(v1), v2);
4706 }
4707 
4709 
4712 inline void accumulate(cdotprecision& dot, const scmatrix_subv& v1, const cvector& v2) {
4713  spf_vv_accu<cdotprecision,scvector,cvector,sparse_cdot>(dot, scvector(v1), v2);
4714 }
4715 
4717 
4720 inline void accumulate(cdotprecision& dot, const scmatrix_subv& v1, const rvector& v2) {
4721  spf_vv_accu<cdotprecision,scvector,rvector,sparse_cdot>(dot, scvector(v1), v2);
4722 }
4723 
4725 
4728 inline void accumulate(cdotprecision& dot, const srmatrix_subv& v1, const cvector& v2) {
4729  spf_vv_accu<cdotprecision,srvector,cvector,sparse_cdot>(dot, srvector(v1), v2);
4730 }
4731 
4733 
4736 inline void accumulate(cdotprecision& dot, const scmatrix_subv& v1, const cvector_slice& v2) {
4737  spf_vv_accu<cdotprecision,scvector,cvector_slice,sparse_cdot>(dot, scvector(v1), v2);
4738 }
4739 
4741 
4744 inline void accumulate(cdotprecision& dot, const scmatrix_subv& v1, const rvector_slice& v2) {
4745  spf_vv_accu<cdotprecision,scvector,rvector_slice,sparse_cdot>(dot, scvector(v1), v2);
4746 }
4747 
4749 
4752 inline void accumulate(cdotprecision& dot, const srmatrix_subv& v1, const cvector_slice& v2) {
4753  spf_vv_accu<cdotprecision,srvector,cvector_slice,sparse_cdot>(dot, srvector(v1), v2);
4754 }
4755 
4757 
4760 inline void accumulate(cdotprecision& dot, const scvector& v1, const scmatrix_subv& v2) {
4761  spsp_vv_accu<cdotprecision,scvector,scvector,sparse_cdot>(dot, v1, scvector(v2));
4762 }
4763 
4765 
4768 inline void accumulate(cdotprecision& dot, const scvector& v1, const srmatrix_subv& v2) {
4769  spsp_vv_accu<cdotprecision,scvector,srvector,sparse_cdot>(dot, v1, srvector(v2));
4770 }
4771 
4773 
4776 inline void accumulate(cdotprecision& dot, const srvector& v1, const scmatrix_subv& v2) {
4777  spsp_vv_accu<cdotprecision,srvector,scvector,sparse_cdot>(dot, v1, scvector(v2));
4778 }
4779 
4781 
4784 inline void accumulate(cdotprecision& dot, const scvector_slice& v1, const scmatrix_subv& v2) {
4785  slsp_vv_accu<cdotprecision,scvector_slice,scvector,sparse_cdot>(dot, v1, scvector(v2));
4786 }
4787 
4789 
4792 inline void accumulate(cdotprecision& dot, const scvector_slice& v1, const srmatrix_subv& v2) {
4793  slsp_vv_accu<cdotprecision,scvector_slice,srvector,sparse_cdot>(dot, v1, srvector(v2));
4794 }
4795 
4797 
4800 inline void accumulate(cdotprecision& dot, const srvector_slice& v1, const scmatrix_subv& v2) {
4801  slsp_vv_accu<cdotprecision,srvector_slice,scvector,sparse_cdot>(dot, v1, scvector(v2));
4802 }
4803 
4805 
4808 inline void accumulate(cdotprecision& dot, const cvector& v1, const scmatrix_subv& v2) {
4809  fsp_vv_accu<cdotprecision,cvector,scvector,sparse_cdot>(dot, v1, scvector(v2));
4810 }
4811 
4813 
4816 inline void accumulate(cdotprecision& dot, const cvector& v1, const srmatrix_subv& v2) {
4817  fsp_vv_accu<cdotprecision,cvector,srvector,sparse_cdot>(dot, v1, srvector(v2));
4818 }
4819 
4821 
4824 inline void accumulate(cdotprecision& dot, const rvector& v1, const scmatrix_subv& v2) {
4825  fsp_vv_accu<cdotprecision,rvector,scvector,sparse_cdot>(dot, v1, scvector(v2));
4826 }
4827 
4829 
4832 inline void accumulate(cdotprecision& dot, const cvector_slice& v1, const scmatrix_subv& v2) {
4833  fsp_vv_accu<cdotprecision,cvector_slice,scvector,sparse_cdot>(dot, v1, scvector(v2));
4834 }
4835 
4837 
4840 inline void accumulate(cdotprecision& dot, const cvector_slice& v1, const srmatrix_subv& v2) {
4841  fsp_vv_accu<cdotprecision,cvector_slice,srvector,sparse_cdot>(dot, v1, srvector(v2));
4842 }
4843 
4845 
4848 inline void accumulate(cdotprecision& dot, const rvector_slice& v1, const scmatrix_subv& v2) {
4849  fsp_vv_accu<cdotprecision,rvector_slice,scvector,sparse_cdot>(dot, v1, scvector(v2));
4850 }
4851 
4853 
4858 inline void accumulate_approx(cdotprecision& dot, const scmatrix_subv& v1, const scmatrix_subv& v2) {
4859  spsp_vv_accuapprox<cdotprecision,scvector,scvector,sparse_cdot>(dot, scvector(v1), scvector(v2));
4860 }
4861 
4863 
4868 inline void accumulate_approx(cdotprecision& dot, const scmatrix_subv& v1, const srmatrix_subv& v2) {
4869  spsp_vv_accuapprox<cdotprecision,scvector,srvector,sparse_cdot>(dot, scvector(v1), srvector(v2));
4870 }
4871 
4873 
4878 inline void accumulate_approx(cdotprecision& dot, const srmatrix_subv& v1, const scmatrix_subv& v2) {
4879  spsp_vv_accuapprox<cdotprecision,srvector,scvector,sparse_cdot>(dot, srvector(v1), scvector(v2));
4880 }
4881 
4883 
4888 inline void accumulate_approx(cdotprecision& dot, const scmatrix_subv& v1, const scvector& v2) {
4889  spsp_vv_accuapprox<cdotprecision,scvector,scvector,sparse_cdot>(dot, scvector(v1), v2);
4890 }
4891 
4893 
4898 inline void accumulate_approx(cdotprecision& dot, const scmatrix_subv& v1, const srvector& v2) {
4899  spsp_vv_accuapprox<cdotprecision,scvector,srvector,sparse_cdot>(dot, scvector(v1), v2);
4900 }
4901 
4903 
4908 inline void accumulate_approx(cdotprecision& dot, const srmatrix_subv& v1, const scvector& v2) {
4909  spsp_vv_accuapprox<cdotprecision,srvector,scvector,sparse_cdot>(dot, srvector(v1), v2);
4910 }
4911 
4913 
4918 inline void accumulate_approx(cdotprecision& dot, const scmatrix_subv& v1, const scvector_slice& v2) {
4919  spsl_vv_accuapprox<cdotprecision,scvector,scvector_slice,sparse_cdot>(dot, scvector(v1), v2);
4920 }
4921 
4923 
4928 inline void accumulate_approx(cdotprecision& dot, const scmatrix_subv& v1, const srvector_slice& v2) {
4929  spsl_vv_accuapprox<cdotprecision,scvector,srvector_slice,sparse_cdot>(dot, scvector(v1), v2);
4930 }
4931 
4933 
4938 inline void accumulate_approx(cdotprecision& dot, const srmatrix_subv& v1, const scvector_slice& v2) {
4939  spsl_vv_accuapprox<cdotprecision,srvector,scvector_slice,sparse_cdot>(dot, srvector(v1), v2);
4940 }
4941 
4943 
4948 inline void accumulate_approx(cdotprecision& dot, const scmatrix_subv& v1, const cvector& v2) {
4949  spf_vv_accuapprox<cdotprecision,scvector,cvector,sparse_cdot>(dot, scvector(v1), v2);
4950 }
4951 
4953 
4958 inline void accumulate_approx(cdotprecision& dot, const scmatrix_subv& v1, const rvector& v2) {
4959  spf_vv_accuapprox<cdotprecision,scvector,rvector,sparse_cdot>(dot, scvector(v1), v2);
4960 }
4961 
4963 
4968 inline void accumulate_approx(cdotprecision& dot, const srmatrix_subv& v1, const cvector& v2) {
4969  spf_vv_accuapprox<cdotprecision,srvector,cvector,sparse_cdot>(dot, srvector(v1), v2);
4970 }
4971 
4973 
4978 inline void accumulate_approx(cdotprecision& dot, const scmatrix_subv& v1, const cvector_slice& v2) {
4979  spf_vv_accuapprox<cdotprecision,scvector,cvector_slice,sparse_cdot>(dot, scvector(v1), v2);
4980 }
4981 
4983 
4988 inline void accumulate_approx(cdotprecision& dot, const scmatrix_subv& v1, const rvector_slice& v2) {
4989  spf_vv_accuapprox<cdotprecision,scvector,rvector_slice,sparse_cdot>(dot, scvector(v1), v2);
4990 }
4991 
4993 
4998 inline void accumulate_approx(cdotprecision& dot, const srmatrix_subv& v1, const cvector_slice& v2) {
4999  spf_vv_accuapprox<cdotprecision,srvector,cvector_slice,sparse_cdot>(dot, srvector(v1), v2);
5000 }
5001 
5003 
5008 inline void accumulate_approx(cdotprecision& dot, const scvector& v1, const scmatrix_subv& v2) {
5009  spsp_vv_accuapprox<cdotprecision,scvector,scvector,sparse_cdot>(dot, v1, scvector(v2));
5010 }
5011 
5013 
5018 inline void accumulate_approx(cdotprecision& dot, const scvector& v1, const srmatrix_subv& v2) {
5019  spsp_vv_accuapprox<cdotprecision,scvector,srvector,sparse_cdot>(dot, v1, srvector(v2));
5020 }
5021 
5023 
5028 inline void accumulate_approx(cdotprecision& dot, const srvector& v1, const scmatrix_subv& v2) {
5029  spsp_vv_accuapprox<cdotprecision,srvector,scvector,sparse_cdot>(dot, v1, scvector(v2));
5030 }
5031 
5033 
5038 inline void accumulate_approx(cdotprecision& dot, const scvector_slice& v1, const scmatrix_subv& v2) {
5039  slsp_vv_accuapprox<cdotprecision,scvector_slice,scvector,sparse_cdot>(dot, v1, scvector(v2));
5040 }
5041 
5043 
5048 inline void accumulate_approx(cdotprecision& dot, const scvector_slice& v1, const srmatrix_subv& v2) {
5049  slsp_vv_accuapprox<cdotprecision,scvector_slice,srvector,sparse_cdot>(dot, v1, srvector(v2));
5050 }
5051 
5053 
5058 inline void accumulate_approx(cdotprecision& dot, const srvector_slice& v1, const scmatrix_subv& v2) {
5059  slsp_vv_accuapprox<cdotprecision,srvector_slice,scvector,sparse_cdot>(dot, v1, scvector(v2));
5060 }
5061 
5063 
5068 inline void accumulate_approx(cdotprecision& dot, const cvector& v1, const scmatrix_subv& v2) {
5069  fsp_vv_accuapprox<cdotprecision,cvector,scvector,sparse_cdot>(dot, v1, scvector(v2));
5070 }
5071 
5073 
5078 inline void accumulate_approx(cdotprecision& dot, const cvector& v1, const srmatrix_subv& v2) {
5079  fsp_vv_accuapprox<cdotprecision,cvector,srvector,sparse_cdot>(dot, v1, srvector(v2));
5080 }
5081 
5083 
5088 inline void accumulate_approx(cdotprecision& dot, const rvector& v1, const scmatrix_subv& v2) {
5089  fsp_vv_accuapprox<cdotprecision,rvector,scvector,sparse_cdot>(dot, v1, scvector(v2));
5090 }
5091 
5093 
5098 inline void accumulate_approx(cdotprecision& dot, const cvector_slice& v1, const scmatrix_subv& v2) {
5099  fsp_vv_accuapprox<cdotprecision,cvector_slice,scvector,sparse_cdot>(dot, v1, scvector(v2));
5100 }
5101 
5103 
5108 inline void accumulate_approx(cdotprecision& dot, const cvector_slice& v1, const srmatrix_subv& v2) {
5109  fsp_vv_accuapprox<cdotprecision,cvector_slice,srvector,sparse_cdot>(dot, v1, srvector(v2));
5110 }
5111 
5113 
5118 inline void accumulate_approx(cdotprecision& dot, const rvector_slice& v1, const scmatrix_subv& v2) {
5119  fsp_vv_accuapprox<cdotprecision,rvector_slice,scvector,sparse_cdot>(dot, v1, scvector(v2));
5120 }
5121 
5123 
5126 inline void accumulate(cidotprecision& dot, const scmatrix_subv& v1, const scmatrix_subv& v2) {
5127  cdotprecision tmp(0.0);
5128  tmp.set_k(dot.get_k());
5129  accumulate(tmp,scvector(v1),scvector(v2));
5130  dot += tmp;
5131 }
5132 
5134 
5137 inline void accumulate(cidotprecision& dot, const scmatrix_subv& v1, const srmatrix_subv& v2) {
5138  cdotprecision tmp(0.0);
5139  tmp.set_k(dot.get_k());
5140  accumulate(tmp,scvector(v1),srvector(v2));
5141  dot += tmp;
5142 }
5143 
5145 
5148 inline void accumulate(cidotprecision& dot, const srmatrix_subv& v1, const scmatrix_subv& v2) {
5149  cdotprecision tmp(0.0);
5150  tmp.set_k(dot.get_k());
5151  accumulate(tmp,srvector(v1),scvector(v2));
5152  dot += tmp;
5153 }
5154 
5156 
5159 inline void accumulate(cidotprecision& dot, const scmatrix_subv& v1, const scvector& v2) {
5160  cdotprecision tmp(0.0);
5161  tmp.set_k(dot.get_k());
5162  accumulate(tmp,scvector(v1),v2);
5163  dot += tmp;
5164 }
5165 
5167 
5170 inline void accumulate(cidotprecision& dot, const scmatrix_subv& v1, const srvector& v2) {
5171  cdotprecision tmp(0.0);
5172  tmp.set_k(dot.get_k());
5173  accumulate(tmp,scvector(v1),v2);
5174  dot += tmp;
5175 }
5176 
5178 
5181 inline void accumulate(cidotprecision& dot, const srmatrix_subv& v1, const scvector& v2) {
5182  cdotprecision tmp(0.0);
5183  tmp.set_k(dot.get_k());
5184  accumulate(tmp,srvector(v1),v2);
5185  dot += tmp;
5186 }
5187 
5189 
5192 inline void accumulate(cidotprecision& dot, const scmatrix_subv& v1, const scvector_slice& v2) {
5193  cdotprecision tmp(0.0);
5194  tmp.set_k(dot.get_k());
5195  accumulate(tmp,scvector(v1),v2);
5196  dot += tmp;
5197 }
5198 
5200 
5203 inline void accumulate(cidotprecision& dot, const scmatrix_subv& v1, const srvector_slice& v2) {
5204  cdotprecision tmp(0.0);
5205  tmp.set_k(dot.get_k());
5206  accumulate(tmp,scvector(v1),v2);
5207  dot += tmp;
5208 }
5209 
5211 
5214 inline void accumulate(cidotprecision& dot, const srmatrix_subv& v1, const scvector_slice& v2) {
5215  cdotprecision tmp(0.0);
5216  tmp.set_k(dot.get_k());
5217  accumulate(tmp,srvector(v1),v2);
5218  dot += tmp;
5219 }
5220 
5222 
5225 inline void accumulate(cidotprecision& dot, const scmatrix_subv& v1, const cvector& v2) {
5226  cdotprecision tmp(0.0);
5227  tmp.set_k(dot.get_k());
5228  accumulate(tmp,scvector(v1),v2);
5229  dot += tmp;
5230 }
5231 
5233 
5236 inline void accumulate(cidotprecision& dot, const scmatrix_subv& v1, const rvector& v2) {
5237  cdotprecision tmp(0.0);
5238  tmp.set_k(dot.get_k());
5239  accumulate(tmp,scvector(v1),v2);
5240  dot += tmp;
5241 }
5242 
5244 
5247 inline void accumulate(cidotprecision& dot, const srmatrix_subv& v1, const cvector& v2) {
5248  cdotprecision tmp(0.0);
5249  tmp.set_k(dot.get_k());
5250  accumulate(tmp,srvector(v1),v2);
5251  dot += tmp;
5252 }
5253 
5255 
5258 inline void accumulate(cidotprecision& dot, const scmatrix_subv& v1, const cvector_slice& v2) {
5259  cdotprecision tmp(0.0);
5260  tmp.set_k(dot.get_k());
5261  accumulate(tmp,scvector(v1),v2);
5262  dot += tmp;
5263 }
5264 
5266 
5269 inline void accumulate(cidotprecision& dot, const scmatrix_subv& v1, const rvector_slice& v2) {
5270  cdotprecision tmp(0.0);
5271  tmp.set_k(dot.get_k());
5272  accumulate(tmp,scvector(v1),v2);
5273  dot += tmp;
5274 }
5275 
5277 
5280 inline void accumulate(cidotprecision& dot, const srmatrix_subv& v1, const cvector_slice& v2) {
5281  cdotprecision tmp(0.0);
5282  tmp.set_k(dot.get_k());
5283  accumulate(tmp,srvector(v1),v2);
5284  dot += tmp;
5285 }
5286 
5288 
5291 inline void accumulate(cidotprecision& dot, const scvector& v1, const scmatrix_subv& v2) {
5292  cdotprecision tmp(0.0);
5293  tmp.set_k(dot.get_k());
5294  accumulate(tmp,v1,scvector(v2));
5295  dot += tmp;
5296 }
5297 
5299 
5302 inline void accumulate(cidotprecision& dot, const scvector& v1, const srmatrix_subv& v2) {
5303  cdotprecision tmp(0.0);
5304  tmp.set_k(dot.get_k());
5305  accumulate(tmp,v1,srvector(v2));
5306  dot += tmp;
5307 }
5308 
5310 
5313 inline void accumulate(cidotprecision& dot, const srvector& v1, const scmatrix_subv& v2) {
5314  cdotprecision tmp(0.0);
5315  tmp.set_k(dot.get_k());
5316  accumulate(tmp,v1,scvector(v2));
5317  dot += tmp;
5318 }
5319 
5321 
5324 inline void accumulate(cidotprecision& dot, const scvector_slice& v1, const scmatrix_subv& v2) {
5325  cdotprecision tmp(0.0);
5326  tmp.set_k(dot.get_k());
5327  accumulate(tmp,v1,scvector(v2));
5328  dot += tmp;
5329 }
5330 
5332 
5335 inline void accumulate(cidotprecision& dot, const scvector_slice& v1, const srmatrix_subv& v2) {
5336  cdotprecision tmp(0.0);
5337  tmp.set_k(dot.get_k());
5338  accumulate(tmp,v1,srvector(v2));
5339  dot += tmp;
5340 }
5341 
5343 
5346 inline void accumulate(cidotprecision& dot, const srvector_slice& v1, const scmatrix_subv& v2) {
5347  cdotprecision tmp(0.0);
5348  tmp.set_k(dot.get_k());
5349  accumulate(tmp,v1,scvector(v2));
5350  dot += tmp;
5351 }
5352 
5354 
5357 inline void accumulate(cidotprecision& dot, const cvector& v1, const scmatrix_subv& v2) {
5358  cdotprecision tmp(0.0);
5359  tmp.set_k(dot.get_k());
5360  accumulate(tmp,v1,scvector(v2));
5361  dot += tmp;
5362 }
5363 
5365 
5368 inline void accumulate(cidotprecision& dot, const cvector& v1, const srmatrix_subv& v2) {
5369  cdotprecision tmp(0.0);
5370  tmp.set_k(dot.get_k());
5371  accumulate(tmp,v1,srvector(v2));
5372  dot += tmp;
5373 }
5374 
5376 
5379 inline void accumulate(cidotprecision& dot, const rvector& v1, const scmatrix_subv& v2) {
5380  cdotprecision tmp(0.0);
5381  tmp.set_k(dot.get_k());
5382  accumulate(tmp,v1,scvector(v2));
5383  dot += tmp;
5384 }
5385 
5387 
5390 inline void accumulate(cidotprecision& dot, const cvector_slice& v1, const scmatrix_subv& v2) {
5391  cdotprecision tmp(0.0);
5392  tmp.set_k(dot.get_k());
5393  accumulate(tmp,v1,scvector(v2));
5394  dot += tmp;
5395 }
5396 
5398 
5401 inline void accumulate(cidotprecision& dot, const cvector_slice& v1, const srmatrix_subv& v2) {
5402  cdotprecision tmp(0.0);
5403  tmp.set_k(dot.get_k());
5404  accumulate(tmp,v1,srvector(v2));
5405  dot += tmp;
5406 }
5407 
5409 
5412 inline void accumulate(cidotprecision& dot, const rvector_slice& v1, const scmatrix_subv& v2) {
5413  cdotprecision tmp(0.0);
5414  tmp.set_k(dot.get_k());
5415  accumulate(tmp,v1,scvector(v2));
5416  dot += tmp;
5417 }
5418 
5419 } //namespace cxsc;
5420 
5421 #include "sparsematrix.inl"
5422 
5423 #endif
5424 
scmatrix operator()(const intmatrix &P)
Performs a row permutation using the permutation matrix P. Faster than explicitly computing the produ...
Definition: scmatrix.hpp:673
scmatrix_slice & operator-=(const srmatrix &M)
Assigns the element wise difference of the sparse slice and M to the slice.
Definition: scmatrix.hpp:2204
scmatrix_slice & operator=(const srmatrix &C)
Assing C to the slice.
Definition: scmatrix.hpp:2030
scmatrix_subv & operator-=(const srvector &)
Assign the difference of the subvector with a vector to the subvector.
Definition: scmatrix.hpp:4261
std::vector< complex > & values()
Returns a reference to the vector containing the stored values (the array)
Definition: scmatrix.hpp:92
scmatrix_slice & operator+=(const srmatrix &M)
Assigns the element wise sum of the sparse slice and M to the slice.
Definition: scmatrix.hpp:2156
const complex operator()(const int i, const int j) const
Returns a copy of the element (i,j) of the matrix.
Definition: scmatrix.hpp:2244
scmatrix & operator+=(const rmatrix &B)
Add B to the sparse matrix and assign the result to it.
Definition: scmatrix.hpp:689
scmatrix(const cmatrix &A)
Creates a sparse matrix out of a dense matrix A. Only the non zero elements of A are stored explicitl...
Definition: scmatrix.hpp:390
The Data Type rmatrix_slice.
Definition: rmatrix.hpp:1442
scmatrix & operator=(const real &A)
Assigns a real value to all elements of the matrix (resulting in a dense matrix!)
Definition: scmatrix.hpp:497
scmatrix_subv & operator=(const scvector_slice &v)
Assigns a vector to a subvector.
Definition: scmatrix.hpp:3428
cmatrix_subv & operator=(const scmatrix_subv &rv)
Implementation of standard assigning operator.
Definition: scmatrix.hpp:4386
cmatrix_subv & operator+=(const scmatrix_subv &rv)
Implementation of standard assigning operator.
Definition: scmatrix.hpp:4326
scmatrix_subv & operator=(const real &v)
Assigns v to all elements of the subvector.
Definition: scmatrix.hpp:3403
friend srmatrix Im(const scmatrix_slice &)
Returns the imaginary part of the slice.
Definition: scmatrix.hpp:2379
scmatrix & operator+=(const srmatrix &B)
Add B to the sparse matrix and assign the result to it.
Definition: scmatrix.hpp:709
scmatrix_slice & operator-=(const srmatrix_slice &M)
Assigns the element wise difference of the sparse slice and M to the slice.
Definition: scmatrix.hpp:2192
Helper class for slices of sparse vectors.
Definition: scvector.hpp:1245
cmatrix & operator=(const complex &r)
Implementation of standard assigning operator.
Definition: cmatrix.inl:438
scmatrix & operator=(const cmatrix &A)
Assigns a dense matrix to the sparse matrix. Only the non zero entries of the dense matrix slice are ...
Definition: scmatrix.hpp:512
scmatrix & 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: scmatrix.hpp:517
scmatrix_subv & operator+=(const srvector &)
Assign the sum of the subvector with a vector to the subvector.
Definition: scmatrix.hpp:4241
friend scmatrix Sup(const scimatrix &)
Returns the Supremum of the matrix A.
Definition: scimatrix.hpp:1437
scmatrix_slice & operator/=(const complex &r)
Assigns the component wise division of the sparse slice and M to the slice.
Definition: scmatrix.hpp:2138
cmatrix()
Constructor of class cmatrix.
Definition: cmatrix.inl:31
The Data Type intmatrix.
Definition: intmatrix.hpp:313
friend srmatrix Re(const scmatrix &)
Returns the real part of the sparse matrix A.
Definition: scmatrix.hpp:1004
friend scvector operator-(const scmatrix_subv &)
Unary negation operator.
Definition: scmatrix.hpp:3667
scmatrix_slice & operator-=(const rmatrix_slice &M)
Assigns the element wise difference of the sparse slice and M to the slice.
Definition: scmatrix.hpp:2228
scmatrix_subv operator[](const int)
Returns a row of the matrix.
Definition: scmatrix.hpp:3605
friend int Ub(const scmatrix_subv &)
Returns the upper index bound of the subvector.
Definition: scmatrix.hpp:3535
scmatrix_subv & operator=(const rvector &v)
Assigns a vector to a subvector.
Definition: scmatrix.hpp:3433
friend std::istream & operator>>(std::istream &, scmatrix_slice &)
Standard input operator for sparse matrix slice.
Definition: scmatrix.hpp:3332
int Lb(const cimatrix &rm, const int &i)
Returns the lower bound index.
Definition: cimatrix.inl:1156
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition: cdot.cpp:29
cmatrix_slice & operator *=(const srmatrix &m)
Implementation of addition and assignment operator.
Definition: scmatrix.hpp:1768
complex & element(const int i, const int j)
Returns a reference to the element (i,j) of the matrix.
Definition: scmatrix.hpp:2258
scmatrix_slice & operator=(const srmatrix_slice &C)
Assing C to the slice.
Definition: scmatrix.hpp:2060
scmatrix_slice & operator-=(const scmatrix &M)
Assigns the element wise difference of the sparse slice and M to the slice.
Definition: scmatrix.hpp:2210
scmatrix & operator+=(const scmatrix &B)
Add B to the sparse matrix and assign the result to it.
Definition: scmatrix.hpp:714
void set_k(unsigned int i)
Set precision for computation of dot products.
Definition: cdot.hpp:93
scmatrix_slice & operator=(const scmatrix_slice &C)
Assing C to the slice.
Definition: scmatrix.hpp:2066
scmatrix & operator=(const rmatrix &A)
Assigns a dense matrix to the sparse matrix. Only the non zero entries of the dense matrix slice are ...
Definition: scmatrix.hpp:507
cimatrix & SetLb(cimatrix &m, const int &i, const int &j)
Sets the lower bound index.
Definition: cimatrix.inl:1184
int VecLen(const scimatrix_subv &S)
Returns the length of the subvector.
Definition: scimatrix.hpp:9966
friend srmatrix abs(const scmatrix &)
Returns the componentwise absolute value of the sparse matrix A.
Definition: scmatrix.hpp:1040
friend int Lb(const scmatrix_slice &, const int)
Returns the lower index bound of the rows (if i==ROW) or columns (if i==COL) of the slice.
Definition: scmatrix.hpp:2364
const complex operator()(int i, int j) const
Returns a copy of the element in row i and column j.
Definition: scmatrix.hpp:558
scmatrix & operator-=(const cmatrix_slice &B)
Subtract B from the sparse matrix and assign the result to it.
Definition: scmatrix.hpp:734
scmatrix(const rmatrix &A)
Creates a sparse matrix out of a dense matrix A. Only the non zero elements of A are stored explicitl...
Definition: scmatrix.hpp:367
scmatrix_slice & operator/=(const real &r)
Assigns the component wise division of the sparse slice and M to the slice.
Definition: scmatrix.hpp:2132
scmatrix_slice & operator-=(const cmatrix_slice &M)
Assigns the element wise difference of the sparse slice and M to the slice.
Definition: scmatrix.hpp:2234
friend void SetLb(scmatrix &, const int, const int)
Sets the lower index bound of the rows (i==ROW) or columns (i==COL) to j.
Definition: scmatrix.hpp:926
friend srmatrix Re(const scmatrix_slice &)
Return the real part of the slice.
Definition: scmatrix.hpp:2374
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
scmatrix_slice & operator+=(const rmatrix_slice &M)
Assigns the element wise sum of the sparse slice and M to the slice.
Definition: scmatrix.hpp:2180
A sparse real matrix.
Definition: srmatrix.hpp:77
const std::vector< int > & column_pointers() const
Returns a constant reference to the vector containing the column pointers (the array)
Definition: scmatrix.hpp:97
friend int Lb(const scmatrix &, int)
Returns the lower index bound for the rows or columns of A.
Definition: scmatrix.hpp:956
scmatrix(const int ms, const int ns, const cmatrix &A)
Constructor for banded matrices.
Definition: scmatrix.hpp:416
scmatrix & operator=(const complex &A)
Assigns a complex value to all elements of the matrix (resulting in a dense matrix!...
Definition: scmatrix.hpp:502
scmatrix & operator=(const cmatrix_slice &A)
Assigns a dense matrix slice to the sparse matrix. Only the non zero entries of the dense matrix slic...
Definition: scmatrix.hpp:522
const complex operator[](const int i) const
Returns a copy of the i-th element of the subvector.
Definition: scmatrix.hpp:3386
scmatrix_subv & operator=(const cvector_slice &v)
Assigns a vector to a subvector.
Definition: scmatrix.hpp:3448
A sparse complex interval vector.
Definition: scivector.hpp:62
scmatrix & operator+=(const cmatrix &B)
Add B to the sparse matrix and assign the result to it.
Definition: scmatrix.hpp:694
complex & element(int i, int j)
Returns a reference to the element (i,j) of the matrix.
Definition: scmatrix.hpp:579
scmatrix & operator-=(const scmatrix &B)
Subtract B from the sparse matrix and assign the result to it.
Definition: scmatrix.hpp:744
The Data Type cidotprecision.
Definition: cidot.hpp:57
friend scmatrix Inf(const scimatrix &)
Returns the Infimum of the matrix A.
Definition: scimatrix.hpp:1419
friend int RowLen(const scmatrix_slice &)
Returns the number columns of the matrix slice.
Definition: scmatrix.hpp:2328
A sparse real vector.
Definition: srvector.hpp:58
friend int RowLen(const scmatrix &)
Returns the number of columns of the matrix.
Definition: scmatrix.hpp:979
cimatrix Id(const cimatrix &A)
Returns the Identity matrix.
Definition: cimatrix.cpp:61
The Data Type rvector_slice.
Definition: rvector.hpp:1063
A sparse complex vector.
Definition: scvector.hpp:58
The Data Type cvector.
Definition: cvector.hpp:57
scmatrix & operator+=(const rmatrix_slice &B)
Add B to the sparse matrix and assign the result to it.
Definition: scmatrix.hpp:699
friend int Ub(const scmatrix_slice &, const int)
Returns the upper index bound of the rows (if i==ROW) or columns (if i==COL) of the slice.
Definition: scmatrix.hpp:2369
scmatrix_slice & operator-=(const scmatrix_slice &M)
Assigns the element wise difference of the sparse slice and M to the slice.
Definition: scmatrix.hpp:2198
cimatrix & SetUb(cimatrix &m, const int &i, const int &j)
Sets the upper bound index.
Definition: cimatrix.inl:1191
friend scmatrix transp(const scmatrix &)
Returns the transpose of A.
Definition: scmatrix.hpp:886
scmatrix & operator+=(const cmatrix_slice &B)
Add B to the sparse matrix and assign the result to it.
Definition: scmatrix.hpp:704
friend int ColLen(const scmatrix &)
Returns the number of rows of the matrix.
Definition: scmatrix.hpp:984
The Data Type rvector.
Definition: rvector.hpp:57
The Data Type cmatrix.
Definition: cmatrix.hpp:513
cmatrix_slice & operator+=(const srmatrix &m)
Implementation of addition and assignment operator.
Definition: scmatrix.hpp:1736
civector operator *(const cimatrix_subv &rv, const cinterval &s)
Implementation of multiplication operation.
Definition: cimatrix.inl:731
scmatrix_subv & operator=(const complex &v)
Assigns v to all elements of the subvector.
Definition: scmatrix.hpp:3408
int get_nnz() const
Returns the number of non-zero entries (including explicitly stored zeros).
Definition: srmatrix.hpp:630
friend int VecLen(const scmatrix_subv &)
Returns the length of the subvector.
Definition: scmatrix.hpp:3543
A slice of a sparse real matrix.
Definition: srmatrix.hpp:1360
The Data Type cdotprecision.
Definition: cdot.hpp:60
scmatrix_slice & operator=(const real &C)
Assing C to all elements of the slice.
Definition: scmatrix.hpp:2020
cmatrix_slice & operator=(const cmatrix &m)
Implementation of standard assigning operator.
Definition: cmatrix.inl:450
A sparse complex matrix.
Definition: scmatrix.hpp:69
scmatrix & operator-=(const rmatrix &B)
Subtract B from the sparse matrix and assign the result to it.
Definition: scmatrix.hpp:719
A slice of a sparse complex matrix.
Definition: scmatrix.hpp:1956
friend scmatrix Id(const scmatrix &)
Return a sparse unity matrix of the same dimension as A.
Definition: scmatrix.hpp:863
void Resize(cimatrix &A)
Resizes the matrix.
Definition: cimatrix.inl:1211
std::vector< int > & row_indices()
Returns a reference to the vector containing the row indices (the array)
Definition: scmatrix.hpp:87
scmatrix(const int m, const int n, const int nnz, const int *rows, const int *cols, const complex *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: scmatrix.hpp:254
scmatrix_subv & operator=(const scmatrix_subv &v)
Assigns a vector to a subvector.
Definition: scmatrix.hpp:3458
scmatrix_subv & operator=(const srvector &v)
Assigns a vector to a subvector.
Definition: scmatrix.hpp:3413
friend srvector Re(const scmatrix_subv &)
Returns the real part of the subvector.
Definition: scmatrix.hpp:3548
scmatrix_slice & operator=(const rmatrix_slice &C)
Assing C to the slice.
Definition: scmatrix.hpp:2050
scmatrix & operator/=(const complex &r)
Divide all elements of the sparse matrix by r and assign the result to it.
Definition: scmatrix.hpp:794
scmatrix_subv operator[](const cxscmatrix_column &)
Returns a column of the matrix as a sparse subvector object.
Definition: scmatrix.hpp:3573
Represents a row or column vector of a sparse matrix.
Definition: srmatrix.hpp:2157
The Data Type rmatrix.
Definition: rmatrix.hpp:470
scvector()
Default constructor, creates an empty vector of size 0.
Definition: scvector.hpp:68
scmatrix_subv & operator=(const cvector &v)
Assigns a vector to a subvector.
Definition: scmatrix.hpp:3438
int ColLen(const cimatrix &)
Returns the column dimension.
Definition: cimatrix.inl:1202
friend scmatrix mid(const scimatrix &)
Returns the componentwise midpoint of the matrix A.
Definition: scimatrix.hpp:1491
scmatrix & operator=(const srmatrix &A)
Assigns a sparse real matrix to the sparse complex matrix.
Definition: scmatrix.hpp:527
The Scalar Type complex.
Definition: complex.hpp:49
scmatrix(const srmatrix &A)
Creates a sparse complex matrix out of a sparse real matrix.
Definition: scmatrix.hpp:359
scmatrix_subv & operator=(const scvector &v)
Assigns a vector to a subvector.
Definition: scmatrix.hpp:3418
scmatrix & operator *=(const cmatrix &B)
Multiply the sparse matrix by B and assign the result to it.
Definition: scmatrix.hpp:749
The Data Type cmatrix_subv.
Definition: cmatrix.hpp:53
const std::vector< complex > & values() const
Returns a constant reference to the vector containing the stored values (the array)
Definition: scmatrix.hpp:107
scmatrix_slice & operator+=(const scmatrix_slice &M)
Assigns the element wise sum of the sparse slice and M to the slice.
Definition: scmatrix.hpp:2150
int get_k() const
Get currently set precision for computation of dot products.
Definition: cidot.hpp:89
int RowLen(const cimatrix &)
Returns the row dimension.
Definition: cimatrix.inl:1199
friend std::ostream & operator<<(std::ostream &, const scmatrix_slice &)
Standard output operator for sparse matrix slice.
Definition: scmatrix.hpp:3322
cmatrix & operator+=(const srmatrix &m)
Implementation of addition and assignment operator.
Definition: scmatrix.hpp:1728
rmatrix CompMat(const cimatrix &A)
Returns Ostrowski's comparison matrix.
Definition: cimatrix.cpp:45
The Data Type cmatrix_slice.
Definition: cmatrix.hpp:1202
void dropzeros()
Drops explicitly stored zeros from the data structure.
Definition: scmatrix.hpp:473
scmatrix_subv & operator=(const srmatrix_subv &v)
Assigns a vector to a subvector.
Definition: scmatrix.hpp:3453
real density() const
Returns the density (the number of non-zeros divided by the number of elements) of the matrix.
Definition: scmatrix.hpp:679
friend scmatrix diam(const scimatrix &)
Returns the componentwise diameter of the matrix A.
Definition: scimatrix.hpp:1509
int get_nnz() const
Returns the number of non-zero entries (including explicitly stored zeros).
Definition: scmatrix.hpp:684
scmatrix_slice & operator+=(const cmatrix_slice &M)
Assigns the element wise sum of the sparse slice and M to the slice.
Definition: scmatrix.hpp:2186
friend srmatrix CompMat(const scmatrix &)
Returns Ostrowskis comparison matrix for A.
Definition: scmatrix.hpp:1050
scmatrix_slice & operator+=(const cmatrix &M)
Assigns the element wise sum of the sparse slice and M to the slice.
Definition: scmatrix.hpp:2174
scmatrix_slice & operator-=(const rmatrix &M)
Assigns the element wise difference of the sparse slice and M to the slice.
Definition: scmatrix.hpp:2216
scmatrix_slice & operator=(const scmatrix &C)
Assing C to the slice.
Definition: scmatrix.hpp:2035
scmatrix(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: scmatrix.hpp:119
friend int Lb(const scmatrix_subv &)
Returns the lower index bound of the subvector.
Definition: scmatrix.hpp:3527
scmatrix operator()(const intmatrix &P, const intmatrix &Q)
Performs row and column permutations using the two permutation matrices P and Q. Faster than explicit...
Definition: scmatrix.hpp:666
cmatrix & operator *=(const srmatrix &m)
Implementation of addition and assignment operator.
Definition: scmatrix.hpp:1760
scmatrix_slice & operator+=(const srmatrix_slice &M)
Assigns the element wise sum of the sparse slice and M to the slice.
Definition: scmatrix.hpp:2144
Represents a row or column vector of a sparse matrix.
Definition: scimatrix.hpp:9628
STORAGE_TYPE
Enumeration depicting the storage type of a sparse matrix (Triplet storage, Compressed column storage...
Definition: srmatrix.hpp:42
Helper class for slices of sparse vectors.
Definition: srvector.hpp:868
scmatrix_slice & operator-=(const cmatrix &M)
Assigns the element wise difference of the sparse slice and M to the slice.
Definition: scmatrix.hpp:2222
scmatrix & operator-=(const cmatrix &B)
Subtract B from the sparse matrix and assign the result to it.
Definition: scmatrix.hpp:724
scmatrix()
Standard constructor, creates an empty matrix of dimension 0x0.
Definition: scmatrix.hpp:112
scmatrix_subv & operator=(const srvector_slice &v)
Assigns a vector to a subvector.
Definition: scmatrix.hpp:3423
scmatrix(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: scmatrix.hpp:128
int Ub(const cimatrix &rm, const int &i)
Returns the upper bound index.
Definition: cimatrix.inl:1163
complex & operator[](const int i)
Returns a reference to the i-th element of the subvector.
Definition: scmatrix.hpp:3366
scmatrix_subv & operator/=(const real &)
Assign the componentwise division of the subvector with a scalar to the subvector.
Definition: scmatrix.hpp:4231
void full(cmatrix &A) const
Creates a full matrix out of the sparse matrix and stores it in A. This should normally be done using...
Definition: scmatrix.hpp:458
Represents a row or column vector of a sparse matrix.
Definition: scmatrix.hpp:3345
scmatrix_subv & operator=(const rvector_slice &v)
Assigns a vector to a subvector.
Definition: scmatrix.hpp:3443
civector operator/(const cimatrix_subv &rv, const cinterval &s)
Implementation of division operation.
Definition: cimatrix.inl:730
cimatrix transp(const cimatrix &A)
Returns the transposed matrix.
Definition: cimatrix.cpp:74
friend void SetUb(scmatrix &, const int, const int)
Sets the upper index bound of the rows (i==ROW) or columns (i==COL) to j.
Definition: scmatrix.hpp:941
scmatrix & operator-=(const rmatrix_slice &B)
Subtract B from the sparse matrix and assign the result to it.
Definition: scmatrix.hpp:729
The Data Type cvector_slice.
Definition: cvector.hpp:844
friend srvector Im(const scmatrix_subv &)
Returns the imaginary part of the subvector.
Definition: scmatrix.hpp:3553
scmatrix_slice & operator+=(const rmatrix &M)
Assigns the element wise sum of the sparse slice and M to the slice.
Definition: scmatrix.hpp:2168
scmatrix(const int m, const int n, const int nnz, const intvector &rows, const intvector &cols, const cvector &values, const enum STORAGE_TYPE t=triplet)
Creates a sparse matrix out of three vectors (arrays) forming a matrix stored in triplet,...
Definition: scmatrix.hpp:143
scmatrix_slice & operator=(const cmatrix &C)
Assing C to the slice.
Definition: scmatrix.hpp:2045
std::vector< int > & column_pointers()
Returns a reference to the vector containing the column pointers (the array)
Definition: scmatrix.hpp:82
A slice of a sparse complex interval matrix.
Definition: scimatrix.hpp:4918
scmatrix_slice & operator=(const complex &C)
Assing C to all elements of the slice.
Definition: scmatrix.hpp:2025
scmatrix_slice & operator *=(const srmatrix_slice &M)
Assigns the product of the sparse slice and M to the slice.
Definition: scmatrix.hpp:2072
scmatrix_subv & operator *=(const real &)
Assign the componentwise product of the subvector with a scalar to the subvector.
Definition: scmatrix.hpp:4221
friend int Ub(const scmatrix &, int)
Returns the upper index bound for the rows or columns of A.
Definition: scmatrix.hpp:969
scmatrix_slice & operator+=(const scmatrix &M)
Assigns the element wise sum of the sparse slice and M to the slice.
Definition: scmatrix.hpp:2162
scmatrix_slice & operator=(const rmatrix &C)
Assing C to the slice.
Definition: scmatrix.hpp:2040
scmatrix operator()(const intvector &pervec, const intvector &q)
Performs a row and column permutation using two permutation vectors.
Definition: scmatrix.hpp:615
The Scalar Type real.
Definition: real.hpp:113
friend std::istream & operator>>(std::istream &, scmatrix_subv &)
Standard input operator for subvectors.
Definition: scmatrix.hpp:3564
The Data Type intvector.
Definition: intvector.hpp:51
const std::vector< int > & row_indices() const
Returns a constant reference to the vector containing the row indices (the array)
Definition: scmatrix.hpp:102
scmatrix_slice & operator=(const cmatrix_slice &C)
Assing C to the slice.
Definition: scmatrix.hpp:2055
friend int ColLen(const scmatrix_slice &)
Returns the number of rows of the matrix slice.
Definition: scmatrix.hpp:2333
friend srmatrix Im(const scmatrix &)
Returns the imaginary part of the sparse matrix A.
Definition: scmatrix.hpp:1022
The Data Type cimatrix.
Definition: cimatrix.hpp:907
ivector abs(const cimatrix_subv &mv)
Returns the absolute value of the matrix.
Definition: cimatrix.inl:737
scmatrix & operator/=(const real &r)
Divide all elements of the sparse matrix by r and assign the result to it.
Definition: scmatrix.hpp:789
scmatrix & operator-=(const srmatrix &B)
Subtract B from the sparse matrix and assign the result to it.
Definition: scmatrix.hpp:739
scmatrix operator()(const intvector &pervec)
Performs a row permutation using a permutation vector.
Definition: scmatrix.hpp:642