All Classes Namespaces Functions 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 
72  Matrix() : memory(), ele(0), n(0), lda(0) {
73  }
74 
75 // template<class Ini=All<AT> >
76 // Matrix(int n_, Ini ini=All<AT>()) : memory(n_*n_), ele((AT*)memory.get()), n(n_), lda(n_) {
77 // init(ini);
78 // }
79 // template<class Ini=All<AT> >
80 // Matrix(int m_, int n_, Ini ini=All<AT>()) : memory(n_*n_), ele((AT*)memory.get()), n(n_), lda(n_) {
81 // init(ini);
82 // }
83 
84  Matrix(int n_, Noinit) : memory(n_*n_), ele((AT*)memory.get()), n(n_), lda(n_) { }
85  Matrix(int n_, Init ini=INIT, const AT &a=0) : memory(n_*n_), ele((AT*)memory.get()), n(n_), lda(n_) { init(a); }
86  Matrix(int n_, Eye ini, const AT &a=1) : memory(n_*n_), ele((AT*)memory.get()), n(n_), lda(n_) { init(ini,a); }
87  Matrix(int m_, int n_, Noinit) : memory(n_*n_), ele((AT*)memory.get()), n(n_), lda(n_) { }
88  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); }
89  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); }
90 
98  Matrix(const Matrix<Symmetric,Ref,Ref,AT> &A) : memory(A.memory), ele(A.ele) , n(A.n), lda(A.lda) { }
99 
100 
105  explicit Matrix(const Matrix<General,Ref,Ref,AT>& A) : memory(A.memory), ele(A.ele) , n(A.n), lda(A.lda) {
106 #ifndef FMATVEC_NO_SIZE_CHECK
107  assert(A.rows() == A.cols());
108 #endif
109  }
110 
111  template<class Type, class Row, class Col>
112  explicit Matrix(const Matrix<Type,Row,Col,AT> &A) : memory(A.rows()*A.cols()), ele((AT*)memory.get()), n(A.cols()), lda(A.cols()) {
113 #ifndef FMATVEC_NO_SIZE_CHECK
114  assert(A.rows() == A.cols());
115 #endif
116 
117  deepCopy(A);
118  }
119 
127  Matrix(int n_, AT* ele_) : memory(), ele(ele_), n(n_), lda(n_) { }
128 
131  ~Matrix() { }
132 
133  Matrix<Symmetric,Ref,Ref,AT>& resize() {
134  n = lda = 0;
135  memory.resize(0);
136  ele = 0;
137  return *this;
138  }
139 
140  Matrix<Symmetric,Ref,Ref,AT>& resize(int n_, Noinit) {
141  n = n_; lda = n;
142  memory.resize(n*n);
143  ele = (AT*)memory.get();
144  return *this;
145  }
146 
147  Matrix<Symmetric,Ref,Ref,AT>& resize(int n, Init ini=INIT, const AT &a=0) { return resize(n,Noinit()).init(a); }
148 
149  Matrix<Symmetric,Ref,Ref,AT>& resize(int n, Eye ini, const AT &a=1) { return resize(n,Noinit()).init(ini,a); }
150 
159  inline Matrix<Symmetric,Ref,Ref,AT>& operator=(const Matrix<Symmetric,Ref,Ref,AT> &A);
160 
161  template<class Type, class Row, class Col>
162  inline Matrix<Symmetric,Ref,Ref,AT>& operator=(const Matrix<Type,Row,Col,AT> &A);
163 
170  template<class Type, class Row, class Col>
171  inline Matrix<Symmetric,Ref,Ref,AT>& operator<<(const Matrix<Type,Row,Col,AT> &A);
172 
179  inline Matrix<Symmetric,Ref,Ref,AT>& operator>>(const Matrix<Symmetric,Ref,Ref,AT> &A);
180 
192  AT& operator()(int i, int j) {
193 #ifndef FMATVEC_NO_BOUNDS_CHECK
194  assert(i>=0);
195  assert(j>=0);
196  assert(i<n);
197  assert(j<n);
198 #endif
199  return e(i,j);
200  };
201 
206  const AT& operator()(int i, int j) const {
207 #ifndef FMATVEC_NO_BOUNDS_CHECK
208  assert(i>=0);
209  assert(j>=0);
210  assert(i<n);
211  assert(j<n);
212 #endif
213 
214  return e(i,j);
215  };
216 
217  AT& ei(int i, int j) {
218  return ele[i+j*lda];
219  };
220 
221  const AT& ei(int i, int j) const {
222  return ele[i+j*lda];
223  };
224 
225  AT& ej(int i, int j) {
226  return ele[i*lda+j];
227  };
228 
229  const AT& ej(int i, int j) const {
230  return ele[i*lda+j];
231  };
232 
233  AT& e(int i, int j) {
234  return j > i ? ej(i,j) : ei(i,j);
235  };
236 
237  const AT& e(int i, int j) const {
238  return j > i ? ej(i,j) : ei(i,j);
239  };
240 
246  AT* operator()() {return ele;};
247 
252  const AT* operator()() const {return ele;};
253 
258  int size() const {return n;};
259 
264  int rows() const {return n;};
265 
270  int cols() const {return n;};
271 
276  int ldim() const {return lda;};
277 
285  const CBLAS_ORDER blasOrder() const {
286  return CblasColMajor;
287  };
288 
296  const CBLAS_UPLO blasUplo() const {
297  return CblasLower;
298  };
299 
305  inline Matrix<Symmetric,Ref,Ref,AT> copy() const;
306 
314  inline Matrix<Symmetric,Ref,Ref,AT>& init(const AT &a);
315  inline Matrix<Symmetric,Ref,Ref,AT>& init(Init, const AT &a=0) { return init(a); };
316  inline Matrix<Symmetric,Ref,Ref,AT>& init(Eye, const AT &a=1);
317  inline Matrix<Symmetric,Ref,Ref,AT>& init(Noinit, const AT &a=0) { return *this; }
318 
330  inline Matrix<General,Ref,Ref,AT> operator()(const Index &I, const Index &J);
331 
336  inline const Matrix<General,Ref,Ref,AT> operator()(const Index &I, const Index &J) const;
337 
349  inline Matrix<Symmetric,Ref,Ref,AT> operator()(const Index &I);
350 
355  inline const Matrix<Symmetric,Ref,Ref,AT> operator()(const Index &I) const;
356 
368  inline Matrix<General,Ref,Ref,AT> operator()(int i1, int j1, int i2, int j2);
369 
374  inline const Matrix<General,Ref,Ref,AT> operator()(int i1, int j1, int i2, int j2) const;
375 
380  inline operator std::vector<std::vector<AT> >();
381  };
382 
383  template <class AT>
385 
386  n=A.n;
387  memory = A.memory;
388  ele = A.ele;
389  lda = A.lda;
390 
391  return *this;
392  }
393 
394  template <class AT>
396 
397  if(!ele) {
398  n = A.size();
399  lda = n;
400  memory.resize(n*n);
401  ele = (AT*)memory.get();
402  } else {
403 #ifndef FMATVEC_NO_SIZE_CHECK
404  assert(n == A.size());
405 #endif
406  }
407 
408  deepCopy(A);
409 
410  return *this;
411  }
412 
413  template <class AT> template< class Type, class Row, class Col>
415 
416  if(!ele) {
417  n = A.rows();
418  lda = n;
419  memory.resize(n*n);
420  ele = (AT*)memory.get();
421  } else {
422 #ifndef FMATVEC_NO_SIZE_CHECK
423  assert(A.rows() == A.cols());
424  assert(n == A.rows());
425 #endif
426  }
427 
428  deepCopy(A);
429 
430  return *this;
431  }
432 
433 
434  template <class AT> template< class Type, class Row, class Col>
436 
437 #ifndef FMATVEC_NO_SIZE_CHECK
438  assert(A.rows() == A.cols());
439 #endif
440 
441  if(n!=A.rows()) {
442  n = A.rows();
443  lda = n;
444  memory.resize(n*n);
445  ele = (AT*)memory.get();
446  }
447 
448  deepCopy(A);
449 
450  return *this;
451  }
452 
453  template <class AT>
455 
456  for(int i=0; i<rows(); i++)
457  for(int j=i; j<cols(); j++)
458  ej(i,j) = val;
459 
460  return *this;
461  }
462 
463  template <class AT>
465 
466  for(int i=0; i<size(); i++) {
467  ej(i,i) = val;
468  for(int j=i+1; j<size(); j++) {
469  ej(i,j) = 0;
470  }
471  }
472  return *this;
473  }
474 
475  template <class AT>
477 
478  Matrix<Symmetric,Ref,Ref,AT> A(n,NONINIT);
479  A.deepCopy(*this);
480 
481  return A;
482  }
483 
484  template <class AT>
485  inline const Matrix<General,Ref,Ref,AT> Matrix<Symmetric,Ref,Ref,AT>::operator()(int i1, int j1, int i2, int j2) const {
486  return operator()(Index(i1,i2),Index(j1,j2));
487  }
488 
489  template <class AT>
491  return operator()(Index(i1,i2),Index(j1,j2));
492  }
493 
494  template <class AT>
496 #ifndef FMATVEC_NO_BOUNDS_CHECK
497  assert(I.end()<n);
498  assert(J.end()<n);
499 #endif
500 
501  if(I.start() >= J.start()) {
502  assert(J.end() <= I.start());
503  return Matrix<General,Ref,Ref,AT>(I.end()-I.start()+1,J.end()-J.start()+1,lda,false,memory,elePtr(I.start(),J.start()));
504  } else {
505  assert(I.end() <= J.start());
506  return Matrix<General,Ref,Ref,AT>(I.end()-I.start()+1,J.end()-J.start()+1,lda,true,memory,elePtr(J.start(),I.start()));
507  }
508  }
509 
510  template <class AT>
512 #ifndef FMATVEC_NO_BOUNDS_CHECK
513  assert(I.end()<n);
514  assert(J.end()<n);
515 #endif
516 
517  if(I.start() >= J.start()) {
518  assert(J.end() <= I.start());
519  return Matrix<General,Ref,Ref,AT>(I.end()-I.start()+1,J.end()-J.start()+1,lda,false,memory,elePtr(I.start(),J.start()));
520  } else {
521  assert(I.end() <= J.start());
522  return Matrix<General,Ref,Ref,AT>(I.end()-I.start()+1,J.end()-J.start()+1,lda,true,memory,elePtr(J.start(),I.start()));
523  }
524  }
525 
526  template <class AT>
528 #ifndef FMATVEC_NO_BOUNDS_CHECK
529  assert(I.end()<n);
530 #endif
531 
532  return Matrix<Symmetric,Ref,Ref,AT>(I.end()-I.start()+1,lda,memory,elePtr(I.start(),I.start()));
533  }
534 
535  template <class AT>
537 #ifndef FMATVEC_NO_BOUNDS_CHECK
538  assert(I.end()<n);
539 #endif
540 
541  return Matrix<Symmetric,Ref,Ref,AT>(I.end()-I.start()+1,lda,memory,elePtr(I.start(),I.start()));
542  }
543 
544  template <class AT>
545  inline Matrix<Symmetric,Ref,Ref,AT>::operator std::vector<std::vector<AT> >() {
546  std::vector<std::vector<AT> > ret(rows());
547  for(int r=0; r<rows(); r++) {
548  ret[r].resize(cols());
549  for(int c=0; c<cols(); c++)
550  ret[r][c]=operator()(r,c);
551  }
552  return ret;
553  }
554 
556 
557  template <class AT>
559  for(int i=0; i<n; i++)
560  for(int j=i; j<n; j++)
561  ej(i,j) = A.ej(i,j);
562  }
563 
564  template <class AT> template <class Type, class Row, class Col>
565  inline void Matrix<Symmetric,Ref,Ref,AT>::deepCopy(const Matrix<Type,Row,Col,AT> &A) {
566  for(int i=0; i<n; i++)
567  for(int j=i; j<n; j++)
568  ej(i,j) = A.e(i,j);
569  }
570 
572 
573 }
574 
575 #endif
Matrix(const Matrix< General, Ref, Ref, AT > &A)
Element operator.
Definition: symmetric_matrix.h:105
AT & operator()(int i, int j)
Element operator.
Definition: symmetric_matrix.h:192
This is the basic matrix class for arbitrary matrices.
Definition: matrix.h:56
Shape class for symmetric matrices.
Definition: types.h:116
AT & operator()(int i, int j)
Standard constructor.
Definition: matrix.h:85
int rows() const
Number of rows.
Definition: general_matrix.h:270
Definition: matrix.h:39
int rows() const
Number of rows.
Definition: symmetric_matrix.h:264
int rows() const
Number of rows.
int end() const
Last element.
Definition: index.h:85
int cols() const
Number of columns.
Definition: types.h:62
const CBLAS_UPLO blasUplo() const
Symmetry convention.
Definition: symmetric_matrix.h:296
const CBLAS_ORDER blasOrder() const
Storage convention.
Definition: symmetric_matrix.h:285
Matrix< Symmetric, Ref, Ref, AT > & init(const AT &a)
Initialization.
Definition: symmetric_matrix.h:454
int cols() const
Number of columns.
Definition: symmetric_matrix.h:270
Matrix(int n_, AT *ele_)
Regular Constructor.
Definition: symmetric_matrix.h:127
This is a matrix class for general matrices.
Definition: general_matrix.h:40
int size() const
Size.
Definition: symmetric_matrix.h:258
AT * operator()()
Pointer operator.
Definition: symmetric_matrix.h:246
std::istream & operator>>(std::istream &is, Matrix< Type, Row, Col, AT > &A)
Matrix input.
Definition: matrix.h:170
const AT & operator()(int i, int j) const
Element operator.
Definition: symmetric_matrix.h:206
This is an index class for creating submatrices.
Definition: index.h:34
int ldim() const
Leading dimension.
Definition: symmetric_matrix.h:276
Matrix(const Matrix< Symmetric, Ref, Ref, AT > &A)
Copy Constructor.
Definition: symmetric_matrix.h:98
Definition: matrix.h:38
This is a matrix class for symmetric matrices.
Definition: symmetric_matrix.h:39
int start() const
First element.
Definition: index.h:79
const AT * operator()() const
Pointer operator.
Definition: symmetric_matrix.h:252
Matrix()
Standard constructor.
Definition: symmetric_matrix.h:72
~Matrix()
Destructor.
Definition: symmetric_matrix.h:131
int cols() const
Number of columns.
Definition: general_matrix.h:276
Definition: matrix.h:40

Impressum / Disclaimer / Datenschutz Generated by doxygen 1.8.5 Valid HTML