TBCI Numerical high perf. C++ Library 2.8.0
crmatrix.h
Go to the documentation of this file.
1
5
6//----------------------------------------------------------
7// version 1.1, may 1998
8//
9// Marco Wiedenhaus
10// Basis: TCRMatrixComplex by Jens Lenge
11//
12// $Id: crmatrix.h,v 1.11.2.19 2019/05/28 11:13:02 garloff Exp $
13//----------------------------------------------------------
14
15#ifndef TBCI_CRMATRIX_H
16#define TBCI_CRMATRIX_H
17
18#include "tbci/matrix_sig.h"
19#include "tbci/matrix.h"
20#include "tbci/vector.h"
21
22//#include <cstdlib> // for use of null macro
23
24// Avoid -fguiding-decls
25#if !defined(NO_GD) && !defined(AUTO_DECL)
26# include "tbci/crmatrix_gd.h"
27#endif
28
29NAMESPACE_TBCI
30
31// exception class
32#ifndef TBCI_DISABLE_EXCEPT
33//# include "except.h" is included in basics.h
34
35class CRMatErr : public NumErr
36{
37 public:
38 CRMatErr()
39 : NumErr("Error in crmatrix library") {}
40 CRMatErr(const char* t, const long i = 0)
41 : NumErr(t, i) {}
42 CRMatErr(const CRMatErr& ce)
43 : NumErr(ce) {}
44 virtual ~CRMatErr() throw() {}
45};
46#endif
47
48
49#ifdef PRAGMA_I
50# pragma interface "crmatrix.h"
51#endif
52
53
54//----------
55// CRMatrix
56//----------
57
62template <typename T>
63class CRMatrix : public Matrix_Sig<T>
64{
65protected:
66 unsigned int n_rows, n_cols;
67 unsigned int* rowSize; // number of elements actually stored in each row
68 unsigned int** colIndex; // column index of each stored element
69 T** comp;
70 mutable T dummy;
71 // used by constructor and resize()
72 void allocate(unsigned int rows, unsigned int columns);
73 // used by destructor and "=" operator
74 void destroy();
75 // used by copy constr. and "="
76 void copy(const CRMatrix<T>& m);
77
78public:
79 typedef T value_type;
80 typedef T element_type;
81 typedef T aligned_value_type TALIGN(MIN_ALIGN2);
82
83 CRMatrix() { allocate(1,1); }
84 CRMatrix(unsigned int rows, unsigned int columns) { allocate(rows,columns); }
85 CRMatrix(const CRMatrix<T>& m) { copy(m); }
86 CRMatrix(const Matrix/*_Sig*/<T>& mat);
87 ~CRMatrix() { destroy(); }
88
89 // allow instantiation (Matrix_Sig)
90 /*virtual*/ static const char* mat_info () { return ("CRMatrix"); }
91
93 CRMatrix<T>& clear ();
95 CRMatrix<T>& fill (const T&);
97 CRMatrix<T>& setunit (const T& = (T)1);
98
100 inline const T& operator () (unsigned int row, unsigned int column) const;
101 inline const T& get (unsigned int row, unsigned int column) const
102 { return (*this) (row, column); }
103
104 CRMatrix<T>& setval (const T& z, unsigned int row, unsigned int column);
105 T& setval (unsigned int row, unsigned int column);
106 T& operator () (unsigned int row, unsigned int column)
107 { return setval (row, column); }
108
109 // rowwise access
110 TVector<T> get_row(unsigned int row) const;
111 void set_row(const Vector<T>& v, unsigned int row);
112 void set_row(const TVector<T>& tv, unsigned int row);
113 void set_row(const TSVector<T>& tsv, unsigned int row);
114 // replace a certain row using pre-compressed array storage
115 void set_row( T* value, unsigned int row, unsigned int count, unsigned int* index);
116
117 // query matrix dimensions
118 unsigned int rows() const {return n_rows;}
119 unsigned int columns() const {return n_cols;}
120 unsigned int size() const; // number of non-zero elements
121
122 // change matrix dimensions
123 inline void resize(unsigned int newRows, unsigned int newColumns);
124
125 // matrix-matrix assignment and comparison
126 CRMatrix<T>& operator= (const CRMatrix<T>& m);
127 bool operator == (const CRMatrix<T>& m) const;
128 bool operator != (const CRMatrix<T>& m) const { return !(this->operator == (m)); }
129
130 // matrix negation
131 CRMatrix<T> operator - () const;
132
133 // matrix - vector multiplication
134 TVector<T> operator * (const Vector<T>& v) const;
135 TVector<T> operator * (const TVector<T>& tv) const;
136 TVector<T> operator * (TSVector<T>& tsv) const;
137
138 // For friend void MatVecMult FGD (Vector<T>& res, const CRMatrix<T>& m, const Vector<T>& v);
139 void MatVecMult (Vector<T>& res, /*ourself ,*/ const Vector<T>& v) const;
140 void MatVecMult (T* v, T* res); // res <- A*v , for ARPACK++
141
142 // transpose-vector multiplication
143 TVector<T> transMult(const Vector<T>& v) const;
144 TVector<T> transMult(const TVector<T>& tv) const;
145 TVector<T> transMult(const TSVector<T>& tsv) const;
146
147 CRMatrix<T> operator* (const T& z) const; // matrix-scalar multiplication
148 CRMatrix<T>& operator*= (const T& z); // fast in-place matrix-scalar multiplication
149 CRMatrix<T> operator/ (const T& z) const; // matrix-scalar division
150 CRMatrix<T>& operator/= (const T& z); // fast in-place matrix-scalar division
151
152 // for friend CrMatrix,T> operator* (const T&, const CRMatrix<T>&);
153 CRMatrix<T> mult (const T& z) const;
154
155 CRMatrix<T>& swap (CRMatrix<T>&);
157 CRMatrix<T> transposed_copy() const;
158 CRMatrix<T>& transpose();
159
160 // Output operations
161 friend STD__ ostream& operator<< FGD (STD__ ostream& stream, const CRMatrix<T>& m);
162};
163
164
165//-------------------------------
166// CRMatrix function definitions
167//-------------------------------
168
169
170template <typename T>
171CRMatrix<T>::CRMatrix(const Matrix/*_Sig*/<T>& mat)
172{
173 allocate (mat.rows(), mat.columns());
174 for (unsigned i = 0; i < n_rows; ++i)
175 for (unsigned j = 0; j < n_cols; ++j)
176 (*this) (i, j) = mat( i, j);
177}
178
179// element access (read)
180template <typename T>
181inline const T& CRMatrix<T>::operator() (unsigned int row, unsigned int column) const
182{
183 BCHK(row >= n_rows || column >= n_cols, CRMatErr, Illegal index, (row<column ? row : column), dummy = (T)0);
184 for (unsigned int j=0; j<rowSize[row]; j++)
185 if (colIndex[row][j] == column)
186 return comp[row][j];
187 return (dummy = (T)0);
188}
189
190
191// number of non-zero elements
192template <typename T>
193unsigned int CRMatrix<T>::size() const
194{
195 unsigned int erg=0;
196 for (unsigned int i=0; i<n_rows; ++i)
197 erg += rowSize[i];
198 return erg;
199}
200
201
202
203// change matrix dimensions
204template <typename T>
205inline void CRMatrix<T>::resize(unsigned int newRows, unsigned int newColumns)
206{
207 if (newRows != n_rows || newColumns != n_cols)
208 destroy();
209 allocate(newRows, newColumns);
210}
211
212
213// element access (write)
214template <typename T>
215CRMatrix<T>& CRMatrix<T>::setval(const T& z, unsigned int row, unsigned int column)
216{
217 BCHK(row >= n_rows, CRMatErr, Illegal row index, row, *this);
218 BCHK(column >= n_cols, CRMatErr, Illegal col index, column, *this);
219 unsigned int j;
220 for (j=0; j<rowSize[row]; ++j)
221 if (colIndex[row][j] == column) {
222 comp[row][j] = z;
223 return *this;
224 }
225
226 // element not found -- dynamically growing required...
227 if (z==0.0)
228 return *this; // no need for growing
229 // create new column...
230 unsigned int* newIndex = new unsigned int[rowSize[row]+1];
231 BCHK( newIndex==NULL, CRMatErr, Out of Memory, rowSize[row]+1, *this);
232 T* newComp = new T[rowSize[row]+1];
233 BCHK( newComp==NULL, CRMatErr, Out of Memory, rowSize[row]+1, *this);
234 // store "old" elements before new element
235 for (j=0; j<rowSize[row]; ++j)
236 if (colIndex[row][j] < column) {
237 newIndex[j] = colIndex[row][j];
238 newComp[j] = comp[row][j];
239 } else
240 break;
241 // store new element
242 newIndex[j] = column;
243 newComp[j] = z;
244 // store "old" elementes behind new element
245 for (unsigned int i=j; i<rowSize[row]; i++) {
246 newIndex[i+1] = colIndex[row][i];
247 newComp[i+1] = comp[row][i];
248 }
249 // kill old column
250 if (rowSize[row] != 0) {
251 delete[] colIndex[row];
252 delete[] comp[row];
253 }
254 rowSize[row]++;
255 colIndex[row] = newIndex;
256 comp[row] = newComp;
257 return *this;
258}
259
260template <typename T>
261T& CRMatrix<T>::setval(unsigned int row, unsigned int column)
262{
263 BCHK(row >= n_rows || column >= n_cols, CRMatErr, "Illegal index", (row<column ? row : column), comp[0][0]);
264 unsigned int j;
265 for (j=0; j<rowSize[row]; j++)
266 if (colIndex[row][j] == column)
267 {
268 return comp[row][j];
269 }
270
271 // element not found -- dynamically growing required...
272 // create new column...
273 unsigned int* newIndex = new unsigned int[rowSize[row]+1];
274 BCHK( newIndex==NULL, CRMatErr, Out of Memory, rowSize[row]+1, comp[0][0]);
275 T* newComp = new T[rowSize[row]+1];
276 BCHK( newComp==NULL, CRMatErr, Out of Memory, rowSize[row]+1, comp[0][0]);
277 // store "old" elements before new element
278 for (j=0; j<rowSize[row]; j++)
279 if (colIndex[row][j] < column)
280 {
281 newIndex[j] = colIndex[row][j];
282 newComp[j] = comp[row][j];
283 }
284 else break;
285 // store new element
286 unsigned int pos_store(j);
287 newIndex[j] = column;
288// newComp[j] = z;
289 // store "old" elementes behind new element
290 for (unsigned int i=j; i<rowSize[row]; i++)
291 {
292 newIndex[i+1] = colIndex[row][i];
293 newComp[i+1] = comp[row][i];
294 }
295 // kill old column
296 if (rowSize[row] != 0)
297 {
298 delete[] colIndex[row];
299 delete[] comp[row];
300 }
301 rowSize[row]++;
302 colIndex[row] = newIndex;
303 comp[row] = newComp;
304
305 // Init with zero!!
306 comp[row][pos_store]=0;
307 return comp[row][pos_store];
308}
309
310// extract a certain row
311template <typename T>
312TVector<T> CRMatrix<T>::get_row(unsigned int row) const
313{
314 TVector<T> RowVec(n_cols);
315 BCHK(row >= n_rows, CRMatErr, Illegal index, row, RowVec);
316 for (unsigned int j=0; j<n_cols; j++)
317 RowVec.set((*this)(row,j),j);
318 return RowVec;
319}
320
321// replace a certain row
322template <typename T>
323void CRMatrix<T>::set_row( const Vector<T>& v, unsigned int row)
324{
325 BCHKNR(row >= n_rows, CRMatErr, Illegal index, row);
326 unsigned int count = 0;
327 unsigned int j;
328 for (j=0; j<v.size(); j++)
329 if (v(j) != 0.0) count++;
330 if (rowSize[row] != count)
331 {
332 if (rowSize[row]!=0) // kill old row...
333 {
334 delete[] colIndex[row];
335 delete[] comp[row];
336 }
337 if ((rowSize[row]=count) != 0) // create new row...
338 {
339 colIndex[row] = new unsigned int[count];
340 BCHKNR(colIndex[row]==NULL, CRMatErr, Out of memory, count);
341 comp[row] = new T[row];
342 BCHKNR(comp[row]==NULL, CRMatErr, Out of memory, row);
343 }
344 else
345 {
346 colIndex[row] = NULL;
347 comp[row] = NULL;
348 }
349 }
350 count = 0;
351 for (j=0; j<v.size(); j++) // copy values...
352 if (v(j)!=0.0)
353 {
354 colIndex[row][count] = j;
355 comp[row][count] = v(j);
356 count++;
357 }
358}
359
360template <typename T>
361void CRMatrix<T>::set_row( const TVector<T>& tv, unsigned int row)
362{
363 Vector<T> v(tv);
364 set_row(v,row);
365}
366
367template <typename T>
368void CRMatrix<T>::set_row( const TSVector<T>& tsv, unsigned int row)
369{
370 Vector<T> v(tsv);
371 set_row(v,row);
372}
373
374// replace a certain row
375// using pre-compressed array storage
376template <typename T>
377void CRMatrix<T>::set_row( T* value,unsigned int row, unsigned int count, unsigned int* index)
378{
379 BCHKNR(row >= n_rows || count >= n_cols, CRMatErr, Illegal index, (row<count ? row : count));
380 if (rowSize[row] != count)
381 {
382 if (rowSize[row]!=0) // kill old row...
383 {
384 delete[] colIndex[row];
385 delete[] comp[row];
386 }
387 if ((rowSize[row]=count) != 0) // create new row...
388 {
389 colIndex[row] = new unsigned int[count];
390 BCHKNR(colIndex[row]==NULL, CRMatErr, Out of memory, count);
391 comp[row] = new T[row];
392 BCHKNR(comp[row]==NULL, CRMatErr, Out of memory, row);
393 }
394 else
395 {
396 colIndex[row] = NULL;
397 comp[row] = NULL;
398 }
399 }
400 for (unsigned int j=0; j<count; j++) // copy values...
401 {
402 colIndex[row][j] = index[j];
403 comp[row][j] = value[j];
404 }
405}
406
407// matrix-matrix assignment
408template <typename T>
409CRMatrix<T>& CRMatrix<T>::operator= (const CRMatrix<T>& m)
410{
411 if (this == &m) return *this;
412 if (n_rows==m.n_rows) // perform rowwise copy
413 {
414 for (unsigned int i=0; i<n_rows; i++)
415 {
416 if (rowSize[i]!=m.rowSize[i]) // different row sizes
417 {
418 if (rowSize[i] != 0) // kill old row...
419 {
420 delete[] colIndex[i];
421 delete[] comp[i];
422 }
423 if ((rowSize[i]=m.rowSize[i]) != 0) // create new row...
424 {
425 colIndex[i] = new unsigned int[rowSize[i]];
426 BCHK(colIndex[i]==NULL, CRMatErr, Out of memory, rowSize[i], *this);
427 comp[i] = new T[rowSize[i]];
428 BCHK(comp[i]==NULL, CRMatErr, Out of memory, rowSize[i], *this);
429 }
430 else
431 {
432 colIndex[i] = NULL;
433 comp[i] = NULL;
434 }
435 }
436 for (unsigned int j=0; j<rowSize[i]; j++) // copy values...
437 {
438 colIndex[i][j] = m.colIndex[i][j];
439 comp[i][j] = m.comp[i][j];
440 }
441 } // for each row...
442 return *this;
443 } // end of rowwise copy...
444 destroy();
445 copy(m);
446 return *this;
447}
448
449
450// KG 980721
451template <typename T>
452bool CRMatrix<T>::operator == (const CRMatrix<T>& cm) const
453{
454 if (n_rows != cm.n_rows || n_cols != cm.n_cols)
455 return false;
456 for (unsigned int i=0; i<n_rows; ++i)
457 for (unsigned int j=0; j<n_cols; ++j)
458 if ((*this)(i,j) != cm(i,j))
459 return false;
460 return true;
461}
462
463
464// matrix negation
465template <typename T>
466CRMatrix<T> CRMatrix<T>::operator - () const
467{
468 CRMatrix<T> res(n_rows, n_cols);
469 for (unsigned int i=0; i<n_rows; ++i)
470 for (unsigned int j=0; j<n_cols; ++j)
471 if ((*this)(i,j) != (T)0)
472 res.setval(-(*this)(i,j),i,j);
473 return res; // JENS: avoid copy constructor soon
474}
475
476
477// matrix-vector multiplication
478template <typename T>
479TVector<T> CRMatrix<T>::operator * (const Vector<T>& v) const
480{
481 TVector<T> res(n_rows);
482 BCHK(n_cols != v.size(), CRMatErr, Dimension conflict, v.size(), res);
483 for (unsigned int i=0; i<n_rows; ++i) {
484 T var = 0.0;
485 for (unsigned int j=0; j<rowSize[i]; ++j)
486 var += comp[i][j] * v(colIndex[i][j]);
487 res.set(var,i);
488 }
489 return res;
490}
491
492
493template <typename T>
494inline TVector<T> CRMatrix<T>::operator * (const TVector<T>& tv) const
495{
496 Vector<T> v(tv);
497 return (this->operator * (v));
498}
499
500
501template <typename T>
502inline TVector<T> CRMatrix<T>::operator * (TSVector<T>& tsv) const
503{
504 TVector<T> res(n_rows);
505 BCHK(n_cols != tsv.size(), CRMatErr, Dimension conflict, tsv.size(), res);
506 for (unsigned int i=0; i<n_rows; ++i) {
507 T var = 0.0;
508 for (unsigned int j=0; j<rowSize[i]; ++j)
509 var += comp[i][j] * tsv.get(colIndex[i][j]);
510 res.set(var,i);
511 }
512 tsv.destroy ();
513 return res;
514}
515
516
517// matrix-vector multiplication
518template <typename T>
519void CRMatrix<T>::MatVecMult(Vector<T>& res, const Vector<T>& v) const
520{
521 BCHKNR(n_cols != v.size(), CRMatErr, Dimension conflict, v.size());
522 for (unsigned int i=0; i<n_rows; ++i) {
523 T var = 0.0;
524 for (unsigned int j=0; j<rowSize[i]; ++j)
525 var += comp[i][j] * v(colIndex[i][j]);
526 res(i) = var;
527 }
528}
529
530INST(template <typename T> class CRMatrix friend void MatVecMult(Vector<T>&, const CRMatrix<T>&, const Vector<T>&);)
531template <typename T>
532inline void MatVecMult(Vector<T>& res, const CRMatrix<T>& m, const Vector<T>& v)
533{ m.MatVecMult (res, v); }
534
535template <typename T>
536void CRMatrix<T>::MatVecMult (T* v, T* res) // res <- A*v , for ARPACK++
537{
538 // No Boundcheck possible!!
539 for (unsigned int i=0; i<this->n_rows; ++i) {
540 T var = 0.0;
541 for (unsigned int j=0; j<this->rowSize[i]; ++j)
542 var += this->comp[i][j] * v[this->colIndex[i][j]];
543 res[i] = var;
544 }
545}
546
547
548// transpose-vector multiplication
549template <typename T>
550TVector<T> CRMatrix<T>::transMult(const Vector<T>& v) const
551{
552 // Optimzed version: 7.8.98 a. Ahland
553 TVector<T> res(0, n_cols); // Zero initialisation is very important!!
554 BCHK(n_rows != v.size(), CRMatErr, Dimension conflict, v.size(), res);
555
556 unsigned int pos;
557 for (unsigned int i=0; i<n_rows; ++i) {
558 for (unsigned int j=0; j<rowSize[i]; ++j) {
559 pos = colIndex[i][j];
560 res.setval(pos) += comp[i][j] * v(i);
561 }
562 }
563 return res;
564}
565
566
567template <typename T>
568TVector<T> CRMatrix<T>::transMult(const TVector<T>& tv) const
569{
570 Vector<T> v(tv);
571 return(this->transMult(v));
572}
573
574
575template <typename T>
576TVector<T> CRMatrix<T>::transMult(const TSVector<T>& tsv) const
577{
578 Vector<T> v(tsv);
579 return(this->transMult(v));
580}
581
582
583// matrix-scalar multiplication
584template <typename T>
585CRMatrix<T> CRMatrix<T>:: operator* (const T& z) const
586{
587 CRMatrix<T> res(*this);
588 for (unsigned int i=0; i<this->n_rows; ++i)
589 for (unsigned int j=0; j<this->rowSize[i]; ++j)
590 res.comp[i][j] = this->comp[i][j] * z;
591 return res; // JENS: avoid copy constructor soon!
592}
593
594
595// matrix-scalar multiplication
596template <typename T>
597CRMatrix<T> CRMatrix<T>::mult (const T& z) const
598{
599 CRMatrix<T> res(*this);
600 for (unsigned int i=0; i<n_rows; ++i)
601 for (unsigned int j=0; j<rowSize[i]; ++j)
602 res.comp[i][j] = z * comp[i][j];
603 return res;
604}
605
606INST(template <typename T> class CRMatrix friend CRMatrix<T> operator * (const T&, const CRMatrix<T>&);)
607template <typename T>
608inline CRMatrix<T> operator* (const T& z, const CRMatrix<T>& m)
609{ return m.mult (z); }
610
611
612// fast in-place matrix-scalar multiplication
613template <typename T>
614CRMatrix<T>& CRMatrix<T>::operator*= (const T& z)
615{
616 for (unsigned int i=0; i<n_rows; ++i)
617 for (unsigned int j=0; j<rowSize[i]; ++j)
618 comp[i][j] *= z;
619 return *this;
620}
621
622// matrix-scalar division
623template <typename T>
624CRMatrix<T> CRMatrix<T>::operator/ (const T& z) const
625{
626 CRMatrix<T> res(*this);
627 for (unsigned int i=0; i<n_rows; ++i)
628 for (unsigned int j=0; j<this->rowSize[i]; ++j)
629 res.comp[i][j] = this->comp[i][j] / z;
630 return res; // JENS: avoid copy constructor soon!
631}
632
633
634// fast in-place matrix-scalar division
635template <typename T>
636CRMatrix<T>& CRMatrix<T>::operator/= (const T& z)
637{
638 for (unsigned int i=0; i<n_rows; i++)
639 for (unsigned int j=0; j<rowSize[i]; j++)
640 comp[i][j] /= z;
641 return *this;
642}
643
644// allocate memory space
645// used by constructor and resize()
646template <typename T>
647void CRMatrix<T>::allocate(unsigned int rows, unsigned int columns)
648{
649 BCHKNR(rows <= 0 || columns <= 0, CRMatErr, Zero space allocating, 0);
650 n_rows = rows; dummy = 0;
651 n_cols = columns;
652 rowSize = new unsigned int[n_rows];
653 BCHKNR(rowSize==NULL, CRMatErr, Out of memory, n_rows);
654 colIndex = new unsigned int*[n_rows];
655 BCHKNR(colIndex==NULL, CRMatErr, Out of memory, n_rows);
656 comp = new T*[n_rows];
657 BCHKNR(comp==NULL, CRMatErr, Out of memory, n_rows);
658 for (unsigned int i=0; i<n_rows; i++)
659 {
660 rowSize[i] = 0;
661 colIndex[i] = NULL;
662 comp[i] = NULL;
663 }
664}
665
666// delete memory space
667// used by "=" operator and destructor
668template <typename T>
669void CRMatrix<T>::destroy()
670{
671 for (unsigned int i=0; i<n_rows; ++i) {
672 if (colIndex[i]!=NULL) delete[] colIndex[i];
673 if (comp[i]!=NULL) delete[] comp[i];
674 }
675 if (rowSize != NULL)
676 delete[] rowSize;
677 if (comp != NULL)
678 delete[] comp;
679 if (colIndex != NULL)
680 delete[] colIndex;
681}
682
683// copy to new memory space
684// used by "=" operator and copy constructor
685template <typename T>
686void CRMatrix<T>::copy(const CRMatrix<T>& m)
687{
688 n_rows = m.n_rows;
689 n_cols = m.n_cols; dummy = 0;
690 rowSize = new unsigned int[n_rows];
691 BCHKNR(rowSize==NULL, CRMatErr, Out of memory, n_rows);
692 colIndex = new unsigned int*[n_rows];
693 BCHKNR(colIndex==NULL, CRMatErr, Out of memory, n_rows);
694 comp = new T*[n_rows];
695 BCHKNR(comp==NULL, CRMatErr, Out of memory, n_rows);
696 for (unsigned int i=0; i<n_rows; ++i) {
697 if ((rowSize[i]=m.rowSize[i]) != 0) {
698 colIndex[i] = new unsigned int[rowSize[i]];
699 BCHKNR(colIndex[i]==NULL, CRMatErr, Out of memory, rowSize[i]);
700 comp[i] = new T[rowSize[i]];
701 BCHKNR(comp[i]==NULL, CRMatErr, Out of memory, rowSize[i]);
702 for (unsigned int j=0; j<rowSize[i]; ++j) {
703 colIndex[i][j] = m.colIndex[i][j];
704 comp[i][j] = m.comp[i][j];
705 }
706 } else {
707 colIndex[i] = NULL;
708 comp[i] = NULL;
709 }
710 }
711}
712
713
714
715template <typename T>
716inline CRMatrix<T>& CRMatrix<T>::fill (const T& val)
717{
718 for (unsigned int i=0; i<n_rows; ++i)
719 for (unsigned int j=0; j<rowSize[i]; ++j)
720 comp[i][j] = val;
721 return *this;
722}
723
724template <typename T>
725inline CRMatrix<T>& CRMatrix<T>::clear ()
726{ return fill ((T)0); }
727
728template <typename T>
729inline CRMatrix<T>& CRMatrix<T>::setunit (const T& val)
730{
731 BCHK (n_rows != n_cols, CRMatErr, setunit is meaningless for non-square matrices, n_cols, *this);
732 clear ();
733 for (unsigned int i = 0; i < n_cols; ++i)
734 setval (val, i, i);
735 return *this;
736}
737
738template <typename T>
739CRMatrix<T>& CRMatrix<T>::swap(CRMatrix<T>& m)
740{
741 SWAP(n_rows, m.n_rows); SWAP(n_cols, m.n_cols);
742 SWAP(rowSize, m.rowSize); SWAP(colIndex, m.colIndex);
743 SWAP(comp, m.comp);
744 return *this;
745}
746
747template <typename T>
748CRMatrix<T> CRMatrix<T>::transposed_copy() const
749{
750 CRMatrix<T> m(columns(), rows());
751 for (unsigned r = 0; r < rows(); ++r)
752 for (unsigned c = 0; c < columns(); ++c)
753 m(c,r) = (*this)(r,c);
754 return m;
755}
756
757template <typename T>
758inline CRMatrix<T>& CRMatrix<T>::transpose()
759{
760 CRMatrix<T> tr(this->transposed_copy());
761 return swap(tr);
762}
763
764INST(template <typename T> class CRMatrix friend CRMatrix<T> transpose(const CRMatrix<T>& crm);)
765template <typename T>
766inline CRMatrix<T> transpose(const CRMatrix<T>& crm)
767{
768 return crm.transposed_copy();
769}
770
771
772template <typename T>
773STD__ ostream& operator<< (STD__ ostream& stream, const CRMatrix<T>& m)
774{
775 for (unsigned int i = 0; i < m.rows(); ++i) {
776 stream << "|\t";
777 for (unsigned int j = 0; j < m.columns(); ++j)
778 stream << m(i,j) << "\t";
779 stream << "|\n";
780 }
781 return stream << STD__ flush;
782}
783
784NAMESPACE_END
785
786#endif /* TBCI_CRMATRIX_H */
index
long int Vector< T > & index
Definition LM_fit.h:69
i
int i
Definition LM_fit.h:71
NULL
#define NULL
Definition basics.h:250
STD__
#define STD__
Definition basics.h:338
BCHK
#define BCHK(cond, exc, txt, ind, rtval)
Definition basics.h:575
MIN_ALIGN2
#define MIN_ALIGN2
Definition basics.h:424
NAMESPACE_END
#define NAMESPACE_END
Definition basics.h:323
INST
#define INST(x)
Definition basics.h:238
NAMESPACE_TBCI
#define NAMESPACE_TBCI
Definition basics.h:317
BCHKNR
#define BCHKNR(cond, exc, txt, ind)
Definition basics.h:586
FGD
#define FGD
Definition basics.h:144
T
#define T
Definition bdmatlib.cc:20
CRMatErr::CRMatErr
CRMatErr()
Definition crmatrix.h:38
CRMatErr::CRMatErr
CRMatErr(const CRMatErr &ce)
Definition crmatrix.h:42
CRMatErr::~CRMatErr
virtual ~CRMatErr()
Definition crmatrix.h:44
CRMatErr::CRMatErr
CRMatErr(const char *t, const long i=0)
Definition crmatrix.h:40
CRMatrix
C++ class for sparse matrices using compressed row storage.
Definition crmatrix.h:64
CRMatrix::transpose
CRMatrix< T > & transpose()
Definition crmatrix.h:758
CRMatrix::setval
CRMatrix< T > & setval(const T &z, unsigned int row, unsigned int column)
element access (write)
Definition crmatrix.h:215
CRMatrix::size
unsigned int size() const
Definition crmatrix.h:193
CRMatrix::set_row
void set_row(const Vector< T > &v, unsigned int row)
Definition crmatrix.h:323
CRMatrix::columns
unsigned int columns() const
Definition crmatrix.h:119
CRMatrix::setunit
CRMatrix< T > & setunit(const T &=(T) 1)
set CSCMatrix to val times the unit matrix
Definition crmatrix.h:729
CRMatrix::TALIGN
T aligned_value_type TALIGN(MIN_ALIGN2)
Definition crmatrix.h:81
CRMatrix::clear
CRMatrix< T > & clear()
set all elements defined to zero
Definition crmatrix.h:725
CRMatrix::operator==
bool operator==(const CRMatrix< T > &m) const
Definition crmatrix.h:452
CRMatrix::copy
void copy(const CRMatrix< T > &m)
Definition crmatrix.h:686
CRMatrix::fill
CRMatrix< T > & fill(const T &)
set all defined element to a val
Definition crmatrix.h:716
CRMatrix::get_row
TVector< T > get_row(unsigned int row) const
Definition crmatrix.h:312
CRMatrix::operator=
CRMatrix< T > & operator=(const CRMatrix< T > &m)
Definition crmatrix.h:409
CRMatrix::MatVecMult
void MatVecMult(Vector< T > &res, const Vector< T > &v) const
Definition crmatrix.h:519
CRMatrix::resize
void resize(unsigned int newRows, unsigned int newColumns)
Definition crmatrix.h:205
CRMatrix::operator*=
CRMatrix< T > & operator*=(const T &z)
Definition crmatrix.h:614
CRMatrix::colIndex
unsigned int ** colIndex
Definition crmatrix.h:68
CRMatrix::value_type
T value_type
Definition crmatrix.h:79
CRMatrix::operator-
CRMatrix< T > operator-() const
Definition crmatrix.h:466
CRMatrix::destroy
void destroy()
Definition crmatrix.h:669
CRMatrix::n_rows
unsigned int n_rows
Definition crmatrix.h:66
CRMatrix::comp
T ** comp
Definition crmatrix.h:69
CRMatrix::transMult
TVector< T > transMult(const Vector< T > &v) const
Definition crmatrix.h:550
CRMatrix::operator/=
CRMatrix< T > & operator/=(const T &z)
Definition crmatrix.h:636
CRMatrix::swap
CRMatrix< T > & swap(CRMatrix< T > &)
Definition crmatrix.h:739
CRMatrix::allocate
void allocate(unsigned int rows, unsigned int columns)
Definition crmatrix.h:647
CRMatrix::~CRMatrix
~CRMatrix()
Definition crmatrix.h:87
CRMatrix::rows
unsigned int rows() const
Definition crmatrix.h:118
CRMatrix::operator!=
bool operator!=(const CRMatrix< T > &m) const
Definition crmatrix.h:128
CRMatrix::CRMatrix
CRMatrix(const CRMatrix< T > &m)
Definition crmatrix.h:85
CRMatrix::mat_info
static const char * mat_info()
Definition crmatrix.h:90
CRMatrix::CRMatrix
CRMatrix(unsigned int rows, unsigned int columns)
Definition crmatrix.h:84
CRMatrix::operator()
const T & operator()(unsigned int row, unsigned int column) const
element access (read)
Definition crmatrix.h:181
CRMatrix::CRMatrix
CRMatrix()
Definition crmatrix.h:83
CRMatrix::mult
CRMatrix< T > mult(const T &z) const
Definition crmatrix.h:597
CRMatrix::rowSize
unsigned int * rowSize
Definition crmatrix.h:67
CRMatrix::element_type
T element_type
Definition crmatrix.h:80
CRMatrix::operator*
TVector< T > operator*(const Vector< T > &v) const
Definition crmatrix.h:479
CRMatrix::get
const T & get(unsigned int row, unsigned int column) const
Definition crmatrix.h:101
CRMatrix::transposed_copy
CRMatrix< T > transposed_copy() const
Inefficient! Use transMult if possible.
Definition crmatrix.h:748
CRMatrix::n_cols
unsigned int n_cols
Definition crmatrix.h:66
CRMatrix::operator/
CRMatrix< T > operator/(const T &z) const
Definition crmatrix.h:624
Matrix_Sig::Matrix_Sig
Matrix_Sig()
Definition matrix_sig.h:49
NumErr::NumErr
NumErr()
Definition except.h:65
TSVector::destroy
void destroy() const
Definition vector.h:1172
TSVector::get
T get(const unsigned long i) const HOT
Definition vector.h:1072
TSVector::size
unsigned long size() const
Definition vector.h:1120
TVector
Temporary Base Class Idiom: Class TVector is used for temporary variables.
Definition vector.h:73
TVector::size
unsigned long size() const
Definition vector.h:104
TVector::set
T & set(const T &val, const unsigned long i) const
Definition vector.h:194
MatVecMult
void MatVecMult(Vector< T > &res, const CRMatrix< T > &m, const Vector< T > &v)
Definition crmatrix.h:532
transpose
CRMatrix< T > transpose(const CRMatrix< T > &crm)
Definition crmatrix.h:766
operator<<
STD__ ostream & operator<<(STD__ ostream &stream, const CRMatrix< T > &m)
Definition crmatrix.h:773
c
return c
Definition f_matrix.h:760
res
const unsigned TMatrix< T > * res
Definition matrix_kernels.h:26
Generated by doxygen 1.16.1