TBCI Numerical high perf. C++ Library 2.8.0
cscmatrix.h
Go to the documentation of this file.
1
5//----------------------------------------------------------
6// version 1.0, 08.1998
7//
8// Andreas Ahland
9//
10// $Id: cscmatrix.h,v 1.18.2.40 2022/11/03 17:28:10 garloff Exp $
11//----------------------------------------------------------
12
13#ifndef TBCI_CSCMATRIX_H
14#define TBCI_CSCMATRIX_H
15
16#include "tbci/vector.h"
17#include "tbci/matrix_sig.h"
18#include "tbci/f_matrix.h"
19#include "tbci/band_matrix.h"
20
21// Avoid -fguiding-decls
22#if !defined(NO_GD) && !defined(AUTO_DECL)
23# include "tbci/cscmatrix_gd.h"
24#endif
25
27
28// exception class
29#ifndef TBCI_DISABLE_EXCEPT
30//# include "except.h" is include in basics.h
31
32class CSCMatErr : public NumErr
33{
34 public:
36 : NumErr("Error in CSCMatrix") {}
37 CSCMatErr(const char* t, const long i = 0)
38 : NumErr(t, i) {}
40 : NumErr(ce) {}
41 virtual ~CSCMatErr() throw() {}
42};
43#endif
44
45#ifdef PRAGMA_I
46# pragma interface "cscmatrix.h"
47#endif
48
49template <typename T> class CSCMatrix;
50template <typename T> class F_TSMatrix;
51template <typename T> class F_TMatrix;
52template <typename T> class F_Matrix;
53
54//----------
55// CSCMatrix
56//----------
64template <typename T>
65class CSCMatrix : public Matrix_Sig<T>
66{
67 friend class F_TSMatrix<T>;
68
69public:
70 typedef T value_type;
71 typedef T element_type;
72 typedef T aligned_value_type TALIGN(MIN_ALIGN2);
73
75 CSCMatrix() { allocate(1,1); }
76 CSCMatrix(const T& val, const unsigned int rows,
77 const unsigned int columns, const unsigned int nnzeros=1)
78 { allocate(rows, columns, nnzeros);
79 if (val != (T)0) fill(val); }
80 CSCMatrix(const unsigned int rows, const unsigned int columns,
81 const unsigned int nnzeros=1)
82 { allocate(rows, columns, nnzeros); }
84 { copy(m); }
85 CSCMatrix(const F_Matrix<T>& m, const double tol=0)
86 { allocate(m.rows(), m.columns()); fill(m,tol); }
87 CSCMatrix(const BdMatrix<T>& m, const double tol=0)
88 { allocate(m.rows(), m.columns()); fill(m,tol); }
89 CSCMatrix(const Matrix<T>& m, const double tol=0)
90 { allocate(m.rows(), m.columns()); fill(m,tol); }
92 { destroy (); }
93
94 operator F_TMatrix<T> () const;
95
97 /*virtual*/ static const char* mat_info () { return ("CSCMatrix"); }
98
99
101 unsigned int rows() const { return n_rows; }
102 unsigned int columns() const { return n_cols; }
103 unsigned int size() const { return n_size; } // number of non-zero elements
104
105
107 CSCMatrix<T>& resize(const unsigned int newRows,
108 const unsigned int newColumns,
109 const unsigned int nnzeros = 1);
110
114 CSCMatrix<T>& fill (const T&);
116 CSCMatrix<T>& setunit (const T& = (T)1);
117
119 typename tbci_traits<T>::const_refval_type
120 operator() (unsigned int row, unsigned int column) const HOT;
121 typename tbci_traits<T>::const_refval_type
122 get (unsigned int row, unsigned int column) const HOT;
123
125 CSCMatrix<T>& setval(const T& z, unsigned int row, unsigned int column);
126 T& setval(unsigned int row, unsigned int column);
127 T& operator () (unsigned int row, unsigned int column)
128 { return setval (row, column); }
129 T* dataPointer() { return comp; }
130 unsigned int* columnPointer() { return pcol; }
131 unsigned int* rowIndexPointer() { return irow; }
132
134 bool operator == (const CSCMatrix<T>& m) const;
135 bool operator != (const CSCMatrix<T>& m) const
136 { return !(this->operator == (m)); }
138
140 //CSCMatrix<T>& do_import (const Matrix_Sig<T>& M)
141 template <typename MatType>
142 CSCMatrix<T>& do_import /*FGD*/ (const MatType& M)
143 {
144 destroy (); allocate(M.rows(), M.columns(), 20);
145 fill (M);
146 return *this;
147 }
148
150 //void do_export (Matrix_Sig<T>& M)
151 template <typename MatType>
152 void do_export /*FGD*/ (MatType& M)
153 {
154 M.clear (); //Matrix should be cleared!
155 for (unsigned int i=0; i<n_cols; ++i)
156 //for (unsigned int j=pcol[i]; j<pcol[i+1]; ++j) M.setval(comp[j], i,irow[j]);
157 for (unsigned int j=pcol[i]; j<pcol[i+1]; ++j)
158 M(i,irow[j]) = comp[j];
159 }
160
163
167
170
172 //friend F_TMatrix<T> operator* FGD (const F_Matrix<T>& m1, const CSCMatrix<T>& m2);
173 F_TMatrix<T> multf (const F_Matrix<T>&) const;
174 F_TMatrix<T> mult (const F_Matrix<T>&) const;
175 F_TMatrix<T> mult1 (const F_Matrix<T>&) const;
176
178 TVector<T> operator * (const Vector<T>& v) const HOT;
179 TVector<T> operator * (const TVector<T>& tv) const HOT;
181 TVector<T> cscm_vec_mul_exact (const Vector<T>& V) const;
183 void MatVecMult (Vector<T>& res, const Vector<T>& v) const HOT;
184 void MatVecMult (T* v, T* res) HOT; // res <- A*v , for ARPACK++
185
187 TVector<T> transMult(const Vector<T>& v) const HOT;
188 TVector<T> transMult(const TVector<T>& tv) const HOT;
189 TVector<T> transMult(const TSVector<T>& tsv) const HOT;
190
191
192 CSCMatrix<T> operator* (const T& z) const; // matrix-scalar multiplication
194 CSCMatrix<T> mult (const T&) const;
195
196 CSCMatrix<T>& operator*= (const T& z); // fast in-place matrix-scalar multiplication
197 CSCMatrix<T> operator/ (const T& z) const; // matrix-scalar division
198 CSCMatrix<T>& operator/= (const T& z); // fast in-place matrix-scalar division
199
204
206 friend STD__ ostream& operator<< FGDT (STD__ ostream& stream, const CSCMatrix<T>& m);
207
214protected:
216 unsigned int* pcol; // pointer to the beginning of each column
217 unsigned int* irow; // row index of all nonzero elements
218 T* comp; // storage vector
219 mutable T dummy;
220
221 void allocate (unsigned int rows, unsigned int columns, unsigned int nnzeros=1);
222 void destroy (); // used by destructor and "=" operator
223 void copy (const CSCMatrix<T>& m); // used by copy constr. and "="
224 void insert (const unsigned int column, const unsigned int pos); //dynamic inserting
225
227 //void fill (const Matrix_Sig<T>& M, double tol=0)
228 template <typename MatType>
229 void fill /*FGD*/ (const MatType& M, const double tol = 0.0)
230 {
231 unsigned int count = 0, i, j;
232 //count nonzeros
233 for (j = 0; j < n_cols; ++j)
234 for (i = 0; i < n_rows; ++i)
235 if (MATH__ fabs(M(i,j)) > tol)
236 ++count;
237 STD__ cout << "Matrix conversion, NNZ=" << count << STD__ endl;
238 resize(M.rows(), M.columns(), count);
239 for (j = 0; j < n_cols; ++j)
240 for (i = 0; i < n_rows; ++i)
241 if (MATH__ fabs(M(i,j)) > tol)
242 setval(M(i,j), i,j);
243 }
244};
245
246
247template <typename T>
248STD__ ostream& operator<< (STD__ ostream& stream, const CSCMatrix<T>& m)
249{
250 for (unsigned int i = 0; i < m.rows(); ++i) {
251 //stream << "|\t";
252 for (unsigned int j = 0; j < m.columns(); ++j)
253 stream << m(i,j) << " " /* << "\t"*/;
254 stream /*<< "|" */ << "\n";
255 }
256 return stream /* << STD__ flush*/;
257}
258
259//-------------------------------
260// CSCMatrix function definitions
261//-------------------------------
262
264template <typename T>
266{
267 F_TMatrix<T> res((T)0, m1.rows(), columns());
268 BCHK(rows() != m1.columns(), CSCMatErr, Sizes dont match in multf, m1.columns(), res);
269 if (do_exactsum2()) {
270 for (unsigned int i = 0; i < m1.rows(); ++i)
271 for (unsigned int j = 0; j < columns(); ++j) {
272 REGISTER T val = 0;
273 REGISTER T corr = 0;
274 for (unsigned int x = pcol[j]; x < pcol[j+1]; ++x) {
275 const REGISTER T y = m1(i, irow[x] )*comp[x];
276 const REGISTER T t = val + y;
277 corr += (t-val) - y;
278 val = t;
279 }
280 res.setval(val-corr, i, j);
281 }
282 } else {
283 for (unsigned int i = 0; i < m1.rows(); ++i)
284 for (unsigned int j = 0; j < columns(); ++j) {
285 REGISTER T val = 0;
286 for (unsigned int x = pcol[j]; x < pcol[j+1]; ++x)
287 val += m1(i, irow[x] )*comp[x];
288 res.setval(val, i, j);
289 }
290 }
291 return res;
292}
293
294INST(template<typename T> class CSCMatrix<T> friend F_TMatrix<T> operator* (const F_Matrix<T>&, const CSCMatrix<T>&);)
295template <typename T>
296inline F_TMatrix<T> operator * (const F_Matrix<T>& m1, const CSCMatrix<T>& m2)
297{ return m2.multf (m1); }
298
300template <typename T>
302{
303 F_TMatrix<T> res((T)0, rows(), m1.columns());
304 BCHK(columns() != m1.rows(), CSCMatErr, Sizes dont match in mult, columns(), res);
305 if (do_exactsum2()) {
306 for (unsigned int i = 0; i < rows(); ++i)
307 for (unsigned int j = 0; j < m1.columns(); ++j) {
308 REGISTER T val = 0;
309 REGISTER T corr = 0;
310 for (unsigned int x = 0; x < columns(); ++x) {
311 const REGISTER T y = (*this)(i, x) * m1(x, j);
312 const REGISTER T t = val + y;
313 corr += (t-val) - y;
314 val = t;
315 }
316 res.setval(val-corr, i, j);
317 }
318 } else {
319 for (unsigned int i = 0; i < rows(); ++i)
320 for (unsigned int j = 0; j < m1.columns(); ++j) {
321 REGISTER T val = 0;
322 for (unsigned int x = 0; x < columns(); ++x)
323 val += (*this)(i, x) * m1(x, j);
324 res.setval(val, i, j);
325 }
326 }
327 return res;
328}
329
331template <typename T>
333{
334 F_TMatrix<T> res((T)0, rows(), m1.columns());
335 BCHK(columns() != m1.rows(), CSCMatErr, Sizes dont match in mult, columns(), res);
336 for (unsigned int c = 0; c < columns(); ++c)
337 for (unsigned int off = pcol[c]; off < pcol[c+1]; ++off) {
338 /* No exact summing possible */
339 const unsigned int r = irow[off];
340 const T el = comp[off]; // (*this)(r,c)
341 for (unsigned int j = 0; j < m1.columns(); ++j)
342 res.setval(r, j) += el * m1(c, j);
343 }
344 return res;
345}
346
347INST(template<typename T> class CSCMatrix<T> friend F_TMatrix<T> operator* (const CSCMatrix<T>&, const F_Matrix<T>&);)
348template <typename T>
349inline F_TMatrix<T> operator * (const CSCMatrix<T>& m1, const F_Matrix<T>& m2)
350{ return m1.mult1 (m2); }
351
353template <typename T>
355{
356 F_TMatrix<T> res((T)0, rows(), m1.columns());
357 BCHK(columns() != m1.rows(), CSCMatErr, Sizes dont match in mult, columns(), res);
358 if (do_exactsum2()) {
359 for (unsigned int r0 = 0; r0 < rows(); ++r0)
360 for (unsigned int c1 = 0; c1 < m1.columns(); ++c1) {
361 REGISTER T val = 0;
362 REGISTER T corr = 0;
363 for (unsigned int off1 = m1.pcol[c1]; off1 < m1.pcol[c1+1]; ++off1) {
364 const REGISTER T y = (*this)(r0, m1.irow[off1]) * m1.comp[off1];
365 const REGISTER T t = val + y;
366 corr += (t-val) - y;
367 val = t;
368 }
369 res.setval(val-corr, r0, c1);
370 }
371 } else {
372 for (unsigned int r0 = 0; r0 < rows(); ++r0)
373 for (unsigned int c1 = 0; c1 < m1.columns(); ++c1) {
374 REGISTER T val = 0;
375 for (unsigned int off1 = m1.pcol[c1]; off1 < m1.pcol[c1+1]; ++off1)
376 val += (*this)(r0, m1.irow[off1]) * m1.comp[off1];
377 res.setval(val, r0, c1);
378 }
379 }
380 return res;
381}
382
384template <typename T>
386{
388 BCHK(columns() != m1.columns(), CSCMatErr, Columns dont match in add, columns(), res);
389 BCHK(rows() != m1.rows(), CSCMatErr, Rows dont match in add, rows(), res);
390 for (unsigned int c = 0; c < columns(); ++c) {
391 unsigned int off0 = pcol[c];
392 unsigned int off1 = m1.pcol[c];
393 while (off0 < pcol[c+1] || off1 < m1.pcol[c+1]) {
394 // only one is valid ...
395 if (off0 >= pcol[c+1]) {
396 res(m1.irow[off1], c) = m1.comp[off1];
397 ++off1;
398 } else if (off1 >= m1.pcol[c+1]) {
399 res(irow[off0], c) = comp[off0];
400 ++off0;
401 } else if (irow[off0] == m1.irow[off1]) {
402 // both are valid
403 res(irow[off0], c) = comp[off0] + m1.comp[off1];
404 ++off0; ++off1;
405 } else if (irow[off0] < m1.irow[off1]) {
406 res(irow[off0], c) = comp[off0];
407 ++off0;
408 } else {
409 res(m1.irow[off1], c) = m1.comp[off1];
410 ++off1;
411 }
412 }
413 }
414 return res;
415}
416
417template <typename T>
419{
421 BCHK(columns() != m1.columns(), CSCMatErr, Columns dont match in sub, columns(), res);
422 BCHK(rows() != m1.rows(), CSCMatErr, Rows dont match in sub, rows(), res);
423 for (unsigned int c = 0; c < columns(); ++c) {
424 unsigned int off0 = pcol[c];
425 unsigned int off1 = m1.pcol[c];
426 while (off0 < pcol[c+1] || off1 < m1.pcol[c+1]) {
427 // only one is valid ...
428 if (off0 >= pcol[c+1]) {
429 res(m1.irow[off1], c) = -m1.comp[off1];
430 ++off1; continue;
431 }
432 if (off1 >= m1.pcol[c+1]) {
433 res(irow[off0], c) = comp[off0];
434 ++off0; continue;
435 }
436 // both are valid
437 if (irow[off0] == m1.irow[off1]) {
438 res(irow[off0], c) = comp[off0] - m1.comp[off1];
439 ++off0; ++off1; //continue;
440 } else if (irow[off0] < m1.irow[off1]) {
441 res(irow[off0], c) = comp[off0];
442 ++off0; //continue;
443 } else {
444 res(m1.irow[off1], c) = -m1.comp[off1];
445 ++off1; //continue;
446 }
447 }
448 }
449 return res;
450}
451
453template <typename T>
454inline typename tbci_traits<T>::const_refval_type
455 CSCMatrix<T>::get(unsigned int row, unsigned int col) const
456{
457 static T dummy = (T)0;
458/*
459 // Version 1
460 for (unsigned int j=pcol[col]; j<pcol[col+1]; j++)
461 {
462 if (irow[j] == row) return comp[j];
463 }
464
465 // Version 2
466 unsigned int pos = pcol[col];
467 // linear search ...
468 while (pos < pcol[col+1] && irow[pos] < row)
469 ++pos;
470 if (pos < n_size && irow[pos] == row)
471 return comp[pos];
472 // Version 3
473
474 unsigned int pos = pcol[col];
475 // linear search ...
476 while (pos < pcol[col+1]) {
477 if (irow[pos] < row)
478 ++pos;
479 else if (irow[pos] > row)
480 break;
481 else
482 return comp[pos];
483 }
484 // Element not found
485 return (dummy = 0);
486*/
487
488/*
489 // Version 4
490 const unsigned int epos = pcol[col+1];
491 for (unsigned int pos = pcol[col]; pos < epos; ++pos) {
492 if (irow[pos] == row)
493 return comp[pos];
494 else if (irow[pos] > row)
495 break;
496 }
497 return (dummy = 0);
498*/
499 // Version 5
500 const unsigned int spos = pcol[col], epos = pcol[col+1];
501 if (epos == spos)
502 return (dummy = 0);
503 if (irow[spos] > row || irow[epos-1] < row)
504 return (dummy = 0);
505 if (irow[spos] == row)
506 return comp[spos];
507 if (irow[epos-1] == row)
508 return comp[epos-1];
509 /* Binary search */
510 if ((epos - spos) >= 32) {
511 const unsigned long int off = _bin_search(irow, row, spos, epos);
512 if (off != (unsigned long)(-1))
513 return comp[off];
514 else
515 return (dummy = 0);
516 }
517 /* Linear search, be clever at which side to start */
518 if (row <= (irow[spos] + irow[epos-1] + 1) / 2) {
519 for (REGISTER unsigned int pos = spos; pos < epos; ++pos) {
520 if (irow[pos] == row)
521 return comp[pos];
522 else if (irow[pos] > row)
523 break;
524 }
525 } else {
526 for (REGISTER unsigned int pos = epos; pos > spos; --pos) {
527 const unsigned po1 = pos-1;
528 if (irow[po1] == row)
529 return comp[po1];
530 else if (irow[po1] < row)
531 break;
532 }
533 }
534 return (dummy = 0);
535 // TODO: Should we do binary search?
536}
537
538template <typename T>
539inline typename tbci_traits<T>::const_refval_type
540 CSCMatrix<T>::operator() (const unsigned row, const unsigned col) const
541{
542 EXPCHK(row >= n_rows, CSCMatErr, Illegal row index, row, dummy);
543 EXPCHK(col >= n_cols, CSCMatErr, Illegal col index, col, dummy);
544 return get(row, col);
545}
546
548template <typename T>
550{
551 F_TMatrix<T> res((T)0, rows(), columns());
552 for (unsigned int c = 0; c < columns(); ++c)
553 for (unsigned int off = pcol[c]; off < pcol[c+1]; ++off)
554 res.setval(comp[off], irow[off], c);
555 return res;
556}
557
558
560template <typename T>
561CSCMatrix<T>& CSCMatrix<T>::resize(const unsigned int newRows,
562 const unsigned int newColumns,
563 const unsigned int nnzeros)
564{
565 if (UNLIKELY(newRows != n_rows && newColumns != n_cols)) {
566 destroy ();
567 allocate(newRows, newColumns, nnzeros);
568 } else if (UNLIKELY(nnzeros != n_max_size)) {
569 // TODO: Handle size 0
570 unsigned int *newirow = new unsigned int[nnzeros];
571 T *newcomp = new T [nnzeros];
572 unsigned int i;
573 const unsigned int com_sz = MIN(nnzeros, n_max_size);
574 for (i = 0; i < com_sz; ++i) {
575 newirow[i] = irow[i];
576 newcomp[i] = comp[i];
577 }
578 for (; i < nnzeros; ++i)
579 newirow[i] = 0;
580 n_max_size = nnzeros;
581 delete[] irow;
582 delete[] comp;
583 irow = newirow;
584 comp = newcomp;
585 } /* else do nothing */
586 return *this;
587}
588
589
590template <typename T>
591void CSCMatrix<T>::insert(const unsigned int column, const unsigned int pos)
592{
593 if (UNLIKELY(n_size >= n_max_size)) {
594 //not enough space left
595 n_max_size += 16;
596 unsigned int* new_irow = new unsigned int [n_max_size];
597 BCHKNR(new_irow==NULL, CSCMatErr, Out of memory, n_max_size);
598 T* new_comp = new T [n_max_size];
599 BCHKNR(new_comp==NULL, CSCMatErr, Out of memory, n_max_size);
600
601 // if data exist, copy data
602 if (n_size > 0) {
603 unsigned int i;
604 for (i = 0; i < pos; ++i) {
605 // copy data up before
606 new_irow[i] = irow[i];
607 new_comp[i] = comp[i];
608 }
609 for (i = pos; i < n_size; ++i) {
610 // copy data behind insertion
611 new_irow[i+1] = irow[i];
612 new_comp[i+1] = comp[i];
613 }
614 // Following elements are not yet used
615 for (i = n_size+1; i < n_max_size; ++i)
616 new_irow[i] = 0;
617 }
618
619 if (irow != NULL)
620 delete[] irow;
621 if (comp != NULL)
622 delete[] comp;
623 irow = new_irow;
624 comp = new_comp;
625 } else if (n_size > 0) {
626 // if data exist, copy data
627 for (int j = n_size-1; j >= (int)pos; --j) {
628 // copy data behind insertion
629 irow[j+1] = irow[j];
630 comp[j+1] = comp[j];
631 }
632 }
633 for (unsigned int i = column+1; i < n_cols+1; ++i)
634 ++pcol[i];
635 ++n_size;
636}
637
638// element access (write)
639template <typename T>
640CSCMatrix<T>& CSCMatrix<T>::setval(const T& z, unsigned int row, unsigned int column)
641{
642 setval(row, column) = z;
643 return *this;
644}
645
647template <typename T>
648T& CSCMatrix<T>::setval(unsigned int row, unsigned int column)
649{
650 EXPCHKNR(row >= n_rows, CSCMatErr, Illegal row index, row);
651 EXPCHKNR(column >= n_cols, CSCMatErr, Illegal col index, column);
652 unsigned int pos = pcol[column];
653 const unsigned int posn = pcol[column+1];
654 while (pos < posn && irow[pos] < row)
655 ++pos;
656 if (pos < posn && irow[pos] == row)
657 return comp[pos];
658
659 // element not found -- dynamically growing required...
660 insert(column, pos);
661
662 // Write new element
663 irow[pos] = row;
664 comp[pos] = T(0); // Init with zero!!
665
666 return comp[pos];
667}
668
669
670// matrix-matrix assignment
671template <typename T>
673{
674 if (this == &m)
675 return *this;
676 destroy(); copy(m);
677 return *this;
678}
679
680
681// KG 19980721
682template <typename T>
684{
685 if (n_rows != cm.n_rows || n_cols != cm.n_cols)
686 return false;
687 for (unsigned int i = 0; i < n_rows; ++i) {
688 for (unsigned int j = 0; j < n_cols; ++j) {
689 if (operator()(i,j) != cm(i,j))
690 return false;
691 }
692 }
693 return true;
694}
695
696
697// matrix negation
698template <typename T>
700{
702 for (unsigned int i = 0; i < n_size; ++i)
703 res.comp[i] = -comp[i];
704 return res; // avoid copy constructor soon?
705}
706
707
708#if 0
710template <typename T>
712{
713 TVector<T> res(n_rows);
714 BCHK(n_cols != v.size(), CSCMatErr, Dimension conflict, v.size(), res);
715 for (unsigned row = 0; row <= n_rows; ++row) {
716 REGISTER T val(0.0), corr(0.0);
717 /* Search for all elements in row row */
718 for (unsigned int c = 0; c < n_cols; ++c) {
719 const REGISTER T el = this->get(row, c);
720 if (el) {
721 const REGISTER T y = el*v(c);
722 const REGISTER T t = val+y;
723 corr += (t-val) - y;
724 val = t;
725 }
726 }
727 res.setval(row) = val - corr;
728 }
729 return res;
730}
731#else
732template <typename T>
734{
735 Vector<T> res((T)0, n_rows), corr((T)0, n_rows);
736 BCHK(n_cols != v.size(), CSCMatErr, Dimension conflict, v.size(), res);
737 for (unsigned int i = 0; i < n_cols; ++i) {
738 for (unsigned int j = pcol[i]; j < pcol[i+1]; ++j) {
739 const unsigned row = irow[j];
740 const REGISTER T y = comp[j]*v(i);
741 const REGISTER T t = res.get(row) + y;
742 corr.setval(row) += (t-res.get(row)) - y;
743 res.setval(row) = t;
744 }
745 }
746 return res-corr;
747}
748#endif
749
750// matrix-vector multiplication
751template <typename T>
753{
754 if (do_exactsum2())
755 return this->cscm_vec_mul_exact(v);
757 BCHK(n_cols != v.size(), CSCMatErr, Dimension conflict, v.size(), res);
758 /* Exact summing not possible */
759 for (unsigned int i = 0; i < n_cols; ++i) {
760 for (unsigned int j = pcol[i]; j < pcol[i+1]; ++j)
761 res.setval(irow[j]) += comp[j]*v(i);
762 }
763 return res;
764}
765
766template <typename T>
768{
769 Vector<T> v(tv);
770 return (this->operator * (v));
771}
772
773/* TODO: Implement exact summation for CSCM*TSV and MatVecMult */
774template <typename T>
776{
778 BCHK(n_cols != tsv.size(), CSCMatErr, Dimension conflict, tsv.size(), res);
779 /* Exact summing not possible */
780 for (unsigned int i = 0; i < n_cols; ++i) {
781 for (unsigned int j = pcol[i]; j < pcol[i+1]; ++j)
782 res.setval(irow[j]) += comp[j] * tsv.get(i);
783 }
784 tsv.destroy ();
785 return res;
786}
787
788
789// matrix-vector multiplication
790template <typename T>
792{
793 BCHKNR(n_cols != v.size(), CSCMatErr, Dimension conflict, v.size());
794 res = T(0);
795 /* Exact summing not possible */
796 for (unsigned int i = 0; i < n_cols; ++i) {
797 for (unsigned int j = pcol[i]; j < pcol[i+1]; ++j)
798 res.setval(irow[j]) += comp[j]*v(i);
799 }
800}
801
802INST(template <typename T> class CSCMatrix<T> friend void MatVecMult(Vector<T>&, const CSCMatrix<T>&, const Vector<T>&);)
803template <typename T>
804inline void MatVecMult(Vector<T>& res, const CSCMatrix<T>& m, const Vector<T>& v)
805{ m.MatVecMult (res, v); }
806
807
808template <typename T>
809void CSCMatrix<T>::MatVecMult (T* v, T* res) // res <- A*v , for ARPACK++
810{
811 STD__ cout << "CSCMatrix<T>::MatVecMult (T*, T*): Not tested!!" << STD__ endl;
812 // No Boundcheck possible!!
813 unsigned int i;
814 for (i = 0; i < n_cols; ++i)
815 res[i] = 0;
816 /* Exact summing not possible */
817 for (i = 0; i < n_cols; ++i) {
818 for (unsigned int j = pcol[i]; j < pcol[i+1]; ++j)
819 res[irow[j]] += comp[j]*v[i];
820 }
821}
822
823
824// transpose-vector multiplication
825template <typename T>
827{
829 BCHK(n_rows != v.size(), CSCMatErr, Dimension conflict, v.size(), res);
830 if (do_exactsum2()) {
831 for (unsigned int i = 0; i < n_cols; ++i) {
832 REGISTER T val = 0, corr = 0;
833 for (unsigned int j = pcol[i]; j < pcol[i+1]; ++j) {
834 const REGISTER T y = comp[j]*v(irow[j]);
835 const REGISTER T t = val + y;
836 corr += (t-val) - y;
837 val = t;
838 }
839 res.set(val-corr, i);
840 }
841 } else {
842 for (unsigned int i = 0; i < n_cols; ++i) {
843 T val = 0;
844 for (unsigned int j = pcol[i]; j < pcol[i+1]; ++j)
845 val += comp[j]*v(irow[j]);
846 res.set(val, i);
847 }
848 }
849 return res;
850}
851
852
853template <typename T>
855{
856 Vector<T> v(tv);
857 return(this->transMult(v));
858}
859
860
861template <typename T>
863{
864 Vector<T> v(tsv);
865 return(this->transMult(v));
866}
867
868
869// matrix-scalar multiplication
870template <typename T>
872{
873 CSCMatrix<T> res(*this);
874 for (unsigned int i = 0; i < n_size; ++i)
875 res.comp[i] = comp[i]*z;
876 return res; // JENS: avoid copy constructor soon!
877}
878
879
880// matrix-scalar multiplication
881template <typename T>
883{
884 CSCMatrix<T> res(*this);
885 for (unsigned int i = 0; i < n_size; ++i)
886 res.comp[i] = z*comp[i];
887 return res;
888}
889
890
891INST(template <typename T> class CSCMatrix friend CSCMatrix<T> operator* (const T&, const CSCMatrix<T>&);)
892template <typename T>
893inline CSCMatrix<T> operator* (const T& z, const CSCMatrix<T>& m)
894{ return m.mult (z); }
895
896// fast in-place matrix-scalar multiplication
897template <typename T>
899{
900 for (unsigned int i = 0; i < n_size; ++i)
901 comp[i] *= z;
902 return *this;
903}
904
905// matrix-scalar division
906template <typename T>
908{
909 CSCMatrix<T> res(*this);
910 for (unsigned int i = 0; i < n_size; ++i)
911 res.comp[i] = comp[i]/z;
912 return res; // JENS: avoid copy constructor soon!
913}
914
915
916// fast in-place matrix-scalar division
917template <typename T>
919{
920 for (unsigned int i = 0; i < n_size; ++i)
921 comp[i] /= z;
922 return *this;
923}
924
925
926// allocate memory space
927// used by constructor and resize()
928template <typename T>
929void CSCMatrix<T>::allocate(unsigned int rows, unsigned int columns,
930 unsigned int nnzeros)
931{
932 BCHKNR(rows <= 0 || columns <= 0, CSCMatErr, Zero space allocating, 0);
933 n_rows = rows;
934 n_cols = columns;
935 pcol = new unsigned int[n_cols+1];
936 BCHKNR(pcol==NULL, CSCMatErr, Out of memory, n_cols);
937 irow = new unsigned int[nnzeros];
938 BCHKNR(irow==NULL, CSCMatErr, Out of memory, nnzeros);
939 comp = new T[nnzeros];
940 BCHKNR(comp==NULL, CSCMatErr, Out of memory, nnzeros);
941
942 unsigned int i;
943 for (i = 0; i < n_cols+1; ++i)
944 pcol[i] = 0;
945 for (i = 0; i < nnzeros; ++i)
946 irow[i] = 0;
947
948 n_size = 0;
949 n_max_size = nnzeros;
950}
951
952
953// delete memory space
954// used by "=" operator and destructor
955template <typename T>
957{
958 if (irow != NULL)
959 delete[] irow;
960 if (comp != NULL)
961 delete[] comp;
962 if (pcol != NULL)
963 delete[] pcol;
964
965 n_size = 0;
966 n_max_size = 0;
967}
968
969
970// copy to new memory space
971// used by "=" operator and copy constructor
972template <typename T>
974{
975 n_rows = m.n_rows;
976 n_cols = m.n_cols;
977 n_size = m.n_size;
978 n_max_size = m.n_size;
979 pcol = new unsigned int[n_cols+1];
980 BCHKNR(pcol==NULL, CSCMatErr, Out of memory, n_cols);
981 irow = new unsigned int[n_size];
982 BCHKNR(pcol==NULL, CSCMatErr, Out of memory, n_size);
983 comp = new T[n_size];
984 BCHKNR(comp==NULL, CSCMatErr, Out of memory, n_size);
985
986 for (unsigned int j = 0; j < n_cols+1; ++j)
987 pcol[j] = m.pcol[j];
988 for (unsigned int i = 0; i < n_size; ++i) {
989 irow[i] = m.irow[i];
990 comp[i] = m.comp[i];
991 }
992}
993
994
995
996// set all elements to a certain value
997template <typename T>
999{
1000 if (val == (T)0)
1001 return this->clear();
1002 //for (unsigned int i = 0; i < n_size; ++i)
1003 // comp[i] = val;
1004 for (unsigned c = 0; c < columns(); ++c)
1005 for (unsigned r = 0; r < rows(); ++r)
1006 (*this)(r, c) = val;
1007 return *this;
1008}
1009
1010// clear all elements
1011template <typename T>
1013{ CSTD__ memset (comp, 0, sizeof(T) * n_size); return *this; }
1014
1015template <typename T>
1017{
1018 BCHK (n_rows != n_cols, CSCMatErr, setunit is meaningless for non-square matrices, n_cols, *this);
1019 clear ();
1020 for (unsigned int i = 0; i < n_cols; i++)
1021 setval (val, i, i);
1022 return *this;
1023}
1024
1025
1026template <typename T>
1028{
1029 SWAP(n_rows, m.n_rows); SWAP(n_cols, m.n_cols);
1030 SWAP(n_size, m.n_size); SWAP(n_max_size, m.n_max_size);
1031 SWAP(pcol, m.pcol); SWAP(irow, m.irow);
1032 SWAP(comp, m.comp);
1033 return *this;
1034}
1035
1036template <typename T>
1038{
1039 CSCMatrix<T> m(columns(), rows());
1040 for (unsigned r = 0; r < rows(); ++r)
1041 for (unsigned c = 0; c < columns(); ++c)
1042 m(c,r) = (*this)(r,c);
1043 return m;
1044}
1045
1046template <typename T>
1048{
1049 CSCMatrix<T> tr(this->transposed_copy());
1050 return swap(tr);
1051}
1052
1053INST(template <typename T> class CSCMatrix friend CSCMatrix<T> transpose(const CSCMatrix<T>& cscm);)
1054template <typename T>
1056{
1057 return cscm.transposed_copy();
1058}
1059
1061
1062#endif /* TBCI_CSCMATRIX_H */
long int Vector< T > & index
Definition LM_fit.h:69
const Vector< T > const Vector< T > & x
Definition LM_fit.h:97
int i
Definition LM_fit.h:71
doublereal y
Definition TOMS_707.C:27
#define NULL
Definition basics.h:250
#define CSTD__
Definition basics.h:340
#define STD__
Definition basics.h:338
#define HOT
Definition basics.h:495
#define BCHK(cond, exc, txt, ind, rtval)
Definition basics.h:575
#define FGDT
Definition basics.h:145
#define MIN(a, b)
Definition basics.h:655
#define MIN_ALIGN2
Definition basics.h:424
#define NAMESPACE_END
Definition basics.h:323
#define INST(x)
Definition basics.h:238
#define EXPCHK(cond, exc, txt, ind, rtval)
Definition basics.h:630
#define NAMESPACE_TBCI
Definition basics.h:317
#define BCHKNR(cond, exc, txt, ind)
Definition basics.h:586
#define UNLIKELY(expr)
Definition basics.h:101
#define REGISTER
Definition basics.h:108
#define MATH__
Definition basics.h:339
#define EXPCHKNR(cond, exc, txt, ind)
Definition basics.h:631
#define T
Definition bdmatlib.cc:20
The class BdMatrix is an implementation to store and do operations on sparse Matrices with a band str...
unsigned int rows() const
number of rows
unsigned int columns() const
number of columns
CSCMatErr(const CSCMatErr &ce)
Definition cscmatrix.h:39
virtual ~CSCMatErr()
Definition cscmatrix.h:41
CSCMatErr(const char *t, const long i=0)
Definition cscmatrix.h:37
exception class: Use MatErr from matrix.h
Definition cscmatrix.h:66
CSCMatrix(const Matrix< T > &m, const double tol=0)
Definition cscmatrix.h:89
CSCMatrix< T > & transpose()
Definition cscmatrix.h:1047
CSCMatrix< T > & fill(const T &)
set all defined element to a val
Definition cscmatrix.h:998
CSCMatrix< T > operator+(const CSCMatrix< T > &) const
Addition.
Definition cscmatrix.h:385
F_TMatrix< T > operator*(const CSCMatrix< T > &) const
CSCMatrix * CSCMatrix.
Definition cscmatrix.h:354
CSCMatrix< T > operator-() const
matrix negation
Definition cscmatrix.h:699
CSCMatrix< T > & clear()
set all elements defined to zero
Definition cscmatrix.h:1012
CSCMatrix()
constructors
Definition cscmatrix.h:75
unsigned int n_cols
Definition cscmatrix.h:215
T element_type
Definition cscmatrix.h:71
unsigned int * pcol
Definition cscmatrix.h:216
CSCMatrix(const F_Matrix< T > &m, const double tol=0)
Definition cscmatrix.h:85
unsigned int n_max_size
Definition cscmatrix.h:215
F_TMatrix< T > mult1(const F_Matrix< T > &) const
Calculate *this * m1 (smart version).
Definition cscmatrix.h:332
void MatVecMult(Vector< T > &res, const Vector< T > &v) const HOT
for friend void MatVecMult FGD (Vector<T>& res, const CSCMatrix<T>& m, const Vector<T>& v);
Definition cscmatrix.h:791
F_TMatrix< T > multf(const F_Matrix< T > &) const
matrix-matrix multiplication
Definition cscmatrix.h:265
unsigned int n_size
Definition cscmatrix.h:215
void copy(const CSCMatrix< T > &m)
Definition cscmatrix.h:973
CSCMatrix(const T &val, const unsigned int rows, const unsigned int columns, const unsigned int nnzeros=1)
Definition cscmatrix.h:76
CSCMatrix< T > & resize(const unsigned int newRows, const unsigned int newColumns, const unsigned int nnzeros=1)
change matrix dimensions
Definition cscmatrix.h:561
CSCMatrix< T > & setval(const T &z, unsigned int row, unsigned int column)
element access (write)
Definition cscmatrix.h:640
unsigned int * irow
Definition cscmatrix.h:217
CSCMatrix< T > transposed_copy() const
Inefficient! Use transMult if possible.
Definition cscmatrix.h:1037
bool operator==(const CSCMatrix< T > &m) const
matrix-matrix assignment and comparison
Definition cscmatrix.h:683
CSCMatrix< T > & operator/=(const T &z)
Definition cscmatrix.h:918
CSCMatrix< T > & setunit(const T &=(T) 1)
set CSCMatrix to val times the unit matrix
Definition cscmatrix.h:1016
unsigned int size() const
Definition cscmatrix.h:103
void do_export(MatType &M)
Export operation.
Definition cscmatrix.h:152
CSCMatrix< T > & operator*=(const T &z)
Definition cscmatrix.h:898
CSCMatrix(const CSCMatrix< T > &m)
Definition cscmatrix.h:83
void insert(const unsigned int column, const unsigned int pos)
Definition cscmatrix.h:591
CSCMatrix< T > & swap(CSCMatrix< T > &)
Definition cscmatrix.h:1027
TVector< T > transMult(const Vector< T > &v) const HOT
transpose-vector multiplication
Definition cscmatrix.h:826
tbci_traits< T >::const_refval_type get(unsigned int row, unsigned int column) const HOT
element access (read)
Definition cscmatrix.h:455
unsigned int columns() const
Definition cscmatrix.h:102
static const char * mat_info()
allow instantiation (Matrix_Sig)
Definition cscmatrix.h:97
CSCMatrix(const unsigned int rows, const unsigned int columns, const unsigned int nnzeros=1)
Definition cscmatrix.h:80
T aligned_value_type TALIGN(MIN_ALIGN2)
Definition cscmatrix.h:72
F_TMatrix< T > mult(const F_Matrix< T > &) const
Calculate *this * m1 (dumb version).
Definition cscmatrix.h:301
CSCMatrix(const BdMatrix< T > &m, const double tol=0)
Definition cscmatrix.h:87
CSCMatrix< T > & operator=(const CSCMatrix< T > &m)
Definition cscmatrix.h:672
TVector< T > cscm_vec_mul_exact(const Vector< T > &V) const
Definition cscmatrix.h:733
unsigned int n_rows
Storage format: pcol holds the offsets of each column; the length of the column c is pcol[c+1]-pcol[c...
Definition cscmatrix.h:215
bool operator!=(const CSCMatrix< T > &m) const
Definition cscmatrix.h:135
CSCMatrix< T > & do_import(const MatType &M)
Import operation, automatic conversion to sparse matrix.
Definition cscmatrix.h:142
T * dataPointer()
Definition cscmatrix.h:129
CSCMatrix< T > operator/(const T &z) const
Definition cscmatrix.h:907
unsigned int rows() const
query matrix dimensions
Definition cscmatrix.h:101
unsigned int * rowIndexPointer()
Definition cscmatrix.h:131
void allocate(unsigned int rows, unsigned int columns, unsigned int nnzeros=1)
Definition cscmatrix.h:929
tbci_traits< T >::const_refval_type operator()(unsigned int row, unsigned int column) const HOT
element access (read)
void destroy()
Definition cscmatrix.h:956
unsigned int * columnPointer()
Definition cscmatrix.h:130
unsigned int columns() const
Definition f_matrix.h:1213
unsigned int rows() const
Definition f_matrix.h:1214
Temporary Base Class (non referable!) (acc.
Definition f_matrix.h:71
Temporary object for scaled matrices.
Definition f_matrix.h:784
void destroy()
Definition f_matrix.h:914
unsigned int row
Definition f_matrix.h:787
unsigned int col
Definition f_matrix.h:787
T get(const unsigned r, const unsigned c) const
Definition f_matrix.h:829
unsigned int rows() const
TVector< T > transMult(const Vector_Sig< T > &) const
Matrix_Sig< T > & clear()
unsigned int columns() const
NumErr()
Definition except.h:65
unsigned int rows() const
number of rows
Definition matrix.h:317
unsigned int columns() const
number of columns
Definition matrix.h:315
void destroy() const
Definition vector.h:1172
T get(const unsigned long i) const HOT
Definition vector.h:1072
unsigned long size() const
Definition vector.h:1120
Temporary Base Class Idiom: Class TVector is used for temporary variables.
Definition vector.h:73
unsigned long size() const
Definition vector.h:104
T & setval(const unsigned long i) const
Definition vector.h:192
double fabs(const TBCI__ cplx< T > &c)
Definition cplx.h:746
STD__ ostream & operator<<(STD__ ostream &stream, const CSCMatrix< T > &m)
Definition cscmatrix.h:248
CSCMatrix< T > transpose(const CSCMatrix< T > &cscm)
Definition cscmatrix.h:1055
void MatVecMult(Vector< T > &res, const CSCMatrix< T > &m, const Vector< T > &v)
Definition cscmatrix.h:804
return c
Definition f_matrix.h:760
return F_TMatrix< T >(tm)
const unsigned TMatrix< T > * res
unsigned int do_exactsum2()
Definition tbci_param.h:151