TBCI Numerical high perf. C++ Library  2.8.0
f_matrix.h
Go to the documentation of this file.
1 
6 //-------------------------------------------------------------
7 // A M Bilgic, 1/97
8 // last update 97/10/04 AMB
9 // 98/05/07 AMB
10 //
11 // $Id: f_matrix.h,v 1.23.2.59 2019/05/29 06:22:59 garloff Exp $
12 //-------------------------------------------------------------
13 
14 #ifndef TBCI_F_MATRIX_H
15 #define TBCI_F_MATRIX_H
16 
17 #include "tbci/vector.h"
18 #include "tbci/matrix_sig.h"
19 #include "tbci/matrix.h"
20 #include "tbci/cscmatrix.h"
21 
22 // Avoid -fguiding-decls
23 #if !defined(NO_GD) && !defined(AUTO_DECL)
24 # include "tbci/f_matrix_gd.h"
25 #endif
26 
27 NAMESPACE_TBCI
28 
29 #ifdef HAVE_GCC295_FRIEND_BUG
30 # define _VEC getvec()
31 # define _ENDVEC getendvec()
32 # define _DIM size()
33 # define _ROW rows()
34 # define _COL columns()
35 # define _FAC getfac()
36 #else
37 # define _VEC vec
38 # define _ENDVEC endvec
39 # define _DIM dim
40 # define _ROW row
41 # define _COL col
42 # define _FAC fac
43 #endif
44 
45 
47 
48 template <typename T> class CSCMatrix;
49 template <typename T> class F_Matrix;
50 template <typename T> class F_TMatrix;
51 template <typename T> class F_TSMatrix;
52 template <typename T> class BdMatrix;
53 template <typename T> class Tensor;
54 
55 template <typename T> class Matrix;
56 template <typename T> class TMatrix;
57 
58 #ifdef PRAGMA_I
59 # pragma interface "f_matrix.h"
60 #endif
61 
62 //-----------------------------------------------------
67 //-----------------------------------------------------
68 
69 template <typename T>
70 class F_TMatrix : public Matrix_Sig<T>
71 {
72  protected:
73  typedef T mat_t /*TALIGN(sizeof(T))*/;
74  unsigned long dim;
75  unsigned int col, row;
77  T ** mat;
78  T *vec, *endvec;
79 
80  int set_ptrs();
81 
82  public:
83  typedef T value_type;
84  typedef T element_type;
85  typedef T aligned_value_type TALIGN(MIN_ALIGN2);
86 
87  friend class F_Matrix<T>;
88  friend class F_TSMatrix<T>;
89  friend class Tensor<T>;
90  //friend class CSCMatrix<T>;
91 
92 #ifdef HAVE_GCC295_FRIEND_BUG
93  // Provide access to the internals to workaround gcc-2.95 hostility to friends
94  T* const getvec() const { return vec; }
95  T* const getendvec() const { return endvec; }
96 #endif
97 
98  //protected:
99  public:
100  explicit F_TMatrix (const unsigned = 0);
101  F_TMatrix (const unsigned, const unsigned);
102  F_TMatrix (const T&, const unsigned, const unsigned);
103  F_TMatrix (const Vector<T>&, const enum rowcolvec = colvec);
104  //aliasing
105  F_TMatrix (const F_TMatrix<T>&);
106  F_TMatrix (F_TSMatrix<T>&);
107  //copying
108  F_TMatrix (const F_Matrix<T>&);
109  // real destructor
110  void destroy ();
111 
112  //public:
113  // empty destructor: F_TMatrix will never be destroyed but explicitly !!!
114  ~F_TMatrix () {}
115 
116  // conversion
117  operator TMatrix<T>();
118  // allow instantiation (Matrix_Sig)
119  /*virtual*/ static const char* mat_info () { return ("F_TMatrix"); }
120 
121  //public:
122  // assignment
123  F_TMatrix<T>& operator = (const F_Matrix<T>&);
124  F_TMatrix<T>& operator = (const F_TMatrix<T>&);
125  F_TMatrix<T>& operator = (F_TSMatrix<T>);
126  F_TMatrix<T>& operator = (const T&);
127  F_TMatrix<T>& alias (const F_Matrix<T>& m)
128  { destroy (); dim = m.dim; col = m.col; row = m.row;
129  mat = m.mat; vec = m.vec; endvec = m.endvec; return *this; }
130 
131  F_TMatrix<T>& operator += (F_TMatrix<T>); // arg will be destroyed
132  F_TMatrix<T>& operator += (const F_Matrix<T>&); // arg wont be " "
133  F_TMatrix<T>& operator += (const T&);
134  F_TMatrix<T>& operator -= (F_TMatrix<T>);
135  F_TMatrix<T>& operator -= (const F_Matrix<T>&);
136  F_TMatrix<T>& operator -= (const T&);
137 
138  F_TSMatrix<T> operator *= (const T&);
139  F_TSMatrix<T> operator /= (const T&);
140 
141  // operations
142  F_TMatrix<T>& operator - ();
143  F_TMatrix<T>& herm (); //returns the the hermitian adjungated (Tronsposed and conjugate complex)
144  F_TMatrix<T>& conj ();
145  F_TMatrix<T>& trans();
146 
147  F_TMatrix<T>& operator + (F_TMatrix<T>); // arg will be destroyed
148  F_TMatrix<T>& operator + (F_TSMatrix<T>); // arg will be destroyed
149  F_TMatrix<T>& operator + (const F_Matrix<T>&); // arg wont be destroyed
150  F_TMatrix<T>& operator + (const T&);
151  F_TMatrix<T>& operator - (F_TMatrix<T>);
152  F_TMatrix<T>& operator - (F_TSMatrix<T>);
153  F_TMatrix<T>& operator - (const F_Matrix<T>&);
154  F_TMatrix<T>& operator - (const T&);
155  F_TSMatrix<T> operator * (const T&);
156  F_TSMatrix<T> operator / (const T&);
157 
158  F_TMatrix<T> operator * (const F_Matrix<T>&);
159  F_TMatrix<T> operator * (F_TMatrix<T>);
160  F_TMatrix<T> operator * (F_TSMatrix<T>);
161  // TODO: Implement transMult
162 
163  TVector<T> operator * (const Vector<T>& v)
164  { F_Matrix<T> m (*this); return (m * v); }
165  TVector<T> operator * (TVector <T>& tv)
166  { F_Matrix<T> m (*this); return (m * tv); }
167 
168  // friend operations
169 #ifndef HAVE_GCC295_FRIEND_BUG
170  friend NOINST F_TMatrix<T> FRIEND_TBCI__ operator + FGD (const T&, F_TMatrix<T>);
171  friend NOINST F_TMatrix<T> FRIEND_TBCI__ operator + FGD (const T&, const F_Matrix<T>&);
172  friend NOINST F_TMatrix<T> FRIEND_TBCI__ operator - FGD (const T&, F_TMatrix<T>);
173  friend NOINST F_TMatrix<T> FRIEND_TBCI__ operator - FGD (const T&, const F_Matrix<T>&);
174  friend NOINST F_TSMatrix<T> FRIEND_TBCI__ operator * FGD (const T&, F_TMatrix<T>);
175  friend NOINST F_TSMatrix<T> FRIEND_TBCI__ operator * FGD (const T&, const F_Matrix<T>&);
176 #endif
177 
178  // Access
179  T& operator () (const unsigned int, const unsigned int) const;
180  //T operator () (const unsigned, const unsigned) const;
181  T* get_fortran_matrix () const { return vec; }
182  const T& get (const unsigned r, const unsigned c) const { return mat[c][r]; }
183 
184  bool operator == (const F_Matrix<T>& m);
185  bool operator != (const F_Matrix<T>& m)
186  { return !(*this == m); }
187 
188  bool operator == (F_TMatrix<T> tm)
189  { F_Matrix<T> m (tm); return (*this == m); }
190  bool operator != (F_TMatrix<T> tm)
191  { return !(*this == tm); }
192 
193  bool operator == (F_TSMatrix<T>);
194  bool operator != (F_TSMatrix<T> ts)
195  { return !(*this == ts); }
196 
197  //public:
198  unsigned int columns () const { return col; }
199  unsigned int rows () const { return row; }
200  unsigned long size () const { return dim; }
201  //TVector<T> operator () (const unsigned int) const;
202  TVector<T> get_row (const unsigned int) const;
203  TVector<T> get_col (const unsigned int) const;
204  void set_row (const Vector<T>&, const unsigned int);
205  void set_col (const Vector<T>&, const unsigned int);
206  void setval(const T val, const unsigned int r, const unsigned int c)
207  { this->operator() (r, c) = val; }
208  T& setval(const unsigned int r, const unsigned int c)
209  { return this->operator() (r, c); }
210 
211  F_TMatrix<T>& resize (const unsigned int, const unsigned int);
212  F_TMatrix<T>& resize (const unsigned int d) { return this->resize (d, d); }
213  F_TMatrix<T>& resize (const T&, const unsigned int, const unsigned int);
214  F_TMatrix<T>& fill (const T& = 0);
215  F_TMatrix<T>& clear ();
216  F_TMatrix<T>& setunit(const T& = 1);
217 
218  F_TMatrix<T>& swap (F_TMatrix<T>&);
220  F_TMatrix<T> transposed_copy() const;
221  F_TMatrix<T>& transpose();
222 
223  //friend double fabs FGD (const F_Matrix<T>&);
224  //friend double fabs FGD (const F_TMatrix<T>&);
225  double fabs () const;
226 
227  T trace () const;
228 
229 };
230 
231 
232 template <typename T>
233 inline void F_TMatrix<T>::destroy ()
234 { /* Skip if (dim) ? */
235  if (dim) {
236  if (vec) TBCIDELETE(T , vec, dim);
237  if (mat) TBCIDELETE(T*, mat, col);
238  dim = 0;
239  }
240 }
241 
242 template <typename T>
243 F_TMatrix<T>::operator TMatrix<T>()
244 {
245  TMatrix<T> tm(row, col);
246  for (unsigned int c = 0; c < col; ++c)
247  for (unsigned int r = 0; r < row; ++r)
248  tm.setval(this->get(r, c), r, c);
249  destroy();
250  return tm;
251 }
252 
253 template <typename T>
254 inline int F_TMatrix<T>::set_ptrs()
255 {
256  if (mat && vec) {
257  T* ptr = vec;
258  for (unsigned int i=0; i<col; i++) {
259  mat[i] = ptr;
260  ptr += row;
261  }
262  endvec = ptr;
263  return 0;
264  } else {
265  if (vec) { TBCIDELETE(T, vec, dim); vec = (T*) 0; }
266  if (mat) { TBCIDELETE(T*,mat, col); mat = (T**)0; }
267  dim = 0; col = 0; row = 0;
268  endvec = (T*)0;
269  return -1;
270  }
271 }
272 
273 
274 template <typename T>
275 inline F_TMatrix<T>::F_TMatrix (const unsigned d)
276  : dim((unsigned long)d*d), col(d), row(d),
277  mat(d>0?NEW(T*,d):0), vec(d>0?NEW(T,d*d):0)
278 {
279  set_ptrs();
280 }
281 
282 
283 template <typename T>
284 inline F_TMatrix<T>::F_TMatrix (const unsigned r, const unsigned c)
285  : dim(((unsigned long)r)*c), col(c), row(r),
286  mat(dim>0?NEW(T*,c):0), vec(dim>0?NEW(T,dim):0)
287 {
288  set_ptrs();
289 }
290 
291 template <typename T>
292 inline F_TMatrix<T>::F_TMatrix (const T& val, const unsigned r, const unsigned c)
293  : dim(((unsigned long)r)*c), col(c), row(r),
294  mat(dim>0?NEW(T*,c):0), vec(dim>0?NEW(T,dim):0)
295 {
296  if (LIKELY(!set_ptrs()))
297  TBCIFILL (vec, val, T, dim);
298 }
299 
300 
301 template <typename T>
302 inline F_TMatrix<T>::F_TMatrix (const Vector<T>& v, const enum rowcolvec r)
303  : dim(v.dim), col(r==rowvec?1:v.dim), row(r==rowvec?v.dim:1),
304  mat(v.dim?NEW(T*,col):0), vec(v.dim>0?NEW(T,v.dim):0)
305 {
306  if (LIKELY(!set_ptrs()))
307  TBCICOPY (vec, v.vec, T, dim);
308 }
309 
310 // data is duplicated
311 template <typename T>
312 inline F_TMatrix<T>::F_TMatrix (const F_Matrix<T>& m)
313  : dim(m.dim), col(m.col), row(m.row),
314  mat(m.dim>0?NEW(T*,m.col):0),
315  vec(m.dim>0?NEW(T,m.dim):0)
316 {
317  if (LIKELY(!set_ptrs()))
318  TBCICOPY (vec, m.vec, T, dim);
319 }
320 
321 // copy constructor only copies pointers !
322 template <typename T>
323 inline F_TMatrix<T>::F_TMatrix (const F_TMatrix<T>& tm)
324  : dim (tm.dim), col (tm.col), row (tm.row),
325  mat (tm.mat), vec (tm.vec), endvec (tm.endvec)
326 {}
327 
328 template <typename T>
329 inline F_TMatrix<T>::F_TMatrix (F_TSMatrix<T>& ts)
330  : dim (ts.dim), col (ts.col), row (ts.row)
331 {
332  ts.eval (); vec = ts.vec; mat = ts.mat;
333  endvec = ts.endvec; ts.mut = false;
334 }
335 
336 
337 template <typename T>
338 inline T& F_TMatrix<T>::operator () (const unsigned int r, const unsigned int c) const
339 {
340  BCHK(r>=row, MatErr, Illegal row index, r, mat[0][0]);
341  BCHK(c>=col, MatErr, illegal col index, c, mat[0][0]);
342  return mat[c][r];
343 }
344 
345 /*
346 template <typename T>
347 inline T F_TMatrix<T>::operator () (const unsigned i, const unsigned j) const
348 {
349  BCHK(i<0 || i>=row, MatErr, Illegal row index, i, mat[0][0]);
350  BCHK(j<0 || j>=col, MatErr, illegal col index, j, mat[0][0]);
351  return mat[j][i];
352 }
353  */
354 
355 
356 //copy
357 template <typename T>
358 inline F_TMatrix<T>& F_TMatrix<T>::operator = (const F_Matrix<T>& a)
359 {
360  BCHK(dim != a.dim, MatErr, Different sizes in assignment, a.dim, *this );
361  TBCICOPY (vec, a.vec, T, dim);
362  return *this;
363 }
364 
365 //alias (do we need it?)
366 template <typename T>
367 inline F_TMatrix<T>& F_TMatrix<T>::operator = (const F_TMatrix<T>& a)
368 {
369  BCHK(row != a.row || col != a.col, MatErr, Different sizes in assignment, a.dim, *this );
370  if (vec != a.vec) /* Is == possible ? */
371  {
372  this->destroy();
373  dim = a.dim; col = a.col; row = a.row;
374  mat = a.mat; vec = a.vec; endvec = a.endvec;
375  }
376  return *this;
377 }
378 
379 template <typename T>
380 inline F_TMatrix<T>& F_TMatrix<T>::operator = (F_TSMatrix<T> a)
381 {
382  BCHK(col != a.col || row != a.row, MatErr, Different sizes in assignment, a.dim, *this );
383  a.eval (this);
384  return *this;
385 }
386 
387 //fill
388 template <typename T>
389 inline F_TMatrix<T>& F_TMatrix<T>::operator = (const T& val)
390 {
391  return this->fill (val);
392 }
393 
394 
395 template <typename T>
396 inline F_TMatrix<T>& F_TMatrix<T>::setunit (const T& fac)
397 {
398  TBCIFILL (vec, 0, T, dim);
399  BCHK(col!=row, MatErr, setunit makes only sense on quadratic matrices, col, *this);
400  const unsigned iend = MIN(row,col);
401  for (register unsigned i = 0; i < iend; i++) mat[i][i] = fac;
402  return *this;
403 }
404 
405 template <typename T>
406 inline F_TMatrix<T>& F_TMatrix<T>::clear ()
407 {
408  if (dim) TBCIFILL (vec, 0, T, dim);
409  return *this;
410 }
411 
412 #if 0
413 template <typename T>
414 inline TVector<T> F_TMatrix<T>::operator () (const unsigned int i) const
415 {
416  TVector<T> v(col);
417  BCHK(i<0 || i>= row, MatErr, Illegal row index, i, v);
418  TBCICOPY (v.vec, mat[i], T, col);
419  return v;
420 }
421 #endif
422 
423 template <typename T>
424 inline void F_TMatrix<T>::set_col (const Vector<T>& v, const unsigned int c)
425 {
426  BCHKNR(c >= col, MatErr, Illegal col index in set_col, c);
427  BCHKNR(v.dim != (unsigned long)row, MatErr, Different sizes in set_row, v.dim);
428  TBCICOPY(mat[c], v.vec, T, v.dim);
429 }
430 
431 template <typename T>
432 inline TVector<T> F_TMatrix<T>::get_col (const unsigned int c) const
433 {
434  TVector<T> v(row);
435  //BCHK(c >= col, MatErr, Illegal col index in get_col, c, v);
436  TBCICOPY (v.vec, mat[c], T, row);
437  return v;
438 }
439 
440 template <typename T>
441 inline void F_TMatrix<T>::set_row (const Vector<T>& v, const unsigned int r)
442 {
443  BCHKNR(r >= row, MatErr, Illegal row index in set_row, r);
444  for (register unsigned int j=0; j<col; j++)
445  mat[j][r] = v.vec[j];
446 }
447 
448 template <typename T>
449 inline TVector<T> F_TMatrix<T>::get_row (const unsigned int r) const
450 {
451  TVector<T> v(col);
452  BCHK(r >= row, MatErr, Illegal row index in get_row, r, v);
453  for (register unsigned int j=0; j<col; j++)
454  v.vec[j] = mat[j][r];
455  return v;
456 }
457 
458 
459 #ifndef LAPACK_INLINE
460 # define LAPACK_INLINE
461 #endif
462 template <typename T>
463 LAPACK_INLINE F_TMatrix<T>& F_TMatrix<T>::resize (const unsigned int r, const unsigned int c)
464 {
465  if (r == row && c == col)
466  return *this;
467  if (mat)
468  TBCIDELETE(T*,mat,col);
469  T* vv = vec;
470  unsigned long dd =dim;
471  dim = ((unsigned long)r)*c;
472  col = c; row = r;
473  if (dim) {
474  mat = NEW(T*,c); vec = NEW(T,dim);
475  } else {
476  mat = (T**)0; vec = (T*)0; endvec = (T*)0;
477  }
478  if (!set_ptrs())
479  TBCICOPY(vec, vv, T, MIN(dim,dd));
480  if (dd && vv)
481  TBCIDELETE(T, vv, dd);
482  return *this;
483 }
484 
485 template <typename T>
486 F_TMatrix<T>& F_TMatrix<T>::resize (const T& v, const unsigned int r, const unsigned int c)
487 {
488  if (r == row && c == col) {
489  TBCIFILL(vec, v, T, dim);
490  return *this;
491  }
492  if (mat) {
493  TBCIDELETE(T*,mat,col);
494  TBCIDELETE(T,vec, dim);
495  }
496  dim = r*c;
497  row = r; col = c;
498  if (dim) {
499  mat = NEW(T*,c); vec = NEW(T,dim);
500  } else {
501  mat = (T**)0; vec = (T*)0; endvec = (T*)0;
502  }
503  if (LIKELY(!set_ptrs()))
504  TBCIFILL (vec, v, T, dim);
505  return *this;
506 }
507 
508 
509 template <typename T>
510 inline F_TMatrix<T>& F_TMatrix<T>::fill (const T& val)
511 {
512  TBCIFILL(vec, val, T, dim);
513  return *this;
514 }
515 
516 template <typename T>
517 inline T F_TMatrix<T>::trace () const
518 {
519  ALIGN3(T tr, vec[0], MIN_ALIGN2);
520  BCHK (row != col, MatErr, Trace over non quadratic matrix, 0, tr=0);
521  for (register unsigned i = 1; i < row; i++)
522  tr += mat[i][i];
523  return tr;
524 }
525 
526 /************* Arithmetics *****************/
527 
528 
529 template <typename T>
530 F_TMatrix<T>& F_TMatrix<T>::trans()
531 {
532  if (rows() == columns())
533  {
534  for(unsigned int i=0; i<rows(); i++)
535  for(unsigned int j=i+1; j<columns(); j++)
536  SWAP(mat[i][j], mat[j][i]);
537  } else
538  {
539  STD__ cout << "Not implemented!" << STD__ endl;
540  }
541  return (*this);
542 }
543 
544 
545 template <typename T>
546 F_TMatrix<T>& F_TMatrix<T>::herm ()
547 {
548  STD__ cout << "WRONG: F_TMatrix<T>::herm()" << STD__ endl;
549  if (rows() == columns())
550  {
551  for(unsigned int i=0; i<rows(); i++)
552  {
553 // mat[i][i]=(mat[i][i]).conj();
554  for(unsigned int j=i+1; j<columns(); j++)
555  {
556 // mat[i][j]=(mat[i][j]).conj();
557 // mat[j][i]=(mat[j][i]).conj();
558  SWAP(mat[i][j], mat[j][i]);
559  }
560  }
561  } else
562  {
563  STD__ cout << "Not implemented!" << STD__ endl;
564  }
565  return (*this);
566 }
567 
568 
569 template <typename T>
570 F_TMatrix<T>& F_TMatrix<T>::conj ()
571 {
572  STD__ cout << "Wrong: F_TMatrix<T>::conf()" << STD__ endl;
573 // for (register T* ptr=vec; ptr<endvec; ptr++) *ptr = (*ptr).conj();
574  return *this;
575 }
576 
577 // unary -
578 template <typename T>
579 inline F_TMatrix<T>& F_TMatrix<T>::operator - ()
580 {
581  for (register T* ptr=vec; ptr<endvec; ptr++) *ptr = -(*ptr);
582  return *this;
583 }
584 
585 template <typename T>
586 bool F_TMatrix<T>::operator == (const F_Matrix<T>& m)
587 {
588  F_Matrix<T> th (*this);
589  if (row != m.row || col != m.col) return false;
590  if (vec == m.vec || dim == 0) return true;
591  if (TBCICOMP (vec, m.vec, T, dim)) return false;
592  /*else*/ return true;
593 }
594 
595 template <typename T>
596 bool F_TMatrix<T>::operator == (F_TSMatrix<T> ts)
597 {
598  F_Matrix<T> th (*this);
599  if (row != ts.row || col != ts.col) return (ts.destroy(), false);
600  if (dim == 0) { /*ts.destroy();*/ return true; }
601  if (ts.fac == (T)1)
602  {
603  if (vec == ts.vec) return true;
604  if (TBCICOMP (vec, ts.vec, T, dim)) { ts.destroy (); return false; }
605  }
606  else
607  {
608  if (vec == ts.vec) return false;
609  for (register T *p1 = vec, *p2 = ts.vec; p1 < endvec; p1++, p2++)
610  if (*p1 != *p2 * ts.fac) { ts.destroy (); return false; }
611  }
612  ts.destroy (); return true;
613 }
614 
615 
616 #define F_TMFORALL_M(op) \
617 template <typename T> \
618 inline F_TMatrix<T>& F_TMatrix<T>::operator op (const F_Matrix<T>& a) \
619 { \
620  BCHK(row != a.row, MatErr, number of rows differ in op, a.row, *this ); \
621  BCHK(col != a.col, MatErr, number of cols differ in op, a.col, *this ); \
622  for (register T *p1 = vec, *p2 = a.vec; p1 < endvec; p1++, p2++) *p1 op *p2; \
623  return *this; \
624 }
625 
626 F_TMFORALL_M(+=)
627 F_TMFORALL_M(-=)
628 
629 
630 // Changing a F_TMatrix to a F_Matrix destroys it ...
631 #define F_TMFORALL_TM(op) \
632 template <typename T> \
633 inline F_TMatrix<T>& F_TMatrix<T>::operator op (F_TMatrix<T> a)\
634 { \
635  return this->operator op (F_Matrix<T>(a)); \
636 }
637 
638 F_TMFORALL_TM(+=)
639 F_TMFORALL_TM(-=)
640 
641 
642 #define F_TMFORALL_TS(op) \
643 template <typename T> \
644 inline F_TMatrix<T>& F_TMatrix<T>::operator op (F_TSMatrix<T> ts) \
645 { \
646  BCHK(row != ts.row, MatErr, number of rows differ in op, ts.row, (ts.destroy(),*this) ); \
647  BCHK(col != ts.col, MatErr, number of cols differ in op, ts.col, (ts.destroy(),*this) ); \
648  for (register T *p1 = vec, *p2 = ts.vec; p1 < endvec; p1++, p2++) \
649  *p1 op##= *p2 * ts.fac; \
650  ts.destroy (); return *this; \
651 }
652 
653 F_TMFORALL_TS(+)
654 F_TMFORALL_TS(-)
655 
656 
657 //Change these to operate on diag only ?
658 #define F_TMFORALL_T(op) \
659 template <typename T> \
660 inline F_TMatrix<T>& F_TMatrix<T>::operator op (const T& a) \
661 { \
662  for (register T* ptr = vec; ptr < endvec; ptr++) *ptr op a; \
663  return *this; \
664 }
665 
666 F_TMFORALL_T(+=)
667 F_TMFORALL_T(-=)
668 //F_TMFORALL_T(*=)
669 //F_TMFORALL_T(/=)
670 
671 
672 template <typename T>
673 inline F_TSMatrix<T> F_TMatrix<T>::operator *= (const T& a)
674 { return F_TSMatrix<T> (*this, a); }
675 
676 template <typename T>
677 inline F_TSMatrix<T> F_TMatrix<T>::operator /= (const T& a)
678 {
679  BCHK(a==(T)0, MatErr, F_Matrix division by zero, 0, F_TSMatrix<T> (*this));
680  return F_TSMatrix<T> (*this, (T)1/a);
681 }
682 
683 
684 // temporaries can be changed, because we can't refer them
685 #define F_STDDEF_TMM(op) \
686 template <typename T> \
687 inline F_TMatrix<T>& F_TMatrix<T>::operator op (const F_Matrix<T>& a) \
688 { return this->operator op##= (a); }
689 
690 F_STDDEF_TMM(+)
691 F_STDDEF_TMM(-)
692 
693 
694 #define F_STDDEF_TMTM(op) \
695 template <typename T> \
696 inline F_TMatrix<T>& F_TMatrix<T>::operator op (F_TMatrix<T> a) \
697 { return this->operator op##= (F_Matrix<T>(a)); }
698 
699 F_STDDEF_TMTM(+)
700 F_STDDEF_TMTM(-)
701 
702 
703 #define F_STDDEF_TMT(op) \
704 template <typename T> \
705 inline F_TMatrix<T>& F_TMatrix<T>::operator op (const T& a) \
706 { return this->operator op##= (a); }
707 
708 F_STDDEF_TMT(+)
709 F_STDDEF_TMT(-)
710 //STDDEF_TMT(*)
711 
712 template <typename T>
713 inline F_TSMatrix<T> F_TMatrix<T>::operator * (const T& a)
714 { return F_TSMatrix<T> (*this, a); }
715 
716 template <typename T>
717 inline F_TSMatrix<T> F_TMatrix<T>::operator / (const T& b)
718 {
719  BCHK((b == (T)0), MatErr, F_Matrix division by 0, 0, F_TSMatrix<T> (*this));
720  return F_TSMatrix<T> (*this, (T)1/b);
721 }
722 
723 //friends
724 
725 //Change to operate on diag only ?
726 #define F_STDDEF_TTM(op) \
727 INST(template <typename T> class F_TMatrix friend F_TMatrix<T> operator op (const T&, F_TMatrix<T>);) \
728 template <typename T> \
729 inline F_TMatrix<T> operator op (const T& a, F_TMatrix<T> b) \
730 { \
731  for (register T* ptr=b._VEC; ptr<b._ENDVEC; ptr++) \
732  *ptr = a op *ptr; \
733  return b; \
734 }
735 
736 F_STDDEF_TTM(+)
737 F_STDDEF_TTM(-)
738 //STDDEF_TTM(*)
739 //STDDEF_TTM(/)
740 
741 INST(template <typename T> class F_TMatrix friend F_TSMatrix<T> operator * (const T&, F_TMatrix<T>);)
742 template <typename T>
743 inline F_TSMatrix<T> operator * (const T& a, F_TMatrix<T> b)
744 { return F_TSMatrix<T> (b, a); }
745 // Commutativity is assumed here !
746 
747 
748 //Change to operate on diag only ?
749 #define F_STDDEF_TM(op) \
750 INST(template <typename T> class F_TMatrix friend F_TMatrix<T> operator op (const T&, const F_Matrix<T>&);) \
751 template <typename T> \
752 inline F_TMatrix<T> operator op (const T& a, const F_Matrix<T>& b) \
753 { \
754  F_TMatrix<T> c (b._ROW, b._COL); \
755  for (register T *p1=c._VEC, *p2=b._VEC; p1<c._ENDVEC; p1++, p2++) \
756  *p1 = a op *p2; \
757  return c; \
758 }
759 
760 F_STDDEF_TM(+)
761 F_STDDEF_TM(-)
762 //STDDEF_TM(*)
763 //STDDEF_TM(/)
764 
765 
766 INST(template <typename T> class F_TMatrix friend F_TSMatrix<T> operator * (const T&, const F_Matrix<T>&);)
767 template <typename T>
768 inline F_TSMatrix<T> operator * (const T& a, const F_Matrix<T>& b)
769 { return F_TSMatrix<T> (b, a); }
770 // Commutativity is assumed here !
771 
772 
773 template <typename T>
774 inline double F_TMatrix<T>::fabs () const
775 { F_Matrix<T> m (*this); return m.fabs (); }
776 
777 //-------------------------------------------------------------
781 //-------------------------------------------------------------
782 template <typename T>
783 class F_TSMatrix : public Matrix_Sig<T> //: public F_TMatrix
784 {
785  protected:
786  unsigned long dim;
787  unsigned int col, row;
788  T** mat;
789  T *vec, *endvec;
790 
791  void destroy ();
792  T fac ALIGN(MIN_ALIGN); //factor
793  bool mut; //mutability: does vec point to changeable (non-refered) data ?
794 
795  void clone (bool = false, F_TMatrix<T>* = 0);
796 #ifdef HAVE_GCC295_FRIEND_BUG
797  public:
798  T* const getvec () const { return vec; }
799  T* const getendvec () const { return endvec; }
800  T& getfac () { return fac; }
801  const T& getfac () const { return fac; }
802 #endif
803  void detach (F_TMatrix<T>* = 0);
804 
805  public:
806  typedef T value_type;
807  typedef T element_type;
808  typedef T aligned_value_type TALIGN(MIN_ALIGN2);
809 
810  friend class F_TMatrix<T>;
811  friend class F_Matrix<T>;
812 
813  F_TSMatrix () : dim(0), col(0), row(0), mat((T**)0), vec((T*)0), endvec((T*)0) {}
814  // empty destructor: F_TSMatrix will never be destroyed but explicitly !!!
815  ~F_TSMatrix () {}
816 
817  F_TSMatrix (const F_TMatrix<T>& tm, const T& f = (T)1)
818  : dim(tm.dim), col(tm.col), row(tm.row),
819  mat(tm.mat), vec(tm.vec), endvec(tm.endvec), fac(f), mut(true) {}
820  F_TSMatrix (const F_Matrix<T>& m, const T& f = (T)1)
821  : dim(m.dim), col(m.col), row(m.row),
822  mat(m.mat), vec(m.vec), endvec(m.endvec), fac(f), mut(false) {}
823  F_TSMatrix (const F_TSMatrix<T>& ts)
824  : dim(ts.dim), col(ts.col), row(ts.row),
825  mat(ts.mat), vec(ts.vec), endvec(ts.endvec), fac(ts.fac), mut(ts.mut) {}
826 
827  T operator () (const unsigned r, const unsigned c)
828  { T val = mat[c][r] * fac; destroy (); return val; }
829  T get (const unsigned r, const unsigned c) const { return mat[c][r] * fac; }
830 
831  F_TSMatrix<T>& eval (F_TMatrix<T>* = 0);
832  operator TMatrix<T>();
833 
834  // allow instantiation (Matrix_Sig)
835  /*virtual*/ static const char* mat_info () { return ("F_TSMatrix"); }
836 
837  // These operations are very cheap now
838  F_TSMatrix<T>& operator = (const F_TSMatrix<T>& ts)
839  { destroy (); dim = ts.dim; col = ts.col; row = ts.row;
840  mat = ts.mat; vec = ts.vec; endvec = ts.endvec;
841  fac = ts.fac; mut = ts.mut; return *this; }
842  F_TSMatrix<T>& operator = (const F_TMatrix<T>& tm)
843  { destroy (); dim = tm.dim; col = tm.col; row = tm.row;
844  mat = tm.mat; vec = tm.vec; endvec = tm.endvec;
845  fac = (T)1; mut = true; return *this; }
846 
847  F_TSMatrix<T>& operator *= (const T& f) { fac *= f; return *this; }
848  F_TSMatrix<T>& operator /= (const T& f) { fac /= f; return *this; }
849  F_TSMatrix<T>& operator * (const T& f) { fac *= f; return *this; }
850  F_TSMatrix<T>& operator / (const T& f) { fac /= f; return *this; }
851 
852  F_TSMatrix<T>& operator - () { fac *= (T)(-1); return *this; }
853 
854  // Now for the more complicated
855  F_TMatrix<T> operator + (const F_Matrix<T>&);
856  F_TMatrix<T> operator + (const F_TMatrix<T>&);
857  F_TMatrix<T> operator + (F_TSMatrix<T>);
858  F_TMatrix<T> operator + (const T&);
859  F_TMatrix<T> operator - (const F_Matrix<T>&);
860  F_TMatrix<T> operator - (const F_TMatrix<T>&);
861  F_TMatrix<T> operator - (F_TSMatrix<T>);
862  F_TMatrix<T> operator - (const T&);
863 
864 #ifndef HAVE_GCC295_FRIEND_BUG
865  friend NOINST F_TSMatrix<T> FRIEND_TBCI__ operator * FGD (const T&, F_TSMatrix<T>);
866  friend NOINST F_TMatrix<T> FRIEND_TBCI__ operator + FGD (const T&, F_TSMatrix<T>);
867  friend NOINST F_TMatrix<T> FRIEND_TBCI__ operator - FGD (const T&, F_TSMatrix<T>);
868 #endif
869 
870  F_TMatrix<T> operator * (const F_Matrix<T>& m);
871  F_TMatrix<T> operator * (const CSCMatrix<T>& m);
872  TVector<T> operator * (const Vector<T>& v);
873 
874  bool operator == (const F_Matrix<T>&);
875  bool operator != (const F_Matrix<T>& m) { return !(*this == m); }
876 
877  bool operator == (F_TSMatrix<T>);
878  bool operator != (F_TSMatrix<T> ts) { return !(*this == ts); }
879 
880  bool operator == (F_TMatrix<T> tm) { F_Matrix<T> m(tm); return (*this == m); }
881  bool operator != (F_TMatrix<T> tm) { return !(*this == tm); }
882 
883 
884  //Friends
885  //friend double fabs FGD (F_TSMatrix<T>&);
886  double fabs ();
887 };
888 
889 
890 template <typename T>
891 F_TMatrix<T> F_TSMatrix<T>::operator * (const CSCMatrix<T>& m2)
892 {
893  F_TMatrix<T> m1(*this);
894  F_TMatrix<T> erg(0, m1.rows(), m2.columns());
895  for (unsigned int i=0; i<m1.rows(); i++)
896  for (unsigned int j=0; j<m2.columns(); j++)
897  {
898  T val = (T)0;
899  // FIXME: Do we multiply by fac twice here?
900  for (unsigned int x=m2.pcol[j]; x<m2.pcol[j+1]; x++)
901  val += m1(i, m2.irow[x] )*m2.comp[x];
902  erg.setval(fac*val, i, j);
903  }
904 
905  // Really destroy m1!
906  mut = true;
907  destroy();
908  return erg;
909 }
910 
911 
912 
913 template <typename T>
914 inline void F_TSMatrix<T>::destroy ()
915 {
916  if (!mut || !dim) return;
917  if (vec) TBCIDELETE(T, vec,dim); vec = (T*)0; endvec = (T*)0;
918  if (mat) TBCIDELETE(T*,mat,col); dim = 0;
919 }
920 
921 template <typename T>
922 inline void F_TSMatrix<T>::detach (F_TMatrix<T>* tm)
923 {
924  if ((!tm && mut) || !dim) return;
925  if (!tm) //mut = false
926  {
927  vec = NEW(T ,dim); if (!vec) { dim = 0; mat = (T**)0; endvec = (T*)0; return; }
928  mat = NEW(T*,col); if (!mat) { TBCIDELETE(T,vec,dim); dim = 0; vec = (T*)0; endvec = (T*)0; return; }
929  for (register unsigned c = 0; c < row; c++) mat[c] = vec + c*row;
930  endvec = vec + dim;
931  }
932  else //tm = true, mut = ?
933  {
934  vec = tm->vec; mat = tm->mat; endvec = tm->endvec;
935  }
936  mut = true;
937 }
938 
939 template <typename T>
940 void F_TSMatrix<T>::clone (bool evl, F_TMatrix<T>* tm)
941 {
942  if ((mut && !tm)|| !dim) return;
943  T* oldv = vec; T** oldm = mat; bool omut = mut;
944  detach (tm);
945  if (evl && fac != (T)1)
946  {
947  for (register T *p1 = vec, *p2 = oldv; p1 < endvec; p1++, p2++)
948  *p1 = *p2 * fac;
949  fac = (T)1;
950  }
951  else TBCICOPY (vec, oldv, T, dim);
952  if (tm && omut)
953  { TBCIDELETE(T,oldv,dim); TBCIDELETE(T*,oldm,col); }
954 }
955 
956 
957 
958 //creates mutable evaluated F_TSMatrix (can be easily converted to F_TMatrix)
959 template <typename T>
960 inline F_TSMatrix<T>& F_TSMatrix<T>::eval (F_TMatrix<T>* tm)
961 {
962  if (mut && !tm)
963  {
964  if (fac != (T)1)
965  {
966  for (register T* ptr = vec; ptr < endvec; ptr++)
967  *ptr *= fac;
968  fac = (T)1;
969  }
970  }
971  else
972  clone (true, tm);
973  return (*this);
974 }
975 
976 
977 // Arithmetics
978 
979 #define F_STDDEF_TSM(op) \
980 template <typename T> \
981 inline F_TMatrix<T> F_TSMatrix<T>::operator op (const F_Matrix<T>& m) \
982 { \
983  BCHK(row != m.row || col != m.col, MatErr, Op op on diff size Mats, m.row, F_TMatrix<T> (*this)); \
984  F_TSMatrix<T> ts (*this); ts.detach (); \
985  for (register T *p1 = ts.vec, *p2 = vec, *p3 = m.vec; p2 < endvec; p1++, p2++, p3++) \
986  *p1 = *p2 * fac op *p3; \
987  ts.fac = (T)1; return F_TMatrix<T> (ts); \
988 }
989 
990 F_STDDEF_TSM(+)
991 F_STDDEF_TSM(-)
992 
993 
994 #define F_STDDEF_TSTM(op) \
995 template <typename T> \
996 inline F_TMatrix<T> F_TSMatrix<T>::operator op (const F_TMatrix<T>& tm) \
997 { \
998  BCHK(row != tm.row || col != tm.col, MatErr, Op op on diff size Mats, tm.row, tm); \
999  for (register T *p1 = tm.vec, *p2 = vec; p2 < endvec; p1++, p2++) \
1000  *p1 = *p2 * fac op *p1; \
1001  destroy (); return tm; \
1002 }
1003 
1004 F_STDDEF_TSTM(+)
1005 F_STDDEF_TSTM(-)
1006 
1007 
1008 #define F_STDDEF_TSTS(op) \
1009 template <typename T> \
1010 inline F_TMatrix<T> F_TSMatrix<T>::operator op (F_TSMatrix<T> ts) \
1011 { \
1012  BCHK(row != ts.row || col != ts.col, MatErr, Op op on diff size Mats, ts.row, F_TMatrix<T> (*this)); \
1013  F_TSMatrix<T> tm; \
1014  if (!mut && ts.mut) tm = ts; \
1015  else { tm = *this; tm.detach (); } \
1016  for (register T *p1 = tm.vec, *p2 = vec, *p3 = ts.vec; p2 < endvec; p1++, p2++, p3++) \
1017  *p1 = *p2 * fac op *p3 * ts.fac; \
1018  if (!mut && ts.mut) destroy (); else ts.destroy (); \
1019  tm.fac = (T)1; return F_TMatrix<T> (tm); \
1020 }
1021 
1022 F_STDDEF_TSTS(+)
1023 F_STDDEF_TSTS(-)
1024 
1025 
1026 #define F_STDDEF_TST(op) \
1027 template <typename T> \
1028 inline F_TMatrix<T> F_TSMatrix<T>::operator op (const T& a) \
1029 { \
1030  F_TSMatrix<T> ts (*this); ts.detach (); \
1031  for (register T *p1 = ts.vec, *p2 = vec; p2 < endvec; p1++, p2++) \
1032  *p1 = *p2 * fac op a; \
1033  ts.fac = (T)1; return F_TMatrix<T> (ts); \
1034 }
1035 
1036 F_STDDEF_TST(+)
1037 F_STDDEF_TST(-)
1038 
1039 // operate on diag only?
1040 #define F_STDDEF_TTS(op) \
1041 INST(template <typename T> class F_TSMatrix friend F_TMatrix<T> operator op (const T&, F_TSMatrix<T>);) \
1042 template <typename T> \
1043 inline F_TMatrix<T> operator op (const T& a, F_TSMatrix<T> ts) \
1044 { \
1045  F_TSMatrix<T> tm (ts); tm.detach (); \
1046  for (register T *p1 = tm._VEC, *p2 = ts._VEC; p1 < tm._ENDVEC; p1++, p2++) \
1047  *p1 = a op *p2 * ts._FAC; \
1048  tm._FAC = (T)1; return F_TMatrix<T> (tm); \
1049 }
1050 
1051 F_STDDEF_TTS(+)
1052 F_STDDEF_TTS(-)
1053 
1054 
1055 template <typename T>
1056 F_TMatrix<T> F_TSMatrix<T>::operator * (const F_Matrix<T>& m)
1057 {
1058  F_TMatrix<T> ftm (row, m.col);
1059  BCHK(col != m.row, MatErr, Illegal sizes in F_Matrix multiplication, 0, (destroy(), ftm));
1060  register T tmp ALIGN(MIN_ALIGN);
1061  for (unsigned int r=0; r<row; r++) {
1062  for (unsigned int c=0; c<m.col; c++) {
1063  tmp = mat[0][r] * m(0, c);
1064  for (register unsigned int l=1; l<col; l++)
1065  tmp += mat[l][r] * m(l, c);
1066  ftm.setval(tmp*fac, r, c);
1067  }
1068  }
1069  destroy (); return ftm;
1070 }
1071 
1072 
1073 template <typename T>
1074 TVector<T> F_TSMatrix<T>::operator * (const Vector<T>& v)
1075 {
1076  TVector<T> tv (row);
1077  BCHK (v.size() != col, MatErr, multiply F_Matrix w/ Vector: wrong config, v.size(), (destroy(), tv));
1078  /* FIXME: We could call do_mat_vec_transmult<>() here */
1079  for (unsigned int r=0; r<row; r++)
1080  {
1081  ALIGN3(register T el, mat[0][r] * v(0), MIN_ALIGN2);
1082  for (register unsigned int c=1; c<col; c++)
1083  el += mat[c][r] * v(c);
1084  tv.set (el * fac, r);
1085  }
1086  destroy (); return tv;
1087 }
1088 
1089 
1090 template <typename T>
1091 bool F_TSMatrix<T>::operator == (const F_Matrix<T>& m)
1092 {
1093  if (row != m.row || col != m.col) { destroy (); return false; }
1094  if (dim == 0) { /*destroy();*/ return true; }
1095  if (fac == (T)1)
1096  {
1097  if (vec == m.vec) return true;
1098  if (TBCICOMP (vec, m.vec, T, dim)) { destroy (); return false; }
1099  }
1100  else
1101  {
1102  if (vec == m.vec) return false;
1103  for (register T *p1 = vec, *p2 = m.vec; p1 < endvec; p1++, p2++)
1104  if (*p1 * fac != *p2) { destroy (); return false; }
1105  }
1106  destroy (); return true;
1107 }
1108 
1109 
1110 template <typename T>
1111 bool F_TSMatrix<T>::operator == (F_TSMatrix<T> ts)
1112 {
1113  if (row != ts.row || col != ts.col) { destroy (); ts.destroy (); return false; }
1114  if (dim == 0) {/*destroy(); ts.destroy();*/ return true; }
1115  if (fac == ts.fac)
1116  {
1117  if (vec == ts.vec) { destroy (); return true; }
1118  if (TBCICOMP (vec, ts.vec, T, dim)) { destroy (); ts.destroy (); return false; }
1119  }
1120  else
1121  {
1122  if (vec == ts.vec) { destroy (); return false; }
1123  for (register T *p1 = vec, *p2 = ts.vec; p1 < endvec; p1++, p2++)
1124  if (*p1 * fac != *p2) { destroy (); ts.destroy (); return false; }
1125  }
1126  destroy (); ts.destroy (); return true;
1127 }
1128 
1129 template <typename T>
1130 F_TSMatrix<T>::operator TMatrix<T>()
1131 {
1132  TMatrix<T> tm(this->row, this->col);
1133  for (unsigned int c = 0; c < this->col; ++c)
1134  for (unsigned int r = 0; r < this->row; ++r)
1135  tm.setval(this->get(r, c), r, c);
1136  destroy();
1137  return tm;
1138 }
1139 
1140 
1141 template <typename T>
1142 double F_TSMatrix<T>::fabs ()
1143 {
1144  register T* ptr = vec;
1145  double register rv ALIGN(8) = fabssqr (*ptr++);
1146  for (; ptr < endvec; ptr++)
1147  rv += fabssqr (*ptr);
1148  rv *= fabssqr (fac);
1149  destroy (); return MATH__ sqrt (rv);
1150 }
1151 
1152 
1153 INST(template <typename T> class F_TSMatrix friend F_TSMatrix<T> operator * (const T&, F_TSMatrix<T>);)
1154 template <typename T>
1155 inline F_TSMatrix<T> operator * (const T& f, F_TSMatrix<T> ts)
1156 { ts._FAC = f * ts._FAC; return ts; }
1157 
1158 //---------------------------------------------------------------
1159 // Real, user referable, matrix
1160 //---------------------------------------------------------------
1169 template <typename T>
1170 class F_Matrix : public F_TMatrix<T>
1171 {
1172 
1173  public:
1174 
1175  typedef T value_type;
1176  typedef T element_type;
1177  typedef T aligned_value_type TALIGN(MIN_ALIGN2);
1178 
1179  friend class F_TSMatrix<T>;
1180 
1181  // constructor, destructor
1182  explicit F_Matrix (const unsigned d=0) : F_TMatrix<T> (d) {}
1183  F_Matrix (const unsigned r, const unsigned c) : F_TMatrix<T> (r, c) {}
1184  F_Matrix (const T& v, const unsigned r, const unsigned c) : F_TMatrix<T> (v, r, c) {}
1185  //F_Matrix (const Vector<T>& v, const enum rowcolvec r = colvec)
1186  // : F_TMatrix<T> (v, r) {}
1187 
1188  //copy
1189  F_Matrix (const F_Matrix<T>& m) : F_TMatrix<T> (m) {}
1190  F_Matrix (const CSCMatrix<T>& m);
1191 
1192  //alias
1193  F_Matrix (const F_TMatrix<T>& tm) : F_TMatrix<T> (tm) {}
1194  F_Matrix (F_TSMatrix<T>& ts) : F_TMatrix<T> (ts) {}
1195 
1196  ~F_Matrix () { this->destroy (); }
1197 
1198  // Conversion
1199  F_Matrix (const Matrix<T>& m);
1200  operator TMatrix<T>() const;
1201 
1202  // allow instantiation (Matrix_Sig)
1203  /*virtual*/ static const char* mat_info () { return ("F_Matrix"); }
1204 
1205  // Access
1206  T& operator () (const unsigned int, const unsigned int) const;
1207  //T operator () (const unsigned, const unsigned) const;
1208  T* get_fortran_matrix() const { return this->vec; }
1209  const T& get (const unsigned r, const unsigned c) const
1210  { return this->mat[c][r]; }
1211 
1212  unsigned int columns () const { return this->row; }
1213  unsigned int rows () const { return this->col; }
1214  unsigned long size () const { return this->dim; }
1215  //TVector<T> operator () (const unsigned int) const;
1216  TVector<T> get_row (const unsigned int) const;
1217  TVector<T> get_col (const unsigned int) const;
1218  void set_row (const Vector<T>&, const unsigned int);
1219  void set_col (const Vector<T>&, const unsigned int);
1220  void setval(const T val, const unsigned int r, const unsigned int c)
1221  { this->operator() (r, c) = val; }
1222 
1223  // basics
1224  F_Matrix<T>& operator = (const F_Matrix<T>& m)
1225  { this->F_TMatrix<T>::operator = (m); return *this; }
1226  F_Matrix<T>& operator = (const F_TMatrix<T>& tm)
1227  { this->F_TMatrix<T>::operator = (tm); return *this; }
1228  F_Matrix<T>& operator = (F_TSMatrix<T> ts)
1229  { F_TMatrix<T>(this->F_TMatrix<T>::operator = (ts)); return *this; }
1230  F_Matrix<T>& operator = (const T& v)
1231  { this->F_TMatrix<T>::operator = (v); return *this; }
1232 
1233  F_Matrix<T>& resize (const unsigned r, const unsigned c)
1234  { this->F_TMatrix<T>::resize (r, c); return *this; }
1235  F_Matrix<T>& resize (const unsigned d)
1236  { return this->resize (d, d); }
1237  F_Matrix<T>& resize (const T& v, const unsigned r, const unsigned c)
1238  { this->F_TMatrix<T>::resize (v, r, c); return *this; }
1239  F_Matrix<T>& fill (const T& v = 0)
1240  { this->F_TMatrix<T>::fill (v); return *this; }
1241  F_Matrix<T>& setunit (const T& f = 1)
1242  { this->F_TMatrix<T>::setunit (f); return *this; }
1243  F_Matrix<T>& clear ()
1244  { this->F_TMatrix<T>::clear (); return *this; }
1245 
1246  // assignment ops (do we really need em ?)
1247  // Seems we have to make sure, that a F_Matrix and no F_TMatrix
1248  // is returned; otherwise there's danger to be destroyed
1249  F_Matrix<T>& operator += (const F_Matrix<T>& a)
1250  { this->F_TMatrix<T>::operator += (a); return *this; }
1251  F_Matrix<T>& operator -= (const F_Matrix<T>& a)
1252  { this->F_TMatrix<T>::operator -= (a); return *this; }
1253  F_Matrix<T>& operator += (F_TMatrix<T> a)
1254  { this->F_TMatrix<T>::operator += (a); return *this; }
1255  F_Matrix<T>& operator -= (F_TMatrix<T> a)
1256  { this->F_TMatrix<T>::operator -= (a); return *this; }
1257  F_Matrix<T>& operator += (const T& a)
1258  { this->F_TMatrix<T>::operator += (a); return *this; }
1259  F_Matrix<T>& operator -= (const T& a)
1260  { this->F_TMatrix<T>::operator -= (a); return *this; }
1261  F_Matrix<T>& operator *= (const T& a)
1262  { this->F_TMatrix<T>::operator *= (a); return *this; }
1263  F_Matrix<T>& operator /= (const T& a)
1264  { this->F_TMatrix<T>::operator /= (a); return *this; }
1265 
1266 
1267  // basic operations
1268  F_TMatrix<T> operator - () const { return -(F_TMatrix<T> (*this)); }
1269 
1270  F_TMatrix<T> operator + (const F_Matrix<T>&) const;
1271  F_TMatrix<T> operator - (const F_Matrix<T>&) const;
1272  F_TMatrix<T> operator * (const F_Matrix<T>&) const;
1273 
1274  F_TMatrix<T> operator + (F_TMatrix<T>) const;
1275  F_TMatrix<T> operator - (F_TMatrix<T>) const;
1276  F_TMatrix<T> operator * (F_TMatrix<T>) const;
1277 
1278  F_TMatrix<T> operator + (F_TSMatrix<T>) const;
1279  F_TMatrix<T> operator - (F_TSMatrix<T>) const;
1280  F_TMatrix<T> operator * (F_TSMatrix<T>) const;
1281 
1282  F_TMatrix<T> operator + (const T&) const;
1283  F_TMatrix<T> operator - (const T&) const;
1284  F_TSMatrix<T> operator * (const T&) const;
1285  F_TSMatrix<T> operator / (const T&) const;
1286 
1287  TVector<T> operator * (const Vector<T>&) const;
1288  TVector<T> operator * (TVector<T>&) const;
1289 
1290  bool operator == (const F_Matrix<T>& m) const;
1291  bool operator != (const F_Matrix<T>& m) const
1292  { return !(*this == m); }
1293 
1294  bool operator == (F_TMatrix<T>& tm) const
1295  { F_Matrix<T> m (tm); return (*this == m); }
1296  bool operator != (F_TMatrix<T>& tm) const
1297  { return !(*this == tm); }
1298 
1299  bool operator == (F_TSMatrix<T>&) const;
1300  bool operator != (F_TSMatrix<T>& ts) const
1301  { return !(*this == ts); }
1302 
1303  // io-streams
1304 
1305  friend STD__ ostream& operator << FGD (STD__ ostream&, const F_Matrix<T>&);
1306  friend STD__ istream& operator >> FGD (STD__ istream&, F_Matrix<T>&);
1307 
1308  // misc
1309  // returns the the hermitian adjungated (Transposed and conjugate complex)
1310  F_TMatrix<T> herm () const { return (F_TMatrix<T> (*this)).herm(); }
1311  F_TMatrix<T> conj () const { return (F_TMatrix<T> (*this)).conj(); }
1312  F_TMatrix<T> trans() const { return (F_TMatrix<T> (*this)).trans(); }
1313 
1314  double fabs () const;
1315 };
1316 
1317 
1318 // definitions
1319 
1320 template <typename T>
1321 inline F_Matrix<T>::F_Matrix (const CSCMatrix<T>& m)
1322 : F_TMatrix<T>(m.rows(), m.columns())
1323 {
1324  for (unsigned int r=0; r<m.rows(); r++)
1325  for (unsigned int c=0; c<m.columns(); c++)
1326  operator()(r,c) = m(r,c);
1327 }
1328 
1329 template <typename T>
1330 F_Matrix<T>::F_Matrix(const Matrix<T>& m)
1331 : F_TMatrix<T>(m.rows(), m.columns())
1332 {
1333  for (unsigned int r = 0; r < this->row; ++r)
1334  for (unsigned int c = 0; c < this->col; ++c)
1335  setval(m(r, c), r, c);
1336 }
1337 
1338 template <typename T>
1339 F_Matrix<T>::operator TMatrix<T>() const
1340 {
1341  TMatrix<T> tm(this->row, this->col);
1342  for (unsigned int c = 0; c < this->col; ++c)
1343  for (unsigned int r = 0; r < this->row; ++r)
1344  tm.setval((*this)(r, c), r, c);
1345  return tm;
1346 }
1347 
1348 
1349 template <typename T>
1350 STD__ ostream& operator << (STD__ ostream& os, const F_Matrix<T>& m)
1351 {
1352  for (unsigned int r=0; r<m.row; r++) {
1353  for (unsigned int c=0; c<m.col; c++)
1354  os << m(r,c) << " ";
1355  os << "\n";
1356  }
1357  return os;
1358 }
1359 
1360 INST(template <typename T> class F_TMatrix friend STD__ ostream& operator << (STD__ ostream& os, F_TMatrix<T> tm);)
1361 template <typename T>
1362 STD__ ostream& operator << (STD__ ostream& os, F_TMatrix<T> tm)
1363 {
1364  F_Matrix<T> fm(tm);
1365  return os << fm;
1366 }
1367 
1368 #if 1
1369 INST(template <typename T> class F_TSMatrix friend STD__ ostream& operator << (STD__ ostream& os, F_TSMatrix<T> ts);)
1370 template <typename T>
1371 STD__ ostream& operator << (STD__ ostream& os, F_TSMatrix<T> ts)
1372 {
1373  return os << F_Matrix<T>(ts);
1374 }
1375 #endif
1376 
1377 template <typename T>
1378 STD__ istream& operator >> (STD__ istream& in, F_Matrix<T>& m)
1379 {
1380  Vector<T> buf(m.row);
1381 
1382  for (unsigned int i=0; i<m.col; i++)
1383  {
1384  in >> buf;
1385  m.set_col (buf, i);
1386  }
1387  return in;
1388 }
1389 
1390 
1391 /************** Arithmetics ***************/
1392 
1393 //change these to operate on diag only ?
1394 #define F_MSTDDEF_T(op) \
1395 template <typename T> \
1396 inline F_TMatrix<T> F_Matrix<T>::operator op (const T& a) const \
1397 { \
1398  F_TMatrix<T> t (this->row, this->col); \
1399  for (register T *p1=t.vec, *p2=this->vec; p2<this->endvec; ++p1, ++p2) \
1400  *p1 = *p2 op a; \
1401  return t; \
1402 }
1403 
1404 F_MSTDDEF_T(+)
1405 F_MSTDDEF_T(-)
1406 
1407 template <typename T>
1408 inline F_TSMatrix<T> F_Matrix<T>::operator * (const T& a) const
1409 { return F_TSMatrix<T> (*this, a); }
1410 
1411 template <typename T>
1412 inline F_TSMatrix<T> F_Matrix<T>::operator / (const T& a) const
1413 {
1414  BCHK (a==(T)0, MatErr, Divide Mat by 0, 0, *this);
1415  return F_TSMatrix<T> (*this, (T)1/a);
1416 }
1417 
1418 #define F_MSTDDEF_M(op) \
1419 template <typename T> \
1420 inline F_TMatrix<T> F_Matrix<T>::operator op (const F_Matrix<T>& a) const \
1421 { \
1422  F_TMatrix<T> t (this->row, this->col); \
1423  BCHK(this->dim != a.dim, MatErr, Operator op on diff size matrices, a.dim, t); \
1424  for (register T *p1=t.vec, *p2=this->vec, *p3=a.vec; p2<this->endvec; ++p1, ++p2, ++p3) \
1425  *p1 = *p2 op *p3; \
1426  return t; \
1427 }
1428 
1429 F_MSTDDEF_M(+)
1430 F_MSTDDEF_M(-)
1431 
1432 
1433 #define F_MSTDDEF_TM(op) \
1434 template <typename T> \
1435 inline F_TMatrix<T> F_Matrix<T>::operator op (F_TMatrix<T> a) const \
1436 { \
1437  BCHK(this->dim != a.dim, MatErr, Applying op on diff. size matrices, a.dim, a); \
1438  for (register T *p1=a.vec, *p2=this->vec; p2<this->endvec; ++p1, ++p2) \
1439  *p1 = *p2 op *p1; \
1440  return a; \
1441 }
1442 
1443 F_MSTDDEF_TM(+)
1444 F_MSTDDEF_TM(-)
1445 
1446 #define F_MSTDDEF_TS(op) \
1447 template <typename T> \
1448 inline F_TMatrix<T> F_Matrix<T>::operator op (F_TSMatrix<T> ts) const \
1449 { \
1450  BCHK(this->dim != ts.dim, MatErr, Applying op on diff. size matrices, ts.dim, (ts.destroy(),F_TMatrix<T> (*this)) ); \
1451  F_TSMatrix<T> tm (ts); tm.detach (); \
1452  for (register T *p1=tm.vec, *p2=this->vec, *p3=ts.vec; p2<this->endvec; ++p1, ++p2, ++p3) \
1453  *p1 = *p2 op *p3 * ts.fac; \
1454  tm.fac = (T)1; return F_TMatrix<T> (tm); \
1455 }
1456 
1457 F_MSTDDEF_TS(+)
1458 F_MSTDDEF_TS(-)
1459 
1460 
1461 // F_Matrix multiplication
1462 template <typename T>
1463 F_TMatrix<T> F_Matrix<T>::operator * (const F_Matrix<T>& a) const
1464 {
1465  F_TMatrix<T> ftm(this->row, a.col);
1466  BCHK(this->col != a.row, MatErr, Illegal sizes in F_Matrix multiplication, 0, ftm);
1467  for (unsigned int r=0; r<this->row; ++r) {
1468  for (register unsigned int c=0; c<a.col; ++c) {
1469  ALIGN3(register T tmp, (*this)(r, 0) * a(0, c), MIN_ALIGN2);
1470  for (register unsigned int l=1; l<this->col; ++l)
1471  tmp += (*this)(r, l) * a(l, c);
1472  ftm.setval(tmp, r, c);
1473  }
1474  }
1475  return ftm;
1476 }
1477 
1478 template <typename T>
1479 F_TMatrix<T> F_Matrix<T>::operator * (F_TSMatrix<T> a) const
1480 {
1481  F_TMatrix<T> ftm(this->row, a.col);
1482  BCHK(this->col != a.row, MatErr, Illegal sizes in F_Matrix multiplication, 0, ftm);
1483  register T tmp ALIGN(MIN_ALIGN2);
1484  for (unsigned int r=0; r<this->row; ++r) {
1485  for (register unsigned int c=0; c<a.col; ++c) {
1486  tmp = (*this)(c, 0) * a.mat[c][0];
1487  for (register unsigned int l=1; l<this->col; ++l)
1488  tmp += (*this)(r, l) * a.mat[c][l];
1489  ftm.setval(tmp*a.fac, r, c);
1490  }
1491  }
1492  a.destroy (); return ftm;
1493 }
1494 
1495 
1496 template <typename T>
1497 inline F_TMatrix<T> F_Matrix<T>::operator * (F_TMatrix<T> a) const
1498 {
1499  F_TMatrix<T> ftm (this->operator * (F_Matrix<T> (a)));
1500  return ftm;
1501 }
1502 
1503 template <typename T>
1504 inline F_TMatrix<T> F_TMatrix<T>::operator * (const F_Matrix<T>& a)
1505 {
1506  F_Matrix<T> m (*this);
1507  return (m * a);
1508 }
1509 
1510 
1511 template <typename T>
1512 inline F_TMatrix<T> F_TMatrix<T>::operator * (F_TMatrix<T> a)
1513 {
1514  F_Matrix<T> m (*this);
1515  return (m * (F_Matrix<T> (a)));
1516 }
1517 
1518 template <typename T>
1519 inline F_TMatrix<T> F_TMatrix<T>::operator * (F_TSMatrix<T> a)
1520 {
1521  F_Matrix<T> m (*this);
1522  return (m * a);
1523 }
1524 
1525 
1526 template <typename T>
1527 bool F_Matrix<T>::operator == (const F_Matrix<T>& m) const
1528 {
1529  if (this->row != m.row || this->col != m.col)
1530  return false;
1531  if (this->vec == m.vec || this->dim == 0)
1532  return true;
1533  if (TBCICOMP (this->vec, m.vec, T, this->dim))
1534  return false;
1535  else
1536  return true;
1537 }
1538 
1539 template <typename T>
1540 bool F_Matrix<T>::operator == (F_TSMatrix<T>& ts) const
1541 {
1542  if (this->row != ts.row || this->col != ts.col)
1543  return (ts.destroy(), false);
1544  if (this->vec == 0) {
1545  /*ts.destroy();*/ return true;
1546  }
1547  if (ts.fac == (T)1) {
1548  if (this->vec == ts.vec)
1549  return true;
1550  if (TBCICOMP(this->vec, ts.vec, T, this->dim)) {
1551  ts.destroy (); return false;
1552  }
1553  } else {
1554  if (this->vec == ts.vec)
1555  return false;
1556  for (register T *p1 = this->vec, *p2 = ts.vec; p1 < this->endvec; ++p1, ++p2)
1557  if (*p1 != *p2 * ts.fac) {
1558  ts.destroy (); return false;
1559  }
1560  }
1561  ts.destroy ();
1562  return true;
1563 }
1564 
1565 
1566 template<typename T>
1567 TVector<T> F_Matrix<T>::operator * (const Vector<T>& v) const
1568 {
1569  TVector<T> tv (this->row);
1570  for (unsigned int r=0; r<this->row; ++r) {
1571  register T el = 0;
1572  for (register unsigned int c=0; c<this->col; ++c)
1573  el += this->operator() (r, c) * v(c);
1574  tv.set (el, r);
1575  }
1576  return tv;
1577 }
1578 
1579 template<typename T>
1580 inline TVector<T> F_Matrix<T>::operator * (TVector<T>& c) const
1581 {
1582  Vector<T> cn (c);
1583  return (*this * cn);
1584 }
1585 
1586 #if 0
1587 template <typename T>
1588 inline TVector<T> F_Matrix<T>::operator () (const unsigned int i) const
1589 {
1590  TVector<T> v(this->row);
1591  BCHK(i < 0 || i >= this->col, MatErr, Illegal row index in get_row, i, v);
1592  for (register unsigned int j=0; j<this->row; ++j)
1593  v.vec[j] = this->mat[j][i];
1594  return v;
1595 }
1596 #endif
1597 
1598 template <typename T>
1599 inline void F_Matrix<T>::set_col (const Vector<T>& v, const unsigned int c)
1600 {
1601  F_TMatrix<T>::set_col(v, c);
1602 }
1603 
1604 template <typename T>
1605 inline TVector<T> F_Matrix<T>::get_col (const unsigned int c) const
1606 {
1607  return F_TMatrix<T>::get_col(c);
1608 }
1609 
1610 template <typename T>
1611 inline TVector<T> F_Matrix<T>::get_row (const unsigned int r) const
1612 {
1613  return F_TMatrix<T>::get_row(r);
1614 }
1615 
1616 template <typename T>
1617 inline void F_Matrix<T>::set_row (const Vector<T>& v, const unsigned int r)
1618 {
1619  F_TMatrix<T>::set_row(v, r);
1620 }
1621 
1622 template <typename T>
1623 inline T& F_Matrix<T>::operator () (const unsigned int r, const unsigned int c) const
1624 {
1625  BCHK(r>=this->row, MatErr, illegal row index, r, this->mat[0][0]);
1626  BCHK(c>=this->col, MatErr, Illegal col index, c, this->mat[0][0]);
1627  return this->mat[c][r];
1628 }
1629 
1630 /*
1631 template <typename T>
1632 inline T F_Matrix<T>::operator () (const unsigned r, const unsigned c) const
1633 {
1634  BCHK(r<0 || r>=row, MatErr, Illegal row index, r, mat[0][0]);
1635  BCHK(c<0 || c>=col, MatErr, illegal col index, c, mat[0][0]);
1636  return mat[c][r];
1637 }
1638  */
1639 
1640 template <typename T>
1641 double F_Matrix<T>::fabs () const
1642 {
1643  register double rv = 0.0;
1644  for (register T *ptr = this->vec; ptr < this->endvec; ++ptr)
1645  rv += fabssqr (*ptr);
1646  return MATH__ sqrt (rv);
1647 }
1648 
1649 
1650 template <typename T>
1651 F_TMatrix<T> F_TMatrix<T>::transposed_copy() const
1652 {
1653  F_TMatrix<T> trm(this->col, this->row);
1654  for (unsigned r = 0; r < this->row; ++r)
1655  for (unsigned c = 0; c < this->col; ++c)
1656  trm.mat[c][r] = this->mat[r][c];
1657  // memleak if it's really a F_TMatrix :-(
1658  //this->destroy();
1659  //F_Matrix<T> fm(*this);
1660  return trm;
1661 }
1662 
1663 
1664 INST(template <typename T> class F_TMatrix friend F_TMatrix<T> transpose(const F_TMatrix<T>& tm);)
1665 template <typename T>
1666 inline F_TMatrix<T> transpose(const F_TMatrix<T>& ftm)
1667 {
1668  return ftm.transposed_copy();
1669 }
1670 
1671 template <typename T>
1672 F_TMatrix<T>& F_TMatrix<T>::transpose()
1673 {
1674  if (this->col != this->row) {
1675  // for non quadratic shape, doing inplace swap ops is too hard
1676  F_Matrix<T> transmat(this->transposed_copy());
1677  return this->swap(transmat);
1678  } else {
1679  // do in place swap ops
1680  for (unsigned r = 1; r < this->row; ++r)
1681  for (unsigned c = 0; c < r; ++c)
1682  SWAP(this->mat[r][c], this->mat[c][r]);
1683  return *this;
1684  }
1685 }
1686 
1687 template <typename T>
1688 inline F_TMatrix<T>& F_TMatrix<T>::swap (F_TMatrix<T>& m)
1689 {
1690  //BCHK(this->row != m.row || this->col != m.col, MatErr, Swapping Matrices w/ different layout, m.row, *this);
1691  SWAP(this->dim, m.dim);
1692  SWAP(this->row, m.row);
1693  //SWAP(this->col, m.col);
1694  unsigned int tcol = this->col; this->col = m.col; m.col = tcol;
1695  SWAP(this->mat, m.mat);
1696  SWAP(this->vec, m.vec); SWAP(this->endvec, m.endvec);
1697  return *this;
1698 }
1699 
1700 #undef _DIM
1701 #undef _VEC
1702 #undef _ENDVEC
1703 #undef _ROW
1704 #undef _COL
1705 #undef _FAC
1706 
1707 NAMESPACE_END
1708 
1709 NAMESPACE_CSTD
1710 
1711 INST(template <typename T> class TBCI__ F_TMatrix friend double fabs (const TBCI__ F_TMatrix<T>&);)
1712 template <typename T>
1713 inline double fabs (const TBCI__ F_TMatrix<T>& ftm)
1714 { return ftm.fabs (); }
1715 
1716 INST(template <typename T> class TBCI__ F_Matrix friend double fabs (const TBCI__ F_Matrix<T>&);)
1717 template <typename T>
1718 inline double fabs (const TBCI__ F_Matrix<T>& fm)
1719 { return fm.fabs (); }
1720 
1721 INST(template <typename T> class TBCI__ F_TSMatrix friend double fabs (TBCI__ F_TSMatrix<T>&);)
1722 template <typename T>
1723 inline double fabs (/*typename*/ TBCI__ F_TSMatrix<T>& ftsm)
1724 { return ftsm.fabs (); }
1725 
1726 NAMESPACE_CSTD_END
1727 
1728 #endif /* TBCI_F_MATRIX_H */
TBCICOPY
#define TBCICOPY(n, o, t, s)
Definition: basics.h:875
F_TSMatrix::~F_TSMatrix
~F_TSMatrix()
Definition: f_matrix.h:815
F_TMatrix::transposed_copy
F_TMatrix< T > transposed_copy() const
Inefficient! Use transMult if possible.
Definition: f_matrix.h:1651
BCHKNR
#define BCHKNR(cond, exc, txt, ind)
Definition: basics.h:579
a
Matrix< T > a(10, 10)
false
#define false
Definition: bool.h:23
F_TSMatrix::operator*=
F_TSMatrix< T > & operator*=(const T &f)
Definition: f_matrix.h:847
F_TSMatrix::row
unsigned int row
Definition: f_matrix.h:787
MIN_ALIGN
#define MIN_ALIGN
Definition: basics.h:414
ALIGN
#define ALIGN(x)
Definition: basics.h:437
fac
tm fac
Definition: f_matrix.h:1051
F_Matrix::setval
void setval(const T val, const unsigned int r, const unsigned int c)
Definition: f_matrix.h:1220
F_Matrix::fabs
double fabs() const
Definition: f_matrix.h:1641
F_Matrix::operator+
F_TMatrix< T > operator+(const F_Matrix< T > &) const
Definition: f_matrix.h:1429
F_TSMatrix::F_TSMatrix
F_TSMatrix(const F_TMatrix< T > &tm, const T &f=(T) 1)
Definition: f_matrix.h:817
F_TSMatrix::operator-
F_TSMatrix< T > & operator-()
Definition: f_matrix.h:852
F_Matrix::get_col
TVector< T > get_col(const unsigned int) const
Definition: f_matrix.h:1605
F_Matrix::columns
unsigned int columns() const
Definition: f_matrix.h:1212
F_TMatrix::operator!=
bool operator!=(const F_Matrix< T > &m)
Definition: f_matrix.h:185
BdMatrix
The class BdMatrix is an implementation to store and do operations on sparse Matrices with a band str...
Definition: band_matrix.h:103
for
for(register T *p1=c.vec,*p2=b.vec;p1< c.endvec;p1++, p2++)*p1
CSCMatrix::irow
unsigned int * irow
Definition: cscmatrix.h:217
F_Matrix::operator+=
F_Matrix< T > & operator+=(const F_Matrix< T > &a)
Definition: f_matrix.h:1249
c
return c
Definition: f_matrix.h:760
F_TMatrix::operator()
T & operator()(const unsigned int, const unsigned int) const
Definition: f_matrix.h:338
MIN
#define MIN(a, b)
Definition: basics.h:648
FRIEND_TBCI__
#define FRIEND_TBCI__
Definition: basics.h:327
F_Matrix::F_Matrix
F_Matrix(const unsigned d=0)
Definition: f_matrix.h:1182
F_Matrix::get_row
TVector< T > get_row(const unsigned int) const
Definition: f_matrix.h:1611
NAMESPACE_TBCI
#define NAMESPACE_TBCI
Definition: basics.h:310
F_TMatrix::get_row
TVector< T > get_row(const unsigned int) const
Definition: f_matrix.h:449
F_Matrix::operator/
F_TSMatrix< T > operator/(const T &) const
Definition: f_matrix.h:1412
F_TMatrix::trace
T trace() const
Definition: f_matrix.h:517
F_MSTDDEF_T
#define F_MSTDDEF_T(op)
Definition: f_matrix.h:1394
F_Matrix::operator!=
bool operator!=(const F_Matrix< T > &m) const
Definition: f_matrix.h:1291
F_TSMatrix::detach
void detach(F_TMatrix< T > *=0)
Definition: f_matrix.h:922
CSCMatrix::comp
T * comp
Definition: cscmatrix.h:218
Matrix_Sig::rows
unsigned int rows() const
Matrix_Sig
Common interface definition (signature) for all Matrices.
Definition: matrix_sig.h:45
MIN_ALIGN2
#define MIN_ALIGN2
Definition: basics.h:417
F_TMatrix::clear
F_TMatrix< T > & clear()
Definition: f_matrix.h:406
F_Matrix::resize
F_Matrix< T > & resize(const unsigned r, const unsigned c)
Definition: f_matrix.h:1233
BVector::dim
unsigned long dim
Definition: bvector.h:74
F_TSMatrix::dim
unsigned long dim
Definition: f_matrix.h:786
F_Matrix::set_row
void set_row(const Vector< T > &, const unsigned int)
Definition: f_matrix.h:1617
F_TMatrix::rows
unsigned int rows() const
Definition: f_matrix.h:199
BCHK
#define BCHK(cond, exc, txt, ind, rtval)
Definition: basics.h:568
F_Matrix::trans
F_TMatrix< T > trans() const
Definition: f_matrix.h:1312
F_Matrix::get_fortran_matrix
T * get_fortran_matrix() const
Definition: f_matrix.h:1208
F_STDDEF_TTM
#define F_STDDEF_TTM(op)
Definition: f_matrix.h:726
F_Matrix::operator=
F_Matrix< T > & operator=(const F_Matrix< T > &m)
Definition: f_matrix.h:1224
F_TMatrix::element_type
T element_type
Definition: f_matrix.h:84
NAMESPACE_CSTD_END
#define NAMESPACE_CSTD_END
Definition: basics.h:318
F_Matrix::setunit
F_Matrix< T > & setunit(const T &f=1)
Definition: f_matrix.h:1241
CSCMatrix::rows
unsigned int rows() const
query matrix dimensions
Definition: cscmatrix.h:101
F_Matrix::mat_info
static const char * mat_info()
Definition: f_matrix.h:1203
TBCIFILL
#define TBCIFILL(n, v, t, s)
Definition: basics.h:891
F_TSMatrix::clone
void clone(bool=false, F_TMatrix< T > *=0)
Definition: f_matrix.h:940
F_Matrix::TALIGN
T aligned_value_type TALIGN(MIN_ALIGN2)
Definition: f_matrix.h:1177
F_TSMatrix::operator/
F_TSMatrix< T > & operator/(const T &f)
Definition: f_matrix.h:850
F_Matrix::operator*
F_TMatrix< T > operator*(const F_Matrix< T > &) const
Definition: f_matrix.h:1463
F_TMatrix::get_fortran_matrix
T * get_fortran_matrix() const
Definition: f_matrix.h:181
F_TMFORALL_T
#define F_TMFORALL_T(op)
Definition: f_matrix.h:658
x
const Vector< T > const Vector< T > & x
Definition: LM_fit.h:97
F_TSMatrix::TALIGN
T aligned_value_type TALIGN(MIN_ALIGN2)
Definition: f_matrix.h:808
F_Matrix::resize
F_Matrix< T > & resize(const T &v, const unsigned r, const unsigned c)
Definition: f_matrix.h:1237
F_Matrix::operator/=
F_Matrix< T > & operator/=(const T &a)
Definition: f_matrix.h:1263
F_TSMatrix::value_type
T value_type
Definition: f_matrix.h:806
F_Matrix::clear
F_Matrix< T > & clear()
Definition: f_matrix.h:1243
F_TMatrix::resize
F_TMatrix< T > & resize(const unsigned int, const unsigned int)
Definition: f_matrix.h:463
F_STDDEF_TMTM
#define F_STDDEF_TMTM(op)
Definition: f_matrix.h:694
F_Matrix::F_Matrix
F_Matrix(F_TSMatrix< T > &ts)
Definition: f_matrix.h:1194
INST
INST(template< typename T > class F_TMatrix friend F_TMatrix< T > operator+(const T &, F_TMatrix< T >);) template< typename T > inline F_TMatrix< T > operator+(const T &a
F_TSMatrix::F_TSMatrix
F_TSMatrix()
Definition: f_matrix.h:813
F_TMFORALL_M
#define F_TMFORALL_M(op)
Definition: f_matrix.h:616
F_TSMatrix::operator+
F_TMatrix< T > operator+(const F_Matrix< T > &)
Definition: f_matrix.h:990
TVector::set
T & set(const T &val, const unsigned long i) const
Definition: vector.h:194
F_TSMatrix
Temporary object for scaled matrices.
Definition: cscmatrix.h:50
F_TMatrix::destroy
void destroy()
Definition: f_matrix.h:233
TMatrix::setval
T & setval(const T &val, const unsigned int r, const unsigned int c)
Definition: matrix.h:281
F_TMatrix::dim
unsigned long dim
Definition: f_matrix.h:74
F_TMatrix::conj
F_TMatrix< T > & conj()
Definition: f_matrix.h:570
F_TSMatrix::operator/=
F_TSMatrix< T > & operator/=(const T &f)
Definition: f_matrix.h:848
F_Matrix::~F_Matrix
~F_Matrix()
Definition: f_matrix.h:1196
fabs
NAMESPACE_END NAMESPACE_CSTD double fabs(const TBCI__ F_TMatrix< T > &ftm)
Definition: f_matrix.h:1713
F_TMatrix::setval
void setval(const T val, const unsigned int r, const unsigned int c)
Definition: f_matrix.h:206
F_TMatrix::operator-=
F_TMatrix< T > & operator-=(F_TMatrix< T >)
Definition: f_matrix.h:639
F_Matrix::F_Matrix
F_Matrix(const T &v, const unsigned r, const unsigned c)
Definition: f_matrix.h:1184
F_TMatrix::operator/=
F_TSMatrix< T > operator/=(const T &)
Definition: f_matrix.h:677
F_TMatrix::set_ptrs
int set_ptrs()
Definition: f_matrix.h:254
NEW
#define NEW(t, s)
Definition: malloc_cache.h:633
std::sqrt
double sqrt(const int a)
Definition: basics.h:1197
CTensor::fill
CTensor< T > & fill(const T &value)
Definition: tensor.h:140
F_Matrix::fill
F_Matrix< T > & fill(const T &v=0)
Definition: f_matrix.h:1239
F_TMatrix::setunit
F_TMatrix< T > & setunit(const T &=1)
Definition: f_matrix.h:396
F_STDDEF_TSTS
#define F_STDDEF_TSTS(op)
Definition: f_matrix.h:1008
F_TMatrix::operator/
F_TSMatrix< T > operator/(const T &)
Definition: f_matrix.h:717
F_Matrix::resize
F_Matrix< T > & resize(const unsigned d)
Definition: f_matrix.h:1235
F_TMatrix::operator*
F_TSMatrix< T > operator*(const T &)
Definition: f_matrix.h:713
F_TMatrix::mat_t
T mat_t
Definition: f_matrix.h:73
F_TSMatrix::element_type
T element_type
Definition: f_matrix.h:807
F_Matrix::conj
F_TMatrix< T > conj() const
Definition: f_matrix.h:1311
Tensor
Tensor class including arithmetics.
Definition: f_matrix.h:53
ts
F_TSMatrix< T > ts
Definition: f_matrix.h:1051
TBCI__
#define TBCI__
Definition: basics.h:325
F_TSMatrix::operator*
F_TSMatrix< T > & operator*(const T &f)
Definition: f_matrix.h:849
F_MSTDDEF_TS
#define F_MSTDDEF_TS(op)
Definition: f_matrix.h:1446
index
long int Vector< T > & index
Definition: LM_fit.h:69
F_TMatrix::row
unsigned int row
Definition: f_matrix.h:75
F_Matrix::operator*=
F_Matrix< T > & operator*=(const T &a)
Definition: f_matrix.h:1261
F_TMatrix::operator+
F_TMatrix< T > & operator+(F_TMatrix< T >)
Definition: f_matrix.h:699
b
F_TMatrix< T > b
Definition: f_matrix.h:736
F_TSMatrix::mat
T ** mat
Definition: f_matrix.h:788
F_TSMatrix::F_TSMatrix
F_TSMatrix(const F_Matrix< T > &m, const T &f=(T) 1)
Definition: f_matrix.h:820
F_Matrix
Definition: bvector.h:51
F_TMatrix::mat_info
static const char * mat_info()
Definition: f_matrix.h:119
F_STDDEF_TMT
#define F_STDDEF_TMT(op)
Definition: f_matrix.h:703
F_TMatrix::F_TMatrix
F_TMatrix(const unsigned=0)
Definition: f_matrix.h:275
Matrix_Sig::columns
unsigned int columns() const
F_TMatrix::get_col
TVector< T > get_col(const unsigned int) const
Definition: f_matrix.h:432
MatErr
exception class
Definition: matrix.h:53
F_Matrix::rows
unsigned int rows() const
Definition: f_matrix.h:1213
F_TMatrix::value_type
T value_type
Definition: f_matrix.h:83
NAMESPACE_CSTD
#define NAMESPACE_CSTD
Definition: basics.h:312
F_TMatrix::size
unsigned long size() const
Definition: f_matrix.h:200
F_TMatrix::operator==
bool operator==(const F_Matrix< T > &m)
Definition: f_matrix.h:586
F_TSMatrix::mat_info
static const char * mat_info()
Definition: f_matrix.h:835
TBCI::SWAP
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:794
F_TSMatrix::operator!=
bool operator!=(const F_Matrix< T > &m)
Definition: f_matrix.h:875
F_Matrix::operator-=
F_Matrix< T > & operator-=(const F_Matrix< T > &a)
Definition: f_matrix.h:1251
F_TMatrix::alias
F_TMatrix< T > & alias(const F_Matrix< T > &m)
Definition: f_matrix.h:127
F_TSMatrix::ALIGN
T fac ALIGN(MIN_ALIGN)
TVector::size
unsigned long size() const
Definition: vector.h:104
F_Matrix::F_Matrix
F_Matrix(const unsigned r, const unsigned c)
Definition: f_matrix.h:1183
FGD
#define FGD
Definition: basics.h:137
F_TMatrix::columns
unsigned int columns() const
Definition: f_matrix.h:198
F_TMatrix::col
unsigned int col
Definition: f_matrix.h:75
F_TSMatrix::vec
T * vec
Definition: f_matrix.h:789
Matrix
Definition: bvector.h:49
detach
tm detach()
F_TMatrix::transpose
F_TMatrix< T > & transpose()
Definition: f_matrix.h:1672
F_Matrix::operator-
F_TMatrix< T > operator-() const
Definition: f_matrix.h:1268
CSCMatrix::pcol
unsigned int * pcol
Definition: cscmatrix.h:216
F_MSTDDEF_TM
#define F_MSTDDEF_TM(op)
Definition: f_matrix.h:1433
F_TMatrix::trans
F_TMatrix< T > & trans()
Definition: f_matrix.h:530
F_TMatrix
Temporary Base Class (non referable!) (acc.
Definition: bvector.h:50
i
int i
Definition: LM_fit.h:71
F_STDDEF_TM
#define F_STDDEF_TM(op)
Definition: f_matrix.h:749
F_TMatrix::fill
F_TMatrix< T > & fill(const T &=0)
Definition: f_matrix.h:510
F_MSTDDEF_M
#define F_MSTDDEF_M(op)
Definition: f_matrix.h:1418
LAPACK_INLINE
#define LAPACK_INLINE
Definition: f_matrix.h:460
CSCMatrix
exception class: Use MatErr from matrix.h
Definition: cscmatrix.h:49
F_TSMatrix::fabs
double fabs()
Definition: f_matrix.h:1142
F_Matrix::size
unsigned long size() const
Definition: f_matrix.h:1214
STD__
#define STD__
Definition: basics.h:331
F_TSMatrix::destroy
void destroy()
Definition: f_matrix.h:914
F_TSMatrix::operator=
F_TSMatrix< T > & operator=(const F_TSMatrix< T > &ts)
Definition: f_matrix.h:838
TBCIDELETE
#define TBCIDELETE(t, v, sz)
Definition: malloc_cache.h:634
F_TMatrix::set_col
void set_col(const Vector< T > &, const unsigned int)
Definition: f_matrix.h:424
F_TMFORALL_TM
#define F_TMFORALL_TM(op)
Definition: f_matrix.h:631
F_TMatrix::herm
F_TMatrix< T > & herm()
Definition: f_matrix.h:546
F_TMatrix::vec
T * vec
Definition: f_matrix.h:78
F_TMatrix::swap
F_TMatrix< T > & swap(F_TMatrix< T > &)
Definition: f_matrix.h:1688
F_TMatrix::TALIGN
T aligned_value_type TALIGN(MIN_ALIGN2)
Definition: f_matrix.h:85
F_Matrix::value_type
T value_type
Definition: f_matrix.h:1175
TVector
Temporary Base Class Idiom: Class TVector is used for temporary variables.
Definition: bvector.h:52
F_STDDEF_TTS
#define F_STDDEF_TTS(op)
Definition: f_matrix.h:1040
F_TMatrix::operator-
F_TMatrix< T > & operator-()
Definition: f_matrix.h:579
F_Matrix::F_Matrix
F_Matrix(const F_TMatrix< T > &tm)
Definition: f_matrix.h:1193
F_TMatrix::endvec
T * endvec
Definition: f_matrix.h:78
NAMESPACE_END
#define NAMESPACE_END
Definition: basics.h:316
F_TMatrix::operator*=
F_TSMatrix< T > operator*=(const T &)
Definition: f_matrix.h:673
F_TSMatrix::endvec
T * endvec
Definition: f_matrix.h:789
Vector
Definition: bvector.h:54
transpose
F_TMatrix< T > transpose(const F_TMatrix< T > &ftm)
Definition: f_matrix.h:1666
F_TMatrix::operator=
F_TMatrix< T > & operator=(const F_Matrix< T > &)
Definition: f_matrix.h:358
F_Matrix::set_col
void set_col(const Vector< T > &, const unsigned int)
Definition: f_matrix.h:1599
F_TMatrix::setval
T & setval(const unsigned int r, const unsigned int c)
Definition: f_matrix.h:208
BVector::vec
T * vec
Definition: bvector.h:73
F_TSMatrix::F_TSMatrix
F_TSMatrix(const F_TSMatrix< T > &ts)
Definition: f_matrix.h:823
ALIGN3
#define ALIGN3(v, i, x)
Definition: basics.h:435
TBCI::fabssqr
double fabssqr(const double a)
Definition: basics.h:1138
T
#define T
Definition: bdmatlib.cc:20
TBCICOMP
#define TBCICOMP(n, o, t, s)
Definition: basics.h:960
true
#define true
Definition: bool.h:24
F_TSMatrix::col
unsigned int col
Definition: f_matrix.h:787
F_TMatrix::resize
F_TMatrix< T > & resize(const unsigned int d)
Definition: f_matrix.h:212
F_Matrix::operator==
bool operator==(const F_Matrix< T > &m) const
Definition: f_matrix.h:1527
operator>>
STD__ istream & operator>>(STD__ istream &in, F_Matrix< T > &m)
Definition: f_matrix.h:1378
F_TMFORALL_TS
#define F_TMFORALL_TS(op)
Definition: f_matrix.h:642
F_Matrix::element_type
T element_type
Definition: f_matrix.h:1176
F_STDDEF_TSM
#define F_STDDEF_TSM(op)
Definition: f_matrix.h:979
F_TMatrix::~F_TMatrix
~F_TMatrix()
Definition: f_matrix.h:114
CSCMatrix::columns
unsigned int columns() const
Definition: cscmatrix.h:102
F_TSMatrix::mut
bool mut
Definition: f_matrix.h:793
F_TSMatrix::operator==
bool operator==(const F_Matrix< T > &)
Definition: f_matrix.h:1091
MATH__
#define MATH__
Definition: basics.h:332
rowvec
Definition: matrix.h:70
F_TMatrix::mat
T ** mat
Fortran storage layout: mat[col][row].
Definition: f_matrix.h:77
F_Matrix::herm
F_TMatrix< T > herm() const
Definition: f_matrix.h:1310
F_Matrix::F_Matrix
F_Matrix(const F_Matrix< T > &m)
Definition: f_matrix.h:1189
F_TSMatrix::operator()
T operator()(const unsigned r, const unsigned c)
Definition: f_matrix.h:827
F_STDDEF_TST
#define F_STDDEF_TST(op)
Definition: f_matrix.h:1026
colvec
Definition: matrix.h:70
F_STDDEF_TMM
#define F_STDDEF_TMM(op)
Definition: f_matrix.h:685
LIKELY
#define LIKELY(expr)
branch prediction note that we sometimes on purpose mark the unlikely possibility likely and vice ver...
Definition: basics.h:100
F_TMatrix::fabs
double fabs() const
Definition: f_matrix.h:774
F_STDDEF_TSTM
#define F_STDDEF_TSTM(op)
Definition: f_matrix.h:994
F_TSMatrix::eval
F_TSMatrix< T > & eval(F_TMatrix< T > *=0)
Definition: f_matrix.h:960
F_Matrix::operator()
T & operator()(const unsigned int, const unsigned int) const
Definition: f_matrix.h:1623
F_TMatrix::operator+=
F_TMatrix< T > & operator+=(F_TMatrix< T >)
Definition: f_matrix.h:638
F_TMatrix::set_row
void set_row(const Vector< T > &, const unsigned int)
Definition: f_matrix.h:441
rowcolvec
rowcolvec
Definition: matrix.h:70
Generated by   doxygen 1.8.6