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 
32 class CSCMatErr : public NumErr
33 {
34  public:
36  : NumErr("Error in CSCMatrix") {}
37  CSCMatErr(const char* t, const long i = 0)
38  : NumErr(t, i) {}
39  CSCMatErr(const CSCMatErr& ce)
40  : NumErr(ce) {}
41  virtual ~CSCMatErr() throw() {}
42 };
43 #endif
44 
45 #ifdef PRAGMA_I
46 # pragma interface "cscmatrix.h"
47 #endif
48 
49 template <typename T> class CSCMatrix;
50 template <typename T> class F_TSMatrix;
51 template <typename T> class F_TMatrix;
52 template <typename T> class F_Matrix;
53 
54 //----------
55 // CSCMatrix
56 //----------
64 template <typename T>
65 class CSCMatrix : public Matrix_Sig<T>
66 {
67  friend class F_TSMatrix<T>;
68 
69 public:
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 
112  CSCMatrix<T>& clear ();
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 
162  CSCMatrix<T> operator - () const;
163 
165  CSCMatrix<T> operator + (const CSCMatrix<T>&) const;
166  CSCMatrix<T> operator - (const CSCMatrix<T>&) const;
167 
169  F_TMatrix<T> operator * (const CSCMatrix<T>&) const;
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 
214 protected:
215  unsigned int n_rows, n_cols, n_size, n_max_size;
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 
247 template <typename T>
248 STD__ 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 
264 template <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 
294 INST(template<typename T> class CSCMatrix<T> friend F_TMatrix<T> operator* (const F_Matrix<T>&, const CSCMatrix<T>&);)
295 template <typename T>
296 inline F_TMatrix<T> operator * (const F_Matrix<T>& m1, const CSCMatrix<T>& m2)
297 { return m2.multf (m1); }
298 
300 template <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 
331 template <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 
347 INST(template<typename T> class CSCMatrix<T> friend F_TMatrix<T> operator* (const CSCMatrix<T>&, const F_Matrix<T>&);)
348 template <typename T>
349 inline F_TMatrix<T> operator * (const CSCMatrix<T>& m1, const F_Matrix<T>& m2)
350 { return m1.mult1 (m2); }
351 
353 template <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 
384 template <typename T>
386 {
387  CSCMatrix<T> res(rows(), columns(), size());
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 
417 template <typename T>
419 {
420  CSCMatrix<T> res(rows(), columns(), size());
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 
453 template <typename T>
454 inline 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 
538 template <typename T>
539 inline 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 
548 template <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 
560 template <typename T>
561 CSCMatrix<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 
590 template <typename T>
591 void 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)
639 template <typename T>
640 CSCMatrix<T>& CSCMatrix<T>::setval(const T& z, unsigned int row, unsigned int column)
641 {
642  setval(row, column) = z;
643  return *this;
644 }
645 
647 template <typename T>
648 T& 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
671 template <typename T>
673 {
674  if (this == &m)
675  return *this;
676  destroy(); copy(m);
677  return *this;
678 }
679 
680 
681 // KG 19980721
682 template <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
698 template <typename T>
700 {
701  CSCMatrix<T> res(n_rows, n_cols);
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
709 template <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
732 template <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
751 template <typename T>
753 {
754  if (do_exactsum2())
755  return this->cscm_vec_mul_exact(v);
756  TVector<T> res(0, n_rows);
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 
766 template <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 */
774 template <typename T>
776 {
777  TVector<T> res(0, n_rows);
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
790 template <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 
802 INST(template <typename T> class CSCMatrix<T> friend void MatVecMult(Vector<T>&, const CSCMatrix<T>&, const Vector<T>&);)
803 template <typename T>
804 inline void MatVecMult(Vector<T>& res, const CSCMatrix<T>& m, const Vector<T>& v)
805 { m.MatVecMult (res, v); }
806 
807 
808 template <typename T>
809 void 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
825 template <typename T>
827 {
828  TVector<T> res(n_cols);
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 
853 template <typename T>
855 {
856  Vector<T> v(tv);
857  return(this->transMult(v));
858 }
859 
860 
861 template <typename T>
863 {
864  Vector<T> v(tsv);
865  return(this->transMult(v));
866 }
867 
868 
869 // matrix-scalar multiplication
870 template <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
881 template <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 
891 INST(template <typename T> class CSCMatrix friend CSCMatrix<T> operator* (const T&, const CSCMatrix<T>&);)
892 template <typename T>
893 inline CSCMatrix<T> operator* (const T& z, const CSCMatrix<T>& m)
894 { return m.mult (z); }
895 
896 // fast in-place matrix-scalar multiplication
897 template <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
906 template <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
917 template <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()
928 template <typename T>
929 void 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
955 template <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
972 template <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
997 template <typename T>
998 inline CSCMatrix<T>& CSCMatrix<T>::fill (const T& val)
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
1011 template <typename T>
1013 { CSTD__ memset (comp, 0, sizeof(T) * n_size); return *this; }
1014 
1015 template <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 
1026 template <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 
1036 template <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 
1046 template <typename T>
1048 {
1049  CSCMatrix<T> tr(this->transposed_copy());
1050  return swap(tr);
1051 }
1052 
1053 INST(template <typename T> class CSCMatrix friend CSCMatrix<T> transpose(const CSCMatrix<T>& cscm);)
1054 template <typename T>
1056 {
1057  return cscm.transposed_copy();
1058 }
1059 
1061 
1062 #endif /* TBCI_CSCMATRIX_H */
BdMatrix< T > transpose(BdMatrix< T > &mat)
Definition: band_matrix.h:1393
TVector< T > transMult(const Vector_Sig< T > &) const
#define BCHKNR(cond, exc, txt, ind)
Definition: basics.h:586
const Vector< T > const Vector< T > & x
Definition: LM_fit.h:97
unsigned int row
Definition: f_matrix.h:787
CSCMatrix< T > & setval(const T &z, unsigned int row, unsigned int column)
element access (write)
Definition: cscmatrix.h:640
T & setval(const unsigned long i) const
Definition: vector.h:192
~CSCMatrix()
Definition: cscmatrix.h:91
CSCMatErr(const char *t, const long i=0)
Definition: cscmatrix.h:37
bool operator!=(const CSCMatrix< T > &m) const
Definition: cscmatrix.h:135
unsigned int columns() const
Definition: f_matrix.h:1213
The class BdMatrix is an implementation to store and do operations on sparse Matrices with a band str...
Definition: band_matrix.h:103
unsigned int * irow
Definition: cscmatrix.h:217
CSCMatrix(const BdMatrix< T > &m, const double tol=0)
Definition: cscmatrix.h:87
bool operator==(const CSCMatrix< T > &m) const
matrix-matrix assignment and comparison
Definition: cscmatrix.h:683
#define REGISTER
Definition: basics.h:108
return c
Definition: f_matrix.h:760
double fabs(const int a)
Definition: basics.h:1215
#define MIN(a, b)
Definition: basics.h:655
F_TMatrix< T > mult(const F_Matrix< T > &) const
Calculate *this * m1 (dumb version)
Definition: cscmatrix.h:301
#define NAMESPACE_TBCI
Definition: basics.h:317
unsigned int n_size
Definition: cscmatrix.h:215
TVector< T > transMult(const Vector< T > &v) const HOT
transpose-vector multiplication
Definition: cscmatrix.h:826
CSCMatrix(const unsigned int rows, const unsigned int columns, const unsigned int nnzeros=1)
Definition: cscmatrix.h:80
exception base class for the TBCI NumLib
Definition: except.h:58
unsigned int * columnPointer()
Definition: cscmatrix.h:130
T * comp
Definition: cscmatrix.h:218
unsigned int rows() const
Common interface definition (signature) for all Matrices.
Definition: matrix_sig.h:45
#define MIN_ALIGN2
Definition: basics.h:424
F_TMatrix< T > multf(const F_Matrix< T > &) const
matrix-matrix multiplication
Definition: cscmatrix.h:265
TVector< T > cscm_vec_mul_exact(const Vector< T > &V) const
Definition: cscmatrix.h:733
unsigned int columns() const
number of columns
Definition: band_matrix.h:249
#define BCHK(cond, exc, txt, ind, rtval)
Definition: basics.h:575
#define NULL
Definition: basics.h:250
CSCMatErr(const CSCMatErr &ce)
Definition: cscmatrix.h:39
CSCMatrix< T > transposed_copy() const
Inefficient! Use transMult if possible.
Definition: cscmatrix.h:1037
CSCMatrix< T > & operator*=(const T &z)
Definition: cscmatrix.h:898
unsigned int rows() const
query matrix dimensions
Definition: cscmatrix.h:101
F_TMatrix< T > mult1(const F_Matrix< T > &) const
Calculate *this * m1 (smart version)
Definition: cscmatrix.h:332
#define UNLIKELY(expr)
Definition: basics.h:101
CSCMatrix< T > & do_import(const MatType &M)
Import operation, automatic conversion to sparse matrix.
Definition: cscmatrix.h:142
unsigned int size() const
Definition: cscmatrix.h:103
void MatVecMult(Vector< T > &res, const Vector< T > &v) const HOT
for friend void MatVecMult FGD (Vector&lt;T&gt;&amp; res, const CSCMatrix&lt;T&gt;&amp; m, const Vector&lt;T&gt;&amp; v); ...
Definition: cscmatrix.h:791
CSCMatrix< T > & operator=(const CSCMatrix< T > &m)
Definition: cscmatrix.h:672
T & set(const T &val, const unsigned long i) const
Definition: vector.h:194
static const char * mat_info()
allow instantiation (Matrix_Sig)
Definition: cscmatrix.h:97
Temporary object for scaled matrices.
Definition: cscmatrix.h:50
Matrix_Sig< T > & clear()
F_TMatrix< T > operator*(const CSCMatrix< T > &) const
CSCMatrix * CSCMatrix.
Definition: cscmatrix.h:354
CSCMatrix(const CSCMatrix< T > &m)
Definition: cscmatrix.h:83
const Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > & z
Definition: LM_fit.h:172
void destroy() const
Definition: vector.h:1172
CSCMatrix< T > & transpose()
Definition: cscmatrix.h:1047
CSCMatrix< T > operator/(const T &z) const
Definition: cscmatrix.h:907
void setval(const T val, const unsigned int r, const unsigned int c)
Definition: f_matrix.h:206
unsigned int n_cols
Definition: cscmatrix.h:215
T element_type
Definition: cscmatrix.h:71
T get(const unsigned long i) const HOT
Definition: vector.h:1072
for(REGISTER T *p1=c.vec,*p2=b.vec;p1< c.endvec;p1++, p2++)*p1
unsigned long _bin_search(const T *vec, T el, unsigned long start, unsigned long end)
Search for an element el in a sorted vector between start and end-1, returns (unsigned long)-1 if ele...
Definition: basics.h:992
#define CSTD__
Definition: basics.h:340
CSCMatrix< T > & clear()
set all elements defined to zero
Definition: cscmatrix.h:1012
unsigned int columns() const
number of columns
Definition: matrix.h:315
unsigned int * rowIndexPointer()
Definition: cscmatrix.h:131
void insert(const unsigned int column, const unsigned int pos)
Definition: cscmatrix.h:591
long int Vector< T > & index
Definition: LM_fit.h:69
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
CSCMatrix(const F_Matrix< T > &m, const double tol=0)
Definition: cscmatrix.h:85
void destroy()
Definition: cscmatrix.h:956
CSCMatrix< T > operator+(const CSCMatrix< T > &) const
Addition.
Definition: cscmatrix.h:385
void allocate(unsigned int rows, unsigned int columns, unsigned int nnzeros=1)
Definition: cscmatrix.h:929
CSCMatrix< T > & operator/=(const T &z)
Definition: cscmatrix.h:918
void copy(const CSCMatrix< T > &m)
Definition: cscmatrix.h:973
CSCMatrix(const Matrix< T > &m, const double tol=0)
Definition: cscmatrix.h:89
unsigned int columns() const
unsigned int rows() const
Definition: f_matrix.h:1214
Matrix_Sig< T > & setunit(const T &=(T) 1)
CSCMatrix< T > & fill(const T &)
set all defined element to a val
Definition: cscmatrix.h:998
unsigned int rows() const
number of rows
Definition: band_matrix.h:247
void SWAP(T &a, T &b)
SWAP function Note: We could implement a swap function without temporaries: a -= b b += a a -= b a = ...
Definition: basics.h:813
tbci_traits< T >::const_refval_type operator()(unsigned int row, unsigned int column) const HOT
element access (read)
unsigned long size() const
Definition: vector.h:104
Definition: bvector.h:49
#define INST(x)
Definition: basics.h:238
unsigned int * pcol
Definition: cscmatrix.h:216
#define EXPCHK(cond, exc, txt, ind, rtval)
Definition: basics.h:630
CSCMatrix()
constructors
Definition: cscmatrix.h:75
Temporary Base Class (non referable!) (acc.
Definition: bvector.h:50
int i
Definition: LM_fit.h:71
tbci_traits< T >::const_refval_type get(unsigned int row, unsigned int column) const HOT
element access (read)
Definition: cscmatrix.h:455
exception class: Use MatErr from matrix.h
Definition: cscmatrix.h:49
#define STD__
Definition: basics.h:338
void destroy()
Definition: f_matrix.h:914
unsigned int n_max_size
Definition: cscmatrix.h:215
Temporary Base Class Idiom: Class TVector is used for temporary variables.
Definition: bvector.h:52
unsigned long size() const
Definition: vector.h:1120
void MatVecMult(Vector< T > &res, const CRMatrix< T > &m, const Vector< T > &v)
Definition: crmatrix.h:532
#define NAMESPACE_END
Definition: basics.h:323
T value_type
Definition: cscmatrix.h:70
Definition: bvector.h:54
const Vector< T > Vector< T > Vector< T > Vector< T > & y
Definition: LM_fit.h:172
const Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > long int res
Definition: LM_fit.h:199
CSCMatrix< T > operator-() const
matrix negation
Definition: cscmatrix.h:699
T * dataPointer()
Definition: cscmatrix.h:129
void do_export(MatType &M)
Export operation.
Definition: cscmatrix.h:152
#define T
Definition: bdmatlib.cc:20
#define HOT
Definition: basics.h:495
unsigned int do_exactsum2()
Definition: tbci_param.h:151
unsigned int col
Definition: f_matrix.h:787
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 > & swap(CSCMatrix< T > &)
Definition: cscmatrix.h:1027
unsigned int columns() const
Definition: cscmatrix.h:102
#define EXPCHKNR(cond, exc, txt, ind)
Definition: basics.h:631
const Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > long int int char v
&lt; find minimun of func on grid with resolution res
Definition: LM_fit.h:205
#define MATH__
Definition: basics.h:339
T aligned_value_type TALIGN(MIN_ALIGN2)
Definition: cscmatrix.h:72
virtual ~CSCMatErr()
Definition: cscmatrix.h:41
CSCMatrix(const T &val, const unsigned int rows, const unsigned int columns, const unsigned int nnzeros=1)
Definition: cscmatrix.h:76
CSCMatrix< T > & setunit(const T &=(T) 1)
set CSCMatrix to val times the unit matrix
Definition: cscmatrix.h:1016
unsigned int rows() const
number of rows
Definition: matrix.h:317
CSCMatErr()
Definition: cscmatrix.h:35