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.61 2022/11/03 17:28:10 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; }
918  endvec = (T*)0; dim = 0;
919  if (mat) TBCIDELETE(T*,mat,col);
920 }
921 
922 template <typename T>
923 inline void F_TSMatrix<T>::detach (F_TMatrix<T>* tm)
924 {
925  if ((!tm && mut) || !dim) return;
926  if (!tm) //mut = false
927  {
928  vec = NEW(T ,dim); if (!vec) { dim = 0; mat = (T**)0; endvec = (T*)0; return; }
929  mat = NEW(T*,col); if (!mat) { TBCIDELETE(T,vec,dim); dim = 0; vec = (T*)0; endvec = (T*)0; return; }
930  for (REGISTER unsigned c = 0; c < row; c++) mat[c] = vec + c*row;
931  endvec = vec + dim;
932  }
933  else //tm = true, mut = ?
934  {
935  vec = tm->vec; mat = tm->mat; endvec = tm->endvec;
936  }
937  mut = true;
938 }
939 
940 template <typename T>
941 void F_TSMatrix<T>::clone (bool evl, F_TMatrix<T>* tm)
942 {
943  if ((mut && !tm)|| !dim) return;
944  T* oldv = vec; T** oldm = mat; bool omut = mut;
945  detach (tm);
946  if (evl && fac != (T)1)
947  {
948  for (REGISTER T *p1 = vec, *p2 = oldv; p1 < endvec; p1++, p2++)
949  *p1 = *p2 * fac;
950  fac = (T)1;
951  }
952  else TBCICOPY (vec, oldv, T, dim);
953  if (tm && omut)
954  { TBCIDELETE(T,oldv,dim); TBCIDELETE(T*,oldm,col); }
955 }
956 
957 
958 
959 //creates mutable evaluated F_TSMatrix (can be easily converted to F_TMatrix)
960 template <typename T>
961 inline F_TSMatrix<T>& F_TSMatrix<T>::eval (F_TMatrix<T>* tm)
962 {
963  if (mut && !tm)
964  {
965  if (fac != (T)1)
966  {
967  for (REGISTER T* ptr = vec; ptr < endvec; ptr++)
968  *ptr *= fac;
969  fac = (T)1;
970  }
971  }
972  else
973  clone (true, tm);
974  return (*this);
975 }
976 
977 
978 // Arithmetics
979 
980 #define F_STDDEF_TSM(op) \
981 template <typename T> \
982 inline F_TMatrix<T> F_TSMatrix<T>::operator op (const F_Matrix<T>& m) \
983 { \
984  BCHK(row != m.row || col != m.col, MatErr, Op op on diff size Mats, m.row, F_TMatrix<T> (*this)); \
985  F_TSMatrix<T> ts (*this); ts.detach (); \
986  for (REGISTER T *p1 = ts.vec, *p2 = vec, *p3 = m.vec; p2 < endvec; p1++, p2++, p3++) \
987  *p1 = *p2 * fac op *p3; \
988  ts.fac = (T)1; return F_TMatrix<T> (ts); \
989 }
990 
991 F_STDDEF_TSM(+)
992 F_STDDEF_TSM(-)
993 
994 
995 #define F_STDDEF_TSTM(op) \
996 template <typename T> \
997 inline F_TMatrix<T> F_TSMatrix<T>::operator op (const F_TMatrix<T>& tm) \
998 { \
999  BCHK(row != tm.row || col != tm.col, MatErr, Op op on diff size Mats, tm.row, tm); \
1000  for (REGISTER T *p1 = tm.vec, *p2 = vec; p2 < endvec; p1++, p2++) \
1001  *p1 = *p2 * fac op *p1; \
1002  destroy (); return tm; \
1003 }
1004 
1005 F_STDDEF_TSTM(+)
1006 F_STDDEF_TSTM(-)
1007 
1008 
1009 #define F_STDDEF_TSTS(op) \
1010 template <typename T> \
1011 inline F_TMatrix<T> F_TSMatrix<T>::operator op (F_TSMatrix<T> ts) \
1012 { \
1013  BCHK(row != ts.row || col != ts.col, MatErr, Op op on diff size Mats, ts.row, F_TMatrix<T> (*this)); \
1014  F_TSMatrix<T> tm; \
1015  if (!mut && ts.mut) tm = ts; \
1016  else { tm = *this; tm.detach (); } \
1017  for (REGISTER T *p1 = tm.vec, *p2 = vec, *p3 = ts.vec; p2 < endvec; p1++, p2++, p3++) \
1018  *p1 = *p2 * fac op *p3 * ts.fac; \
1019  if (!mut && ts.mut) destroy (); else ts.destroy (); \
1020  tm.fac = (T)1; return F_TMatrix<T> (tm); \
1021 }
1022 
1023 F_STDDEF_TSTS(+)
1024 F_STDDEF_TSTS(-)
1025 
1026 
1027 #define F_STDDEF_TST(op) \
1028 template <typename T> \
1029 inline F_TMatrix<T> F_TSMatrix<T>::operator op (const T& a) \
1030 { \
1031  F_TSMatrix<T> ts (*this); ts.detach (); \
1032  for (REGISTER T *p1 = ts.vec, *p2 = vec; p2 < endvec; p1++, p2++) \
1033  *p1 = *p2 * fac op a; \
1034  ts.fac = (T)1; return F_TMatrix<T> (ts); \
1035 }
1036 
1037 F_STDDEF_TST(+)
1038 F_STDDEF_TST(-)
1039 
1040 // operate on diag only?
1041 #define F_STDDEF_TTS(op) \
1042 INST(template <typename T> class F_TSMatrix friend F_TMatrix<T> operator op (const T&, F_TSMatrix<T>);) \
1043 template <typename T> \
1044 inline F_TMatrix<T> operator op (const T& a, F_TSMatrix<T> ts) \
1045 { \
1046  F_TSMatrix<T> tm (ts); tm.detach (); \
1047  for (REGISTER T *p1 = tm._VEC, *p2 = ts._VEC; p1 < tm._ENDVEC; p1++, p2++) \
1048  *p1 = a op *p2 * ts._FAC; \
1049  tm._FAC = (T)1; return F_TMatrix<T> (tm); \
1050 }
1051 
1052 F_STDDEF_TTS(+)
1053 F_STDDEF_TTS(-)
1054 
1055 
1056 template <typename T>
1057 F_TMatrix<T> F_TSMatrix<T>::operator * (const F_Matrix<T>& m)
1058 {
1059  F_TMatrix<T> ftm (row, m.col);
1060  BCHK(col != m.row, MatErr, Illegal sizes in F_Matrix multiplication, 0, (destroy(), ftm));
1061  REGISTER T tmp ALIGN(MIN_ALIGN);
1062  for (unsigned int r=0; r<row; r++) {
1063  for (unsigned int c=0; c<m.col; c++) {
1064  tmp = mat[0][r] * m(0, c);
1065  for (REGISTER unsigned int l=1; l<col; l++)
1066  tmp += mat[l][r] * m(l, c);
1067  ftm.setval(tmp*fac, r, c);
1068  }
1069  }
1070  destroy (); return ftm;
1071 }
1072 
1073 
1074 template <typename T>
1075 TVector<T> F_TSMatrix<T>::operator * (const Vector<T>& v)
1076 {
1077  TVector<T> tv (row);
1078  BCHK (v.size() != col, MatErr, multiply F_Matrix w/ Vector: wrong config, v.size(), (destroy(), tv));
1079  /* FIXME: We could call do_mat_vec_transmult<>() here */
1080  for (unsigned int r=0; r<row; r++)
1081  {
1082  ALIGN3(REGISTER T el, mat[0][r] * v(0), MIN_ALIGN2);
1083  for (REGISTER unsigned int c=1; c<col; c++)
1084  el += mat[c][r] * v(c);
1085  tv.set (el * fac, r);
1086  }
1087  destroy (); return tv;
1088 }
1089 
1090 
1091 template <typename T>
1092 bool F_TSMatrix<T>::operator == (const F_Matrix<T>& m)
1093 {
1094  if (row != m.row || col != m.col) { destroy (); return false; }
1095  if (dim == 0) { /*destroy();*/ return true; }
1096  if (fac == (T)1)
1097  {
1098  if (vec == m.vec) return true;
1099  if (TBCICOMP (vec, m.vec, T, dim)) { destroy (); return false; }
1100  }
1101  else
1102  {
1103  if (vec == m.vec) return false;
1104  for (REGISTER T *p1 = vec, *p2 = m.vec; p1 < endvec; p1++, p2++)
1105  if (*p1 * fac != *p2) { destroy (); return false; }
1106  }
1107  destroy (); return true;
1108 }
1109 
1110 
1111 template <typename T>
1112 bool F_TSMatrix<T>::operator == (F_TSMatrix<T> ts)
1113 {
1114  if (row != ts.row || col != ts.col) { destroy (); ts.destroy (); return false; }
1115  if (dim == 0) {/*destroy(); ts.destroy();*/ return true; }
1116  if (fac == ts.fac)
1117  {
1118  if (vec == ts.vec) { destroy (); return true; }
1119  if (TBCICOMP (vec, ts.vec, T, dim)) { destroy (); ts.destroy (); return false; }
1120  }
1121  else
1122  {
1123  if (vec == ts.vec) { destroy (); return false; }
1124  for (REGISTER T *p1 = vec, *p2 = ts.vec; p1 < endvec; p1++, p2++)
1125  if (*p1 * fac != *p2) { destroy (); ts.destroy (); return false; }
1126  }
1127  destroy (); ts.destroy (); return true;
1128 }
1129 
1130 template <typename T>
1131 F_TSMatrix<T>::operator TMatrix<T>()
1132 {
1133  TMatrix<T> tm(this->row, this->col);
1134  for (unsigned int c = 0; c < this->col; ++c)
1135  for (unsigned int r = 0; r < this->row; ++r)
1136  tm.setval(this->get(r, c), r, c);
1137  destroy();
1138  return tm;
1139 }
1140 
1141 
1142 template <typename T>
1143 double F_TSMatrix<T>::fabs ()
1144 {
1145  REGISTER T* ptr = vec;
1146  double REGISTER rv ALIGN(8) = fabssqr (*ptr++);
1147  for (; ptr < endvec; ptr++)
1148  rv += fabssqr (*ptr);
1149  rv *= fabssqr (fac);
1150  destroy (); return MATH__ sqrt (rv);
1151 }
1152 
1153 
1154 INST(template <typename T> class F_TSMatrix friend F_TSMatrix<T> operator * (const T&, F_TSMatrix<T>);)
1155 template <typename T>
1156 inline F_TSMatrix<T> operator * (const T& f, F_TSMatrix<T> ts)
1157 { ts._FAC = f * ts._FAC; return ts; }
1158 
1159 //---------------------------------------------------------------
1160 // Real, user referable, matrix
1161 //---------------------------------------------------------------
1170 template <typename T>
1171 class F_Matrix : public F_TMatrix<T>
1172 {
1173 
1174  public:
1175 
1176  typedef T value_type;
1177  typedef T element_type;
1178  typedef T aligned_value_type TALIGN(MIN_ALIGN2);
1179 
1180  friend class F_TSMatrix<T>;
1181 
1182  // constructor, destructor
1183  explicit F_Matrix (const unsigned d=0) : F_TMatrix<T> (d) {}
1184  F_Matrix (const unsigned r, const unsigned c) : F_TMatrix<T> (r, c) {}
1185  F_Matrix (const T& v, const unsigned r, const unsigned c) : F_TMatrix<T> (v, r, c) {}
1186  //F_Matrix (const Vector<T>& v, const enum rowcolvec r = colvec)
1187  // : F_TMatrix<T> (v, r) {}
1188 
1189  //copy
1190  F_Matrix (const F_Matrix<T>& m) : F_TMatrix<T> (m) {}
1191  F_Matrix (const CSCMatrix<T>& m);
1192 
1193  //alias
1194  F_Matrix (const F_TMatrix<T>& tm) : F_TMatrix<T> (tm) {}
1195  F_Matrix (F_TSMatrix<T>& ts) : F_TMatrix<T> (ts) {}
1196 
1197  ~F_Matrix () { this->destroy (); }
1198 
1199  // Conversion
1200  F_Matrix (const Matrix<T>& m);
1201  operator TMatrix<T>() const;
1202 
1203  // allow instantiation (Matrix_Sig)
1204  /*virtual*/ static const char* mat_info () { return ("F_Matrix"); }
1205 
1206  // Access
1207  T& operator () (const unsigned int, const unsigned int) const;
1208  //T operator () (const unsigned, const unsigned) const;
1209  T* get_fortran_matrix() const { return this->vec; }
1210  const T& get (const unsigned r, const unsigned c) const
1211  { return this->mat[c][r]; }
1212 
1213  unsigned int columns () const { return this->row; }
1214  unsigned int rows () const { return this->col; }
1215  unsigned long size () const { return this->dim; }
1216  //TVector<T> operator () (const unsigned int) const;
1217  TVector<T> get_row (const unsigned int) const;
1218  TVector<T> get_col (const unsigned int) const;
1219  void set_row (const Vector<T>&, const unsigned int);
1220  void set_col (const Vector<T>&, const unsigned int);
1221  void setval(const T val, const unsigned int r, const unsigned int c)
1222  { this->operator() (r, c) = val; }
1223 
1224  // basics
1225  F_Matrix<T>& operator = (const F_Matrix<T>& m)
1226  { this->F_TMatrix<T>::operator = (m); return *this; }
1227  F_Matrix<T>& operator = (const F_TMatrix<T>& tm)
1228  { this->F_TMatrix<T>::operator = (tm); return *this; }
1229  F_Matrix<T>& operator = (F_TSMatrix<T> ts)
1230  { F_TMatrix<T>(this->F_TMatrix<T>::operator = (ts)); return *this; }
1231  F_Matrix<T>& operator = (const T& v)
1232  { this->F_TMatrix<T>::operator = (v); return *this; }
1233 
1234  F_Matrix<T>& resize (const unsigned r, const unsigned c)
1235  { this->F_TMatrix<T>::resize (r, c); return *this; }
1236  F_Matrix<T>& resize (const unsigned d)
1237  { return this->resize (d, d); }
1238  F_Matrix<T>& resize (const T& v, const unsigned r, const unsigned c)
1239  { this->F_TMatrix<T>::resize (v, r, c); return *this; }
1240  F_Matrix<T>& fill (const T& v = 0)
1241  { this->F_TMatrix<T>::fill (v); return *this; }
1242  F_Matrix<T>& setunit (const T& f = 1)
1243  { this->F_TMatrix<T>::setunit (f); return *this; }
1244  F_Matrix<T>& clear ()
1245  { this->F_TMatrix<T>::clear (); return *this; }
1246 
1247  // assignment ops (do we really need em ?)
1248  // Seems we have to make sure, that a F_Matrix and no F_TMatrix
1249  // is returned; otherwise there's danger to be destroyed
1250  F_Matrix<T>& operator += (const F_Matrix<T>& a)
1251  { this->F_TMatrix<T>::operator += (a); return *this; }
1252  F_Matrix<T>& operator -= (const F_Matrix<T>& a)
1253  { this->F_TMatrix<T>::operator -= (a); return *this; }
1254  F_Matrix<T>& operator += (F_TMatrix<T> a)
1255  { this->F_TMatrix<T>::operator += (a); return *this; }
1256  F_Matrix<T>& operator -= (F_TMatrix<T> a)
1257  { this->F_TMatrix<T>::operator -= (a); return *this; }
1258  F_Matrix<T>& operator += (const T& a)
1259  { this->F_TMatrix<T>::operator += (a); return *this; }
1260  F_Matrix<T>& operator -= (const T& a)
1261  { this->F_TMatrix<T>::operator -= (a); return *this; }
1262  F_Matrix<T>& operator *= (const T& a)
1263  { this->F_TMatrix<T>::operator *= (a); return *this; }
1264  F_Matrix<T>& operator /= (const T& a)
1265  { this->F_TMatrix<T>::operator /= (a); return *this; }
1266 
1267 
1268  // basic operations
1269  F_TMatrix<T> operator - () const { return -(F_TMatrix<T> (*this)); }
1270 
1271  F_TMatrix<T> operator + (const F_Matrix<T>&) const;
1272  F_TMatrix<T> operator - (const F_Matrix<T>&) const;
1273  F_TMatrix<T> operator * (const F_Matrix<T>&) const;
1274 
1275  F_TMatrix<T> operator + (F_TMatrix<T>) const;
1276  F_TMatrix<T> operator - (F_TMatrix<T>) const;
1277  F_TMatrix<T> operator * (F_TMatrix<T>) const;
1278 
1279  F_TMatrix<T> operator + (F_TSMatrix<T>) const;
1280  F_TMatrix<T> operator - (F_TSMatrix<T>) const;
1281  F_TMatrix<T> operator * (F_TSMatrix<T>) const;
1282 
1283  F_TMatrix<T> operator + (const T&) const;
1284  F_TMatrix<T> operator - (const T&) const;
1285  F_TSMatrix<T> operator * (const T&) const;
1286  F_TSMatrix<T> operator / (const T&) const;
1287 
1288  TVector<T> operator * (const Vector<T>&) const;
1289  TVector<T> operator * (TVector<T>&) const;
1290 
1291  bool operator == (const F_Matrix<T>& m) const;
1292  bool operator != (const F_Matrix<T>& m) const
1293  { return !(*this == m); }
1294 
1295  bool operator == (F_TMatrix<T>& tm) const
1296  { F_Matrix<T> m (tm); return (*this == m); }
1297  bool operator != (F_TMatrix<T>& tm) const
1298  { return !(*this == tm); }
1299 
1300  bool operator == (F_TSMatrix<T>&) const;
1301  bool operator != (F_TSMatrix<T>& ts) const
1302  { return !(*this == ts); }
1303 
1304  // io-streams
1305 
1306  friend STD__ ostream& operator << FGD (STD__ ostream&, const F_Matrix<T>&);
1307  friend STD__ istream& operator >> FGD (STD__ istream&, F_Matrix<T>&);
1308 
1309  // misc
1310  // returns the the hermitian adjungated (Transposed and conjugate complex)
1311  F_TMatrix<T> herm () const { return (F_TMatrix<T> (*this)).herm(); }
1312  F_TMatrix<T> conj () const { return (F_TMatrix<T> (*this)).conj(); }
1313  F_TMatrix<T> trans() const { return (F_TMatrix<T> (*this)).trans(); }
1314 
1315  double fabs () const;
1316 };
1317 
1318 
1319 // definitions
1320 
1321 template <typename T>
1322 inline F_Matrix<T>::F_Matrix (const CSCMatrix<T>& m)
1323 : F_TMatrix<T>(m.rows(), m.columns())
1324 {
1325  for (unsigned int r=0; r<m.rows(); r++)
1326  for (unsigned int c=0; c<m.columns(); c++)
1327  operator()(r,c) = m(r,c);
1328 }
1329 
1330 template <typename T>
1331 F_Matrix<T>::F_Matrix(const Matrix<T>& m)
1332 : F_TMatrix<T>(m.rows(), m.columns())
1333 {
1334  for (unsigned int r = 0; r < this->row; ++r)
1335  for (unsigned int c = 0; c < this->col; ++c)
1336  setval(m(r, c), r, c);
1337 }
1338 
1339 template <typename T>
1340 F_Matrix<T>::operator TMatrix<T>() const
1341 {
1342  TMatrix<T> tm(this->row, this->col);
1343  for (unsigned int c = 0; c < this->col; ++c)
1344  for (unsigned int r = 0; r < this->row; ++r)
1345  tm.setval((*this)(r, c), r, c);
1346  return tm;
1347 }
1348 
1349 
1350 template <typename T>
1351 STD__ ostream& operator << (STD__ ostream& os, const F_Matrix<T>& m)
1352 {
1353  for (unsigned int r=0; r<m.row; r++) {
1354  for (unsigned int c=0; c<m.col; c++)
1355  os << m(r,c) << " ";
1356  os << "\n";
1357  }
1358  return os;
1359 }
1360 
1361 INST(template <typename T> class F_TMatrix friend STD__ ostream& operator << (STD__ ostream& os, F_TMatrix<T> tm);)
1362 template <typename T>
1363 STD__ ostream& operator << (STD__ ostream& os, F_TMatrix<T> tm)
1364 {
1365  F_Matrix<T> fm(tm);
1366  return os << fm;
1367 }
1368 
1369 #if 1
1370 INST(template <typename T> class F_TSMatrix friend STD__ ostream& operator << (STD__ ostream& os, F_TSMatrix<T> ts);)
1371 template <typename T>
1372 STD__ ostream& operator << (STD__ ostream& os, F_TSMatrix<T> ts)
1373 {
1374  return os << F_Matrix<T>(ts);
1375 }
1376 #endif
1377 
1378 template <typename T>
1379 STD__ istream& operator >> (STD__ istream& in, F_Matrix<T>& m)
1380 {
1381  Vector<T> buf(m.row);
1382 
1383  for (unsigned int i=0; i<m.col; i++)
1384  {
1385  in >> buf;
1386  m.set_col (buf, i);
1387  }
1388  return in;
1389 }
1390 
1391 
1392 /************** Arithmetics ***************/
1393 
1394 //change these to operate on diag only ?
1395 #define F_MSTDDEF_T(op) \
1396 template <typename T> \
1397 inline F_TMatrix<T> F_Matrix<T>::operator op (const T& a) const \
1398 { \
1399  F_TMatrix<T> t (this->row, this->col); \
1400  for (REGISTER T *p1=t.vec, *p2=this->vec; p2<this->endvec; ++p1, ++p2) \
1401  *p1 = *p2 op a; \
1402  return t; \
1403 }
1404 
1405 F_MSTDDEF_T(+)
1406 F_MSTDDEF_T(-)
1407 
1408 template <typename T>
1409 inline F_TSMatrix<T> F_Matrix<T>::operator * (const T& a) const
1410 { return F_TSMatrix<T> (*this, a); }
1411 
1412 template <typename T>
1413 inline F_TSMatrix<T> F_Matrix<T>::operator / (const T& a) const
1414 {
1415  BCHK (a==(T)0, MatErr, Divide Mat by 0, 0, *this);
1416  return F_TSMatrix<T> (*this, (T)1/a);
1417 }
1418 
1419 #define F_MSTDDEF_M(op) \
1420 template <typename T> \
1421 inline F_TMatrix<T> F_Matrix<T>::operator op (const F_Matrix<T>& a) const \
1422 { \
1423  F_TMatrix<T> t (this->row, this->col); \
1424  BCHK(this->dim != a.dim, MatErr, Operator op on diff size matrices, a.dim, t); \
1425  for (REGISTER T *p1=t.vec, *p2=this->vec, *p3=a.vec; p2<this->endvec; ++p1, ++p2, ++p3) \
1426  *p1 = *p2 op *p3; \
1427  return t; \
1428 }
1429 
1430 F_MSTDDEF_M(+)
1431 F_MSTDDEF_M(-)
1432 
1433 
1434 #define F_MSTDDEF_TM(op) \
1435 template <typename T> \
1436 inline F_TMatrix<T> F_Matrix<T>::operator op (F_TMatrix<T> a) const \
1437 { \
1438  BCHK(this->dim != a.dim, MatErr, Applying op on diff. size matrices, a.dim, a); \
1439  for (REGISTER T *p1=a.vec, *p2=this->vec; p2<this->endvec; ++p1, ++p2) \
1440  *p1 = *p2 op *p1; \
1441  return a; \
1442 }
1443 
1444 F_MSTDDEF_TM(+)
1445 F_MSTDDEF_TM(-)
1446 
1447 #define F_MSTDDEF_TS(op) \
1448 template <typename T> \
1449 inline F_TMatrix<T> F_Matrix<T>::operator op (F_TSMatrix<T> ts) const \
1450 { \
1451  BCHK(this->dim != ts.dim, MatErr, Applying op on diff. size matrices, ts.dim, (ts.destroy(),F_TMatrix<T> (*this)) ); \
1452  F_TSMatrix<T> tm (ts); tm.detach (); \
1453  for (REGISTER T *p1=tm.vec, *p2=this->vec, *p3=ts.vec; p2<this->endvec; ++p1, ++p2, ++p3) \
1454  *p1 = *p2 op *p3 * ts.fac; \
1455  tm.fac = (T)1; return F_TMatrix<T> (tm); \
1456 }
1457 
1458 F_MSTDDEF_TS(+)
1459 F_MSTDDEF_TS(-)
1460 
1461 
1462 // F_Matrix multiplication
1463 template <typename T>
1464 F_TMatrix<T> F_Matrix<T>::operator * (const F_Matrix<T>& a) const
1465 {
1466  F_TMatrix<T> ftm(this->row, a.col);
1467  BCHK(this->col != a.row, MatErr, Illegal sizes in F_Matrix multiplication, 0, ftm);
1468  for (unsigned int r=0; r<this->row; ++r) {
1469  for (REGISTER unsigned int c=0; c<a.col; ++c) {
1470  ALIGN3(REGISTER T tmp, (*this)(r, 0) * a(0, c), MIN_ALIGN2);
1471  for (REGISTER unsigned int l=1; l<this->col; ++l)
1472  tmp += (*this)(r, l) * a(l, c);
1473  ftm.setval(tmp, r, c);
1474  }
1475  }
1476  return ftm;
1477 }
1478 
1479 template <typename T>
1480 F_TMatrix<T> F_Matrix<T>::operator * (F_TSMatrix<T> a) const
1481 {
1482  F_TMatrix<T> ftm(this->row, a.col);
1483  BCHK(this->col != a.row, MatErr, Illegal sizes in F_Matrix multiplication, 0, ftm);
1484  REGISTER T tmp ALIGN(MIN_ALIGN2);
1485  for (unsigned int r=0; r<this->row; ++r) {
1486  for (REGISTER unsigned int c=0; c<a.col; ++c) {
1487  tmp = (*this)(c, 0) * a.mat[c][0];
1488  for (REGISTER unsigned int l=1; l<this->col; ++l)
1489  tmp += (*this)(r, l) * a.mat[c][l];
1490  ftm.setval(tmp*a.fac, r, c);
1491  }
1492  }
1493  a.destroy (); return ftm;
1494 }
1495 
1496 
1497 template <typename T>
1498 inline F_TMatrix<T> F_Matrix<T>::operator * (F_TMatrix<T> a) const
1499 {
1500  F_TMatrix<T> ftm (this->operator * (F_Matrix<T> (a)));
1501  return ftm;
1502 }
1503 
1504 template <typename T>
1505 inline F_TMatrix<T> F_TMatrix<T>::operator * (const F_Matrix<T>& a)
1506 {
1507  F_Matrix<T> m (*this);
1508  return (m * a);
1509 }
1510 
1511 
1512 template <typename T>
1513 inline F_TMatrix<T> F_TMatrix<T>::operator * (F_TMatrix<T> a)
1514 {
1515  F_Matrix<T> m (*this);
1516  return (m * (F_Matrix<T> (a)));
1517 }
1518 
1519 template <typename T>
1520 inline F_TMatrix<T> F_TMatrix<T>::operator * (F_TSMatrix<T> a)
1521 {
1522  F_Matrix<T> m (*this);
1523  return (m * a);
1524 }
1525 
1526 
1527 template <typename T>
1528 bool F_Matrix<T>::operator == (const F_Matrix<T>& m) const
1529 {
1530  if (this->row != m.row || this->col != m.col)
1531  return false;
1532  if (this->vec == m.vec || this->dim == 0)
1533  return true;
1534  if (TBCICOMP (this->vec, m.vec, T, this->dim))
1535  return false;
1536  else
1537  return true;
1538 }
1539 
1540 template <typename T>
1541 bool F_Matrix<T>::operator == (F_TSMatrix<T>& ts) const
1542 {
1543  if (this->row != ts.row || this->col != ts.col)
1544  return (ts.destroy(), false);
1545  if (this->vec == 0) {
1546  /*ts.destroy();*/ return true;
1547  }
1548  if (ts.fac == (T)1) {
1549  if (this->vec == ts.vec)
1550  return true;
1551  if (TBCICOMP(this->vec, ts.vec, T, this->dim)) {
1552  ts.destroy (); return false;
1553  }
1554  } else {
1555  if (this->vec == ts.vec)
1556  return false;
1557  for (REGISTER T *p1 = this->vec, *p2 = ts.vec; p1 < this->endvec; ++p1, ++p2)
1558  if (*p1 != *p2 * ts.fac) {
1559  ts.destroy (); return false;
1560  }
1561  }
1562  ts.destroy ();
1563  return true;
1564 }
1565 
1566 
1567 template<typename T>
1568 TVector<T> F_Matrix<T>::operator * (const Vector<T>& v) const
1569 {
1570  TVector<T> tv (this->row);
1571  for (unsigned int r=0; r<this->row; ++r) {
1572  REGISTER T el = 0;
1573  for (REGISTER unsigned int c=0; c<this->col; ++c)
1574  el += this->operator() (r, c) * v(c);
1575  tv.set (el, r);
1576  }
1577  return tv;
1578 }
1579 
1580 template<typename T>
1581 inline TVector<T> F_Matrix<T>::operator * (TVector<T>& c) const
1582 {
1583  Vector<T> cn (c);
1584  return (*this * cn);
1585 }
1586 
1587 #if 0
1588 template <typename T>
1589 inline TVector<T> F_Matrix<T>::operator () (const unsigned int i) const
1590 {
1591  TVector<T> v(this->row);
1592  BCHK(i < 0 || i >= this->col, MatErr, Illegal row index in get_row, i, v);
1593  for (REGISTER unsigned int j=0; j<this->row; ++j)
1594  v.vec[j] = this->mat[j][i];
1595  return v;
1596 }
1597 #endif
1598 
1599 template <typename T>
1600 inline void F_Matrix<T>::set_col (const Vector<T>& v, const unsigned int c)
1601 {
1602  F_TMatrix<T>::set_col(v, c);
1603 }
1604 
1605 template <typename T>
1606 inline TVector<T> F_Matrix<T>::get_col (const unsigned int c) const
1607 {
1608  return F_TMatrix<T>::get_col(c);
1609 }
1610 
1611 template <typename T>
1612 inline TVector<T> F_Matrix<T>::get_row (const unsigned int r) const
1613 {
1614  return F_TMatrix<T>::get_row(r);
1615 }
1616 
1617 template <typename T>
1618 inline void F_Matrix<T>::set_row (const Vector<T>& v, const unsigned int r)
1619 {
1620  F_TMatrix<T>::set_row(v, r);
1621 }
1622 
1623 template <typename T>
1624 inline T& F_Matrix<T>::operator () (const unsigned int r, const unsigned int c) const
1625 {
1626  BCHK(r>=this->row, MatErr, illegal row index, r, this->mat[0][0]);
1627  BCHK(c>=this->col, MatErr, Illegal col index, c, this->mat[0][0]);
1628  return this->mat[c][r];
1629 }
1630 
1631 /*
1632 template <typename T>
1633 inline T F_Matrix<T>::operator () (const unsigned r, const unsigned c) const
1634 {
1635  BCHK(r<0 || r>=row, MatErr, Illegal row index, r, mat[0][0]);
1636  BCHK(c<0 || c>=col, MatErr, illegal col index, c, mat[0][0]);
1637  return mat[c][r];
1638 }
1639  */
1640 
1641 template <typename T>
1642 double F_Matrix<T>::fabs () const
1643 {
1644  REGISTER double rv = 0.0;
1645  for (REGISTER T *ptr = this->vec; ptr < this->endvec; ++ptr)
1646  rv += fabssqr (*ptr);
1647  return MATH__ sqrt (rv);
1648 }
1649 
1650 
1651 template <typename T>
1652 F_TMatrix<T> F_TMatrix<T>::transposed_copy() const
1653 {
1654  F_TMatrix<T> trm(this->col, this->row);
1655  for (unsigned r = 0; r < this->row; ++r)
1656  for (unsigned c = 0; c < this->col; ++c)
1657  trm.mat[c][r] = this->mat[r][c];
1658  // memleak if it's really a F_TMatrix :-(
1659  //this->destroy();
1660  //F_Matrix<T> fm(*this);
1661  return trm;
1662 }
1663 
1664 
1665 INST(template <typename T> class F_TMatrix friend F_TMatrix<T> transpose(const F_TMatrix<T>& tm);)
1666 template <typename T>
1667 inline F_TMatrix<T> transpose(const F_TMatrix<T>& ftm)
1668 {
1669  return ftm.transposed_copy();
1670 }
1671 
1672 template <typename T>
1673 F_TMatrix<T>& F_TMatrix<T>::transpose()
1674 {
1675  if (this->col != this->row) {
1676  // for non quadratic shape, doing inplace swap ops is too hard
1677  F_Matrix<T> transmat(this->transposed_copy());
1678  return this->swap(transmat);
1679  } else {
1680  // do in place swap ops
1681  for (unsigned r = 1; r < this->row; ++r)
1682  for (unsigned c = 0; c < r; ++c)
1683  SWAP(this->mat[r][c], this->mat[c][r]);
1684  return *this;
1685  }
1686 }
1687 
1688 template <typename T>
1689 inline F_TMatrix<T>& F_TMatrix<T>::swap (F_TMatrix<T>& m)
1690 {
1691  //BCHK(this->row != m.row || this->col != m.col, MatErr, Swapping Matrices w/ different layout, m.row, *this);
1692  SWAP(this->dim, m.dim);
1693  SWAP(this->row, m.row);
1694  //SWAP(this->col, m.col);
1695  unsigned int tcol = this->col; this->col = m.col; m.col = tcol;
1696  SWAP(this->mat, m.mat);
1697  SWAP(this->vec, m.vec); SWAP(this->endvec, m.endvec);
1698  return *this;
1699 }
1700 
1701 #undef _DIM
1702 #undef _VEC
1703 #undef _ENDVEC
1704 #undef _ROW
1705 #undef _COL
1706 #undef _FAC
1707 
1708 NAMESPACE_END
1709 
1710 NAMESPACE_CSTD
1711 
1712 INST(template <typename T> class TBCI__ F_TMatrix friend double fabs (const TBCI__ F_TMatrix<T>&);)
1713 template <typename T>
1714 inline double fabs (const TBCI__ F_TMatrix<T>& ftm)
1715 { return ftm.fabs (); }
1716 
1717 INST(template <typename T> class TBCI__ F_Matrix friend double fabs (const TBCI__ F_Matrix<T>&);)
1718 template <typename T>
1719 inline double fabs (const TBCI__ F_Matrix<T>& fm)
1720 { return fm.fabs (); }
1721 
1722 INST(template <typename T> class TBCI__ F_TSMatrix friend double fabs (TBCI__ F_TSMatrix<T>&);)
1723 template <typename T>
1724 inline double fabs (/*typename*/ TBCI__ F_TSMatrix<T>& ftsm)
1725 { return ftsm.fabs (); }
1726 
1727 NAMESPACE_CSTD_END
1728 
1729 #endif /* TBCI_F_MATRIX_H */
TBCICOPY
#define TBCICOPY(n, o, t, s)
Definition: basics.h:894
F_TSMatrix::~F_TSMatrix
~F_TSMatrix()
Definition: f_matrix.h:815
transpose
BdMatrix< T > transpose(BdMatrix< T > &mat)
Definition: band_matrix.h:1393
F_TMatrix::transposed_copy
F_TMatrix< T > transposed_copy() const
Inefficient! Use transMult if possible.
Definition: f_matrix.h:1652
BCHKNR
#define BCHKNR(cond, exc, txt, ind)
Definition: basics.h:586
x
const Vector< T > const Vector< T > & x
Definition: LM_fit.h:97
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:421
ALIGN
#define ALIGN(x)
Definition: basics.h:444
fac
tm fac
Definition: f_matrix.h:1052
F_Matrix::setval
void setval(const T val, const unsigned int r, const unsigned int c)
Definition: f_matrix.h:1221
F_Matrix::fabs
double fabs() const
Definition: f_matrix.h:1642
F_Matrix::operator+
F_TMatrix< T > operator+(const F_Matrix< T > &) const
Definition: f_matrix.h:1430
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:1606
F_Matrix::columns
unsigned int columns() const
Definition: f_matrix.h:1213
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
operator>>
STD__ istream & operator>>(STD__ istream &istr, BdMatrix< T > &mat)
Definition: band_matrix.h:2739
CSCMatrix::irow
unsigned int * irow
Definition: cscmatrix.h:217
REGISTER
#define REGISTER
Definition: basics.h:108
F_Matrix::operator+=
F_Matrix< T > & operator+=(const F_Matrix< T > &a)
Definition: f_matrix.h:1250
c
return c
Definition: f_matrix.h:760
std::fabs
double fabs(const int a)
Definition: basics.h:1215
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:655
FRIEND_TBCI__
#define FRIEND_TBCI__
Definition: basics.h:334
F_Matrix::F_Matrix
F_Matrix(const unsigned d=0)
Definition: f_matrix.h:1183
F_Matrix::get_row
TVector< T > get_row(const unsigned int) const
Definition: f_matrix.h:1612
NAMESPACE_TBCI
#define NAMESPACE_TBCI
Definition: basics.h:317
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:1413
F_TMatrix::trace
T trace() const
Definition: f_matrix.h:517
F_MSTDDEF_T
#define F_MSTDDEF_T(op)
Definition: f_matrix.h:1395
F_Matrix::operator!=
bool operator!=(const F_Matrix< T > &m) const
Definition: f_matrix.h:1292
F_TSMatrix::detach
void detach(F_TMatrix< T > *=0)
Definition: f_matrix.h:923
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:424
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:1234
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:1618
F_TMatrix::rows
unsigned int rows() const
Definition: f_matrix.h:199
BCHK
#define BCHK(cond, exc, txt, ind, rtval)
Definition: basics.h:575
F_Matrix::trans
F_TMatrix< T > trans() const
Definition: f_matrix.h:1313
F_Matrix::get_fortran_matrix
T * get_fortran_matrix() const
Definition: f_matrix.h:1209
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:1225
F_TMatrix::element_type
T element_type
Definition: f_matrix.h:84
NAMESPACE_CSTD_END
#define NAMESPACE_CSTD_END
Definition: basics.h:325
F_Matrix::setunit
F_Matrix< T > & setunit(const T &f=1)
Definition: f_matrix.h:1242
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:1204
TBCIFILL
#define TBCIFILL(n, v, t, s)
Definition: basics.h:910
F_TSMatrix::clone
void clone(bool=false, F_TMatrix< T > *=0)
Definition: f_matrix.h:941
F_Matrix::TALIGN
T aligned_value_type TALIGN(MIN_ALIGN2)
Definition: f_matrix.h:1178
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:1464
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
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:1238
F_Matrix::operator/=
F_Matrix< T > & operator/=(const T &a)
Definition: f_matrix.h:1264
F_TSMatrix::value_type
T value_type
Definition: f_matrix.h:806
F_Matrix::clear
F_Matrix< T > & clear()
Definition: f_matrix.h:1244
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:1195
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:991
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:1197
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:1185
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:1216
for
for(REGISTER T *p1=c.vec,*p2=b.vec;p1< c.endvec;p1++, p2++)*p1
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:1240
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:1009
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:1236
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:1312
Tensor
Tensor class including arithmetics.
Definition: f_matrix.h:53
ts
F_TSMatrix< T > ts
Definition: f_matrix.h:1052
TBCI__
#define TBCI__
Definition: basics.h:332
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:1447
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:1262
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_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:1214
F_TMatrix::value_type
T value_type
Definition: f_matrix.h:83
NAMESPACE_CSTD
#define NAMESPACE_CSTD
Definition: basics.h:319
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:813
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:1252
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:1184
FGD
#define FGD
Definition: basics.h:144
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:1673
F_Matrix::operator-
F_TMatrix< T > operator-() const
Definition: f_matrix.h:1269
INST
#define INST(x)
Definition: basics.h:238
CSCMatrix::pcol
unsigned int * pcol
Definition: cscmatrix.h:216
F_MSTDDEF_TM
#define F_MSTDDEF_TM(op)
Definition: f_matrix.h:1434
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:1419
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:1143
F_Matrix::size
unsigned long size() const
Definition: f_matrix.h:1215
STD__
#define STD__
Definition: basics.h:338
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:1689
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:1176
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:1041
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:1194
F_TMatrix::endvec
T * endvec
Definition: f_matrix.h:78
NAMESPACE_END
#define NAMESPACE_END
Definition: basics.h:323
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
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:1600
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:442
TBCI::fabssqr
double fabssqr(const double a)
Definition: basics.h:1157
T
#define T
Definition: bdmatlib.cc:20
TBCICOMP
#define TBCICOMP(n, o, t, s)
Definition: basics.h:979
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:1528
F_TMFORALL_TS
#define F_TMFORALL_TS(op)
Definition: f_matrix.h:642
F_Matrix::element_type
T element_type
Definition: f_matrix.h:1177
F_STDDEF_TSM
#define F_STDDEF_TSM(op)
Definition: f_matrix.h:980
a
const unsigned TMatrix< T > const Matrix< T > * a
Definition: matrix_kernels.h:25
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:1092
v
const Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > long int int char v
&lt; find minimun of func on grid with resolution res
Definition: LM_fit.h:205
MATH__
#define MATH__
Definition: basics.h:339
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:1311
F_Matrix::F_Matrix
F_Matrix(const F_Matrix< T > &m)
Definition: f_matrix.h:1190
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:1027
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:995
F_TSMatrix::eval
F_TSMatrix< T > & eval(F_TMatrix< T > *=0)
Definition: f_matrix.h:961
F_Matrix::operator()
T & operator()(const unsigned int, const unsigned int) const
Definition: f_matrix.h:1624
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
F_TSMatrix::fac
fac(f)
rowcolvec
rowcolvec
Definition: matrix.h:70
Generated by   doxygen 1.8.5