All Classes Namespaces Functions Variables Typedefs Enumerations Pages
symmetric_matrix.h
1 /* Copyright (C) 2003-2005 Martin Förg
2 
3  * This library is free software; you can redistribute it and/or
4  * modify it under the terms of the GNU Lesser General Public
5  * License as published by the Free Software Foundation; either
6  * version 2.1 of the License, or (at your option) any later version.
7  *
8  * This library is distributed in the hope that it will be useful,
9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11  * Lesser General Public License for more details.
12  *
13  * You should have received a copy of the GNU Lesser General Public
14  * License along with this library; if not, write to the Free Software
15  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
16  *
17  * Contact:
18  * martin.o.foerg@googlemail.com
19  *
20  */
21 
22 #ifndef symmetric_matrix_h
23 #define symmetric_matrix_h
24 
25 #include "types.h"
26 #include "matrix.h"
27 #include "_memory.h"
28 
29 namespace fmatvec {
30 
39  template <class AT> class Matrix<Symmetric,Ref,Ref,AT> {
40 
41  protected:
42 
44 
45  Memory<AT> memory;
46  AT *ele;
47  int n;
48  int lda;
49 
50  template <class Type, class Row, class Col> inline void deepCopy(const Matrix<Type,Row,Col,AT> &A);
51  inline void deepCopy(const Matrix<Symmetric,Ref,Ref,AT> &A);
52 
53  const AT* elePtr(int i, int j) const {
54  return j > i ? ele+i*lda+j : ele+i+j*lda;
55  };
56 
57  AT* elePtr(int i, int j) {
58  return j > i ? ele+i*lda+j : ele+i+j*lda;
59  };
60 
61  Matrix(int n_, int lda_, Memory<AT> memory_, const AT* ele_) : memory(memory_), ele((AT*)ele_), n(n_), lda(lda_) {
62  }
63 
65 
66  public:
67 
68  typedef AT AtomicType;
69 
74  Matrix() : memory(), ele(0), n(0), lda(0) {
75  }
76 
77 // template<class Ini=All<AT> >
78 // Matrix(int n_, Ini ini=All<AT>()) : memory(n_*n_), ele((AT*)memory.get()), n(n_), lda(n_) {
79 // init(ini);
80 // }
81 // template<class Ini=All<AT> >
82 // Matrix(int m_, int n_, Ini ini=All<AT>()) : memory(n_*n_), ele((AT*)memory.get()), n(n_), lda(n_) {
83 // init(ini);
84 // }
85 
86  Matrix(int n_, Noinit) : memory(n_*n_), ele((AT*)memory.get()), n(n_), lda(n_) { }
87  Matrix(int n_, Init ini=INIT, const AT &a=0) : memory(n_*n_), ele((AT*)memory.get()), n(n_), lda(n_) { init(a); }
88  Matrix(int n_, Eye ini, const AT &a=1) : memory(n_*n_), ele((AT*)memory.get()), n(n_), lda(n_) { init(ini,a); }
89  Matrix(int m_, int n_, Noinit) : memory(n_*n_), ele((AT*)memory.get()), n(n_), lda(n_) { }
90  Matrix(int m_, int n_, Init ini=INIT, const AT &a=0) : memory(n_*n_), ele((AT*)memory.get()), n(n_), lda(n_) { init(a); }
91  Matrix(int m_, int n_, Eye ini, const AT &a=1) : memory(n_*n_), ele((AT*)memory.get()), n(n_), lda(n_) { init(ini,a); }
92 
100  Matrix(const Matrix<Symmetric,Ref,Ref,AT> &A) : memory(A.memory), ele(A.ele) , n(A.n), lda(A.lda) { }
101 
102 
107  explicit Matrix(const Matrix<General,Ref,Ref,AT>& A) : memory(A.memory), ele(A.ele) , n(A.n), lda(A.lda) {
108 #ifndef FMATVEC_NO_SIZE_CHECK
109  assert(A.rows() == A.cols());
110 #endif
111  }
112 
113  template<class Type, class Row, class Col>
114  explicit Matrix(const Matrix<Type,Row,Col,AT> &A) : memory(A.rows()*A.cols()), ele((AT*)memory.get()), n(A.cols()), lda(A.cols()) {
115 #ifndef FMATVEC_NO_SIZE_CHECK
116  assert(A.rows() == A.cols());
117 #endif
118 
119  deepCopy(A);
120  }
121 
129  Matrix(int n_, AT* ele_) : memory(), ele(ele_), n(n_), lda(n_) { }
130 
133  ~Matrix() { }
134 
135  Matrix<Symmetric,Ref,Ref,AT>& resize() {
136  n = lda = 0;
137  memory.resize(0);
138  ele = 0;
139  return *this;
140  }
141 
142  Matrix<Symmetric,Ref,Ref,AT>& resize(int n_, Noinit) {
143  n = n_; lda = n;
144  memory.resize(n*n);
145  ele = (AT*)memory.get();
146  return *this;
147  }
148 
149  Matrix<Symmetric,Ref,Ref,AT>& resize(int n, Init ini=INIT, const AT &a=0) { return resize(n,Noinit()).init(a); }
150 
151  Matrix<Symmetric,Ref,Ref,AT>& resize(int n, Eye ini, const AT &a=1) { return resize(n,Noinit()).init(ini,a); }
152 
155  void resize(int n, int m) {
156  if(n!=m)
157  throw std::runtime_error("A symmetric matrix cannot have different dimensions for rows and columns.");
158  resize(n);
159  }
160 
162  bool transposed() {
163  return false;
164  }
165 
175 
176  template<class Type, class Row, class Col>
177  inline Matrix<Symmetric,Ref,Ref,AT>& operator=(const Matrix<Type,Row,Col,AT> &A);
178 
185  template<class Type, class Row, class Col>
186  inline Matrix<Symmetric,Ref,Ref,AT>& operator<<(const Matrix<Type,Row,Col,AT> &A);
187 
195 
207  AT& operator()(int i, int j) {
208 #ifndef FMATVEC_NO_BOUNDS_CHECK
209  assert(i>=0);
210  assert(j>=0);
211  assert(i<n);
212  assert(j<n);
213 #endif
214  return e(i,j);
215  };
216 
221  const AT& operator()(int i, int j) const {
222 #ifndef FMATVEC_NO_BOUNDS_CHECK
223  assert(i>=0);
224  assert(j>=0);
225  assert(i<n);
226  assert(j<n);
227 #endif
228 
229  return e(i,j);
230  };
231 
232  AT& ei(int i, int j) {
233  return ele[i+j*lda];
234  };
235 
236  const AT& ei(int i, int j) const {
237  return ele[i+j*lda];
238  };
239 
240  AT& ej(int i, int j) {
241  return ele[i*lda+j];
242  };
243 
244  const AT& ej(int i, int j) const {
245  return ele[i*lda+j];
246  };
247 
248  AT& e(int i, int j) {
249  return j > i ? ej(i,j) : ei(i,j);
250  };
251 
252  const AT& e(int i, int j) const {
253  return j > i ? ej(i,j) : ei(i,j);
254  };
255 
261  AT* operator()() {return ele;};
262 
267  const AT* operator()() const {return ele;};
268 
273  int size() const {return n;};
274 
279  int rows() const {return n;};
280 
285  int cols() const {return n;};
286 
291  int ldim() const {return lda;};
292 
300  const CBLAS_ORDER blasOrder() const {
301  return CblasColMajor;
302  };
303 
311  const CBLAS_UPLO blasUplo() const {
312  return CblasLower;
313  };
314 
320  inline Matrix<Symmetric,Ref,Ref,AT> copy() const;
321 
329  inline Matrix<Symmetric,Ref,Ref,AT>& init(const AT &a);
330  inline Matrix<Symmetric,Ref,Ref,AT>& init(Init, const AT &a=0) { return init(a); };
331  inline Matrix<Symmetric,Ref,Ref,AT>& init(Eye, const AT &a=1);
332  inline Matrix<Symmetric,Ref,Ref,AT>& init(Noinit, const AT &a=0) { return *this; }
333 
345  inline Matrix<General,Ref,Ref,AT> operator()(const Range<Var,Var> &I, const Range<Var,Var> &J);
346 
351  inline const Matrix<General,Ref,Ref,AT> operator()(const Range<Var,Var> &I, const Range<Var,Var> &J) const;
352 
364  inline Matrix<Symmetric,Ref,Ref,AT> operator()(const Range<Var,Var> &I);
365 
370  inline const Matrix<Symmetric,Ref,Ref,AT> operator()(const Range<Var,Var> &I) const;
371 
383  inline Matrix<General,Ref,Ref,AT> operator()(int i1, int j1, int i2, int j2);
384 
389  inline const Matrix<General,Ref,Ref,AT> operator()(int i1, int j1, int i2, int j2) const;
390 
395  inline operator std::vector<std::vector<AT> >();
396  };
397 
398  template <class AT>
400 
401  n=A.n;
402  memory = A.memory;
403  ele = A.ele;
404  lda = A.lda;
405 
406  return *this;
407  }
408 
409  template <class AT>
411 
412  if(!ele) {
413  n = A.size();
414  lda = n;
415  memory.resize(n*n);
416  ele = (AT*)memory.get();
417  } else {
418 #ifndef FMATVEC_NO_SIZE_CHECK
419  assert(n == A.size());
420 #endif
421  }
422 
423  deepCopy(A);
424 
425  return *this;
426  }
427 
428  template <class AT> template< class Type, class Row, class Col>
430 
431  if(!ele) {
432  n = A.rows();
433  lda = n;
434  memory.resize(n*n);
435  ele = (AT*)memory.get();
436  } else {
437 #ifndef FMATVEC_NO_SIZE_CHECK
438  assert(A.rows() == A.cols());
439  assert(n == A.rows());
440 #endif
441  }
442 
443  deepCopy(A);
444 
445  return *this;
446  }
447 
448 
449  template <class AT> template< class Type, class Row, class Col>
451 
452 #ifndef FMATVEC_NO_SIZE_CHECK
453  assert(A.rows() == A.cols());
454 #endif
455 
456  if(n!=A.rows()) {
457  n = A.rows();
458  lda = n;
459  memory.resize(n*n);
460  ele = (AT*)memory.get();
461  }
462 
463  deepCopy(A);
464 
465  return *this;
466  }
467 
468  template <class AT>
470 
471  for(int i=0; i<rows(); i++)
472  for(int j=i; j<cols(); j++)
473  ej(i,j) = val;
474 
475  return *this;
476  }
477 
478  template <class AT>
480 
481  for(int i=0; i<size(); i++) {
482  ej(i,i) = val;
483  for(int j=i+1; j<size(); j++) {
484  ej(i,j) = 0;
485  }
486  }
487  return *this;
488  }
489 
490  template <class AT>
492 
493  Matrix<Symmetric,Ref,Ref,AT> A(n,NONINIT);
494  A.deepCopy(*this);
495 
496  return A;
497  }
498 
499  template <class AT>
500  inline const Matrix<General,Ref,Ref,AT> Matrix<Symmetric,Ref,Ref,AT>::operator()(int i1, int j1, int i2, int j2) const {
501  return operator()(Range<Var,Var>(i1,i2),Range<Var,Var>(j1,j2));
502  }
503 
504  template <class AT>
506  return operator()(Range<Var,Var>(i1,i2),Range<Var,Var>(j1,j2));
507  }
508 
509  template <class AT>
511 #ifndef FMATVEC_NO_BOUNDS_CHECK
512  assert(I.end()<n);
513  assert(J.end()<n);
514 #endif
515 
516  if(I.start() >= J.start()) {
517  assert(J.end() <= I.start());
518  return Matrix<General,Ref,Ref,AT>(I.end()-I.start()+1,J.end()-J.start()+1,lda,false,memory,elePtr(I.start(),J.start()));
519  } else {
520  assert(I.end() <= J.start());
521  return Matrix<General,Ref,Ref,AT>(I.end()-I.start()+1,J.end()-J.start()+1,lda,true,memory,elePtr(J.start(),I.start()));
522  }
523  }
524 
525  template <class AT>
527 #ifndef FMATVEC_NO_BOUNDS_CHECK
528  assert(I.end()<n);
529  assert(J.end()<n);
530 #endif
531 
532  if(I.start() >= J.start()) {
533  assert(J.end() <= I.start());
534  return Matrix<General,Ref,Ref,AT>(I.end()-I.start()+1,J.end()-J.start()+1,lda,false,memory,elePtr(I.start(),J.start()));
535  } else {
536  assert(I.end() <= J.start());
537  return Matrix<General,Ref,Ref,AT>(I.end()-I.start()+1,J.end()-J.start()+1,lda,true,memory,elePtr(J.start(),I.start()));
538  }
539  }
540 
541  template <class AT>
543 #ifndef FMATVEC_NO_BOUNDS_CHECK
544  assert(I.end()<n);
545 #endif
546 
547  return Matrix<Symmetric,Ref,Ref,AT>(I.end()-I.start()+1,lda,memory,elePtr(I.start(),I.start()));
548  }
549 
550  template <class AT>
552 #ifndef FMATVEC_NO_BOUNDS_CHECK
553  assert(I.end()<n);
554 #endif
555 
556  return Matrix<Symmetric,Ref,Ref,AT>(I.end()-I.start()+1,lda,memory,elePtr(I.start(),I.start()));
557  }
558 
559  template <class AT>
560  inline Matrix<Symmetric,Ref,Ref,AT>::operator std::vector<std::vector<AT> >() {
561  std::vector<std::vector<AT> > ret(rows());
562  for(int r=0; r<rows(); r++) {
563  ret[r].resize(cols());
564  for(int c=0; c<cols(); c++)
565  ret[r][c]=operator()(r,c);
566  }
567  return ret;
568  }
569 
571 
572  template <class AT>
574  for(int i=0; i<n; i++)
575  for(int j=i; j<n; j++)
576  ej(i,j) = A.ej(i,j);
577  }
578 
579  template <class AT> template <class Type, class Row, class Col>
580  inline void Matrix<Symmetric,Ref,Ref,AT>::deepCopy(const Matrix<Type,Row,Col,AT> &A) {
581  for(int i=0; i<n; i++)
582  for(int j=i; j<n; j++)
583  ej(i,j) = A.e(i,j);
584  }
585 
587 
588 }
589 
590 #endif
Matrix(const Matrix< General, Ref, Ref, AT > &A)
Element operator.
Definition: symmetric_matrix.h:107
AT & operator()(int i, int j)
Element operator.
Definition: symmetric_matrix.h:207
This is the basic matrix class for arbitrary matrices.
Definition: fmatvec.h:41
Shape class for symmetric matrices.
Definition: types.h:116
AT & operator()(int i, int j)
Standard constructor.
Definition: matrix.h:86
int rows() const
Number of rows.
Definition: general_matrix.h:272
Definition: matrix.h:39
bool transposed()
The storage format of a symmetric matrix is fortran-storage order -&gt; return always false...
Definition: symmetric_matrix.h:162
int rows() const
Number of rows.
Definition: symmetric_matrix.h:279
int rows() const
Number of rows.
int cols() const
Number of columns.
Definition: types.h:62
const CBLAS_UPLO blasUplo() const
Symmetry convention.
Definition: symmetric_matrix.h:311
const CBLAS_ORDER blasOrder() const
Storage convention.
Definition: symmetric_matrix.h:300
Matrix< Symmetric, Ref, Ref, AT > & init(const AT &a)
Initialization.
Definition: symmetric_matrix.h:469
int cols() const
Number of columns.
Definition: symmetric_matrix.h:285
Matrix(int n_, AT *ele_)
Regular Constructor.
Definition: symmetric_matrix.h:129
This is a matrix class for general matrices.
Definition: general_matrix.h:40
int size() const
Size.
Definition: symmetric_matrix.h:273
AT * operator()()
Pointer operator.
Definition: symmetric_matrix.h:261
std::istream & operator>>(std::istream &is, Matrix< Type, Row, Col, AT > &A)
Matrix input.
Definition: matrix.h:171
const AT & operator()(int i, int j) const
Element operator.
Definition: symmetric_matrix.h:221
int ldim() const
Leading dimension.
Definition: symmetric_matrix.h:291
Matrix(const Matrix< Symmetric, Ref, Ref, AT > &A)
Copy Constructor.
Definition: symmetric_matrix.h:100
Definition: matrix.h:38
This is a matrix class for symmetric matrices.
Definition: symmetric_matrix.h:39
int end() const
Last element.
Definition: range.h:95
const AT * operator()() const
Pointer operator.
Definition: symmetric_matrix.h:267
Matrix()
Standard constructor.
Definition: symmetric_matrix.h:74
void resize(int n, int m)
Definition: symmetric_matrix.h:155
This is an index class for creating submatrices.
Definition: range.h:44
~Matrix()
Destructor.
Definition: symmetric_matrix.h:133
int start() const
First element.
Definition: range.h:89
int cols() const
Number of columns.
Definition: general_matrix.h:278
Definition: matrix.h:40

Impressum / Disclaimer / Datenschutz Generated by doxygen 1.8.5 Valid HTML