Imagine++
SymMatrix.h
1 // ===========================================================================
2 // Imagine++ Libraries
3 // Copyright (C) Imagine
4 // For detailed information: http://imagine.enpc.fr/software
5 // ===========================================================================
6 
7 namespace Imagine {
10 
11  template <typename T> class Matrix;
12 
21  template <typename T>
22  class SymMatrix : public Array<T> {
23  typedef Array<T> Base;
24  int _n;
25  size_t sz(int n) { return size_t(n)*(n+1)/2; }
26  public:
28  typedef typename Base::iterator iterator;
36  SymMatrix() : Base(), _n(0) {}
43  explicit SymMatrix(int N) : Base(sz(N)), _n(N) {}
50  SymMatrix(const SymMatrix& A) : Base(A), _n(A._n) {}
62  SymMatrix(T* t, int N,bool handleDelete=false) : Base(t,sz(N),handleDelete), _n(N) { }
69  SymMatrix(const Matrix<T>& A) : Base(sz(A.nrow())), _n(A.nrow()) { // Upper part
70  assert(A.nrow()==A.ncol());
71  for (int j=0; j<A.nrow(); j++)
72  equalize(j+1,A.data()+size_t(j)*A.nrow(),1,this->data()+sz(j),1);
73  }
84  SymMatrix& operator=(const SymMatrix& A) { Base::operator =(A); _n = A._n; return *this; }
91  SymMatrix clone() const
92  {
93  SymMatrix A(nrow());
94  equalize((int)this->size(),this->data(),1,A.data(),1);
95  return A;
96  }
107  void setSize(int n) {
108  Base::setSize(sz(n));
109  _n = n;
110  }
117  int nrow() const { return _n; }
124  int ncol() const { return _n; }
132 
133  SymMatrix& fill(T x) { Base::fill(x); return *this; }
141  static SymMatrix Zero(int n){
142  SymMatrix Z(n);
143  memset(Z.data(),0,Z.size()*sizeof(T));
144  return Z;
145  }
153  static SymMatrix Identity(int n){
155  for (int i=0;i<n;i++)
156  I(i,i)=T(1);
157  return I;
158  }
167  const T& operator()(int i,int j) const {
168  assert(i>=0 && i<_n && j>=0 && j<_n);
169  if(i<=j)
170  return (*this)[i+size_t(j)*(j+1)/2];
171  else
172  return (*this)[j+size_t(i)*(i+1)/2];
173  }
182  T& operator()(int i,int j) {
183  assert(i>=0 && i<_n && j>=0 && j<_n);
184  if(i<=j)
185  return (*this)[i+size_t(j)*(j+1)/2];
186  else
187  return (*this)[j+size_t(i)*(i+1)/2];
188  }
198  Matrix<T> getSubMat(int i0, int m, int j0, int n) const {
199  assert (i0>=0 && i0+m<=nrow() && j0>=0 && j0+n<=ncol());
200  Matrix<T> retMat(m,n);
201  for(int j=0;j<n;j++)
202  for(int i=0;i<m;i++)
203  retMat(i,j)=(*this)(i0+i,j0+j);
204  return retMat;
205  }
213  Vector<T> v(_n);
214  for (int i=0;i<_n;i++)
215  v[i]=(*this)(i,i);
216  return v;
217  }
225  Vector<T> operator *(const Vector<T> &v) const {
226  assert((size_t)_n==v.size());
227  Vector<T> y(_n);
228  y.fill(T(0));
229  int _nn=_n;
230  vectorProduct(_nn,T(1),this->data(),v.data(),T(0),y.data());
231  return y;
232  }
241  friend SymMatrix inverse(const SymMatrix& A) {
242 
243  SymMatrix invA=A.clone();
244  T *pivots=new T[A.nrow()*A.nrow()];
245  int info,n=A.nrow();
246  // LU Inverse
247 
248  symInverse(&n,invA.data(),pivots,&info);
249  Matrix<T> inv(pivots,A.nrow(),A.nrow());
250  SymMatrix<T> invA2(inv);
251  delete[] pivots;
252  if(info!=0) {
253  std::cerr << "Cannot invert symmetric matrix" << std::endl;
254  invA.fill(T(0));
255  return invA;
256  }
257  return invA2;
258  }
267  friend SymMatrix posDefInverse(const SymMatrix& A) {
268 
269  SymMatrix invA=A.clone();
270  T *pivots=new T[A.nrow()*A.nrow()];
271  int info,n=A.nrow();
272  // LU Inverse
273 
274  symDPInverse(&n,invA.data(),pivots,&info);
275  Matrix<T> inv(pivots,A.nrow(),A.nrow());
276  SymMatrix<T> invA2(inv);
277  delete[] pivots;
278  if(info!=0) {
279  std::cerr << "Cannot invert symmetric matrix" << std::endl;
280  invA.fill(T(0));
281  return invA;
282  }
283  return invA2;
284  }
293  friend T det(const SymMatrix& A)
294  {
295  int info,n=A.nrow();
296  SymMatrix symtemp = A.clone();
297  T d =1;
298  d = symDeterminant(&n,symtemp.data(),&info);
299  if(info!=0) {
300  std::cerr << "Cannot compute determinant" << std::endl;
301  return 0;
302  };
303  return d;
304  }
316  friend void EVD(const SymMatrix& A,Matrix<T>& Q, Vector<T>& Lambda) {
317  SymMatrix symtemp = A.clone();
318  Lambda = Vector<T>(A.nrow());
319  Q = Matrix<T>(A.nrow(),A.nrow());
320  int info,n=A.nrow();
321  symEigenValues(&n,symtemp.data(),Lambda.data(),Q.data(),&info);
322  if(info) {
323  std::cerr << "Cannot compute EVD decomposition" << std::endl;
324  return;
325  }
326  }
334  Matrix<T> operator *(const SymMatrix& B) const {
335  // No blas function
336  assert(ncol()==B.nrow());
337  Matrix<T> C(nrow(),B.ncol());
338  for (int j=0;j<B.ncol();j++)
339  for (int i=0;i<nrow();i++) {
340  C(i,j)=0;
341  for (int k=0;k<ncol();k++)
342  C(i,j)+=(*this)(i,k)*B(k,j);
343  }
344  return C;
345  }
353  Matrix<T> operator *(const Matrix<T>& B) const {
354  // No blas function
355  assert(ncol()==B.nrow());
356  Matrix<T> C(nrow(),B.ncol());
357  for (int j=0;j<B.ncol();j++)
358  for (int i=0;i<nrow();i++) {
359  C(i,j)=0;
360  for (int k=0;k<ncol();k++)
361  C(i,j)+=(*this)(i,k)*B(k,j);
362  }
363  return C;
364  }
372  SymMatrix operator +(const SymMatrix& B) const {
373  SymMatrix C=clone();
374  return (C+=B);
375  }
384  assert(nrow()==B.nrow());
385  combine((int)this->size(), T(1), B.data(), this->data());
386  return *this;
387  }
395  SymMatrix operator -(const SymMatrix& B) const {
396  SymMatrix C=clone();
397  return (C-=B);
398  }
407  assert(nrow()==B.nrow());
408  combine((int)this->size(), T(-1), B.data(), this->data());
409  return *this;
410  }
418  return (*this)*T(-1);
419  }
427  SymMatrix operator *(T x) const {
428  SymMatrix C=clone();
429  return (C*=x);
430  }
439  friend inline SymMatrix operator *(T x, const SymMatrix& B) {return B*x;}
447  SymMatrix operator /(T x) const {
448  SymMatrix C=clone();
449  return (C/=x);
450  }
459  multiply((int)this->size(), x, this->data());
460  return *this;
461  }
470  multiply((int)this->size(), T(1)/x, this->data());
471  return *this;
472  }
482  friend Vector<T> linSolve(const SymMatrix& A,const Vector<T>& b)
483  {
484  assert(b.size()==(size_t)A.nrow());
485  SymMatrix invA=A.clone();
486  Vector<T> x=b.clone();
487  // Bunch Kaufmann factorization
488  int *pivots=new int[A.nrow()];
489  int info,n=A.nrow();
490  // Solve linear system using factorization
491  symSolve(&n,invA.data(),pivots,x.data(),&info);
492  if (info) {
493  x.fill(T(0));
494  std::cerr<<"Cannot solve linear system"<<std::endl;
495  }
496  delete[] pivots;
497  return x;
498  }
506  void write(std::ostream& out) const {
507  int n=nrow();
508  out.write((const char*)&n,sizeof(int));
509  out.write((const char*)this->data(),(std::streamsize)this->size()*sizeof(T));
510  }
518  void read(std::istream& in) {
519  int n;
520  in.read((char*)&n,sizeof(int));
521  setSize(n);
522  in.read((char*)this->data(),(std::streamsize)this->size()*sizeof(T));
523  }
533  friend inline std::ostream& operator<<(std::ostream& out,const SymMatrix& A) {
534  out<<A.nrow()<<" "<<A.ncol()<< std::endl;
535  for (int i=0;i<A.nrow();i++) {
536  for (int j=0;j<A.ncol();j++) {
537  out<<A(i,j);
538  if (j<A.ncol()-1) out<<" ";
539  }
540  out<<std::endl;
541  }
542  return out;
543  }
553  friend inline std::istream& operator>>(std::istream& in,SymMatrix& A) {
554  int m,n;
555  in>>m>>n;
556  assert(m==n);
557  A.setSize(n);
558  for (int i=0;i<A.nrow();i++)
559  for (int j=0;j<A.ncol();j++)
560  in>>A(i,j);
561  return in;
562  }
563  };
564 
565 #ifndef DOXYGEN_SHOULD_SKIP_THIS
566  // Declaration in Matrix.h
567  template <typename T>
568  Matrix<T>::Matrix(const SymMatrix<T>& A) : MultiArray<T,2>(A.nrow(),A.ncol()) {
569  for (int j=0; j<ncol(); j++)
570  for (int i=0; i<nrow(); i++)
571  (*this)(i,j)=A(i,j);
572  }
573  // Declaration in Matrix.h
574  template <typename T>
576  {
577  assert(ncol()==B.nrow());
578  Matrix<T> C(nrow(),B.ncol());
579  for (int j=0;j<B.ncol();j++)
580  for (int i=0;i<nrow();i++) {
581  C(i,j)=0;
582  for (int k=0;k<ncol();k++)
583  C(i,j)+=(*this)(i,k)*B(k,j);
584  }
585  return C;
586  }
587 #endif
588 
590 }
friend std::ostream & operator<<(std::ostream &out, const SymMatrix &A)
ASCII write.
Definition: SymMatrix.h:533
Vector & fill(T x)
Filling.
Definition: Vector.h:93
Vector clone() const
Cloning.
Definition: Vector.h:80
Vector of variable size.
Definition: Vector.h:24
SymMatrix(const SymMatrix &A)
Copy constructor.
Definition: SymMatrix.h:50
Vector< T > getDiagonal() const
Get diagonal.
Definition: SymMatrix.h:212
SymMatrix & fill(T x)
Filling.
Definition: SymMatrix.h:133
Matrix< T > getSubMat(int i0, int m, int j0, int n) const
Get sub matrix.
Definition: SymMatrix.h:198
Vector< T > operator*(const Vector< T > &v) const
Product with vector.
Definition: Matrix.h:283
SymMatrix & operator=(const SymMatrix &A)
Assignment.
Definition: SymMatrix.h:84
SymMatrix()
Empty constructor.
Definition: SymMatrix.h:36
static SymMatrix Identity(int n)
Identity.
Definition: SymMatrix.h:153
static SymMatrix Zero(int n)
Zero matrix.
Definition: SymMatrix.h:141
void setSize(int n)
Change sizes.
Definition: SymMatrix.h:107
SymMatrix & operator+=(const SymMatrix &B)
In place Addition.
Definition: SymMatrix.h:383
Array & operator=(const Array &A)
Assignment.
Definition: Array.h:154
~SymMatrix()
Destructor.
Definition: SymMatrix.h:76
size_t size() const
Size.
Definition: Array.h:194
T & operator()(int i, int j)
Write access.
Definition: SymMatrix.h:182
friend SymMatrix posDefInverse(const SymMatrix &A)
Inverse (variant).
Definition: SymMatrix.h:267
int ncol() const
Number of columns.
Definition: Matrix.h:108
friend void EVD(const SymMatrix &A, Matrix< T > &Q, Vector< T > &Lambda)
Eigen values.
Definition: SymMatrix.h:316
friend std::istream & operator>>(std::istream &in, SymMatrix &A)
ASCII read.
Definition: SymMatrix.h:553
int nrow() const
Number of rows.
Definition: SymMatrix.h:117
SymMatrix & operator*=(T x)
Scalar in place multiplication.
Definition: SymMatrix.h:458
int nrow() const
Number of rows.
Definition: Matrix.h:101
int ncol() const
Number of columns.
Definition: SymMatrix.h:124
Base::const_iterator const_iterator
Const iterator type.
Definition: SymMatrix.h:30
void setSize(size_t size)
Change size.
Definition: Array.h:143
SymMatrix clone() const
Cloning.
Definition: SymMatrix.h:91
Matrix()
Empty constructor.
Definition: Matrix.h:36
T * data()
Data pointer (read/write).
Definition: Array.h:217
Base::iterator iterator
Iterator type.
Definition: SymMatrix.h:28
void read(std::istream &in)
Binary read.
Definition: SymMatrix.h:518
const T & operator()(int i, int j) const
Read access.
Definition: SymMatrix.h:167
SymMatrix(int N)
Constructor (known size).
Definition: SymMatrix.h:43
SymMatrix(const Matrix< T > &A)
Conversion from Matrix.
Definition: SymMatrix.h:69
SymMatrix operator/(T x) const
Scalar division.
Definition: SymMatrix.h:447
SymMatrix & operator-=(const SymMatrix &B)
In place Substraction.
Definition: SymMatrix.h:406
void write(std::ostream &out) const
Binary write.
Definition: SymMatrix.h:506
SymMatrix operator+(const SymMatrix &B) const
Addition.
Definition: SymMatrix.h:372
SymMatrix & operator/=(T x)
Scalar in place division.
Definition: SymMatrix.h:469
SymMatrix(T *t, int N, bool handleDelete=false)
Constructor (pre-allocated).
Definition: SymMatrix.h:62
friend Vector< T > linSolve(const SymMatrix &A, const Vector< T > &b)
Linear system.
Definition: SymMatrix.h:482
Array & fill(const T &x)
Filling.
Definition: Array.h:276
Imagine++ namespace.
Definition: Array.h:7
Matrix of variable size.
Definition: Matrix.h:23
RED GREEN BLUE color.
Definition: Color.h:26
const T * const_iterator
Const iterator type.
Definition: Array.h:73
Vector< T > operator*(const Vector< T > &v) const
Product with vector.
Definition: SymMatrix.h:225
friend T det(const SymMatrix &A)
Determinant.
Definition: SymMatrix.h:293
friend SymMatrix inverse(const SymMatrix &A)
Inverse.
Definition: SymMatrix.h:241
SymMatrix operator-() const
Opposite.
Definition: SymMatrix.h:417
Symmetric Matrix of variable size.
Definition: Matrix.h:12