Imagine++
SymMatrix.h
1// ===========================================================================
2// Imagine++ Libraries
3// Copyright (C) Imagine
4// For detailed information: http://imagine.enpc.fr/software
5// ===========================================================================
6
7namespace 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; }
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 }
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 }
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 }
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 }
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 }
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 }
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 }
428 SymMatrix C=clone();
429 return (C*=x);
430 }
439 friend inline SymMatrix operator *(T x, const SymMatrix& B) {return B*x;}
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>
575 Matrix<T> Matrix<T>::operator *(const SymMatrix<T>& B) const
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}
Array of variable size.
Definition: Array.h:20
size_t size() const
Size.
Definition: Array.h:194
Array & fill(const T &x)
Filling.
Definition: Array.h:276
const T * const_iterator
Const iterator type.
Definition: Array.h:73
Array & operator=(const Array &A)
Assignment.
Definition: Array.h:154
void setSize(size_t size)
Change size.
Definition: Array.h:143
T * data()
Data pointer (read/write).
Definition: Array.h:217
Matrix of variable size.
Definition: Matrix.h:23
int nrow() const
Number of rows.
Definition: Matrix.h:101
Vector< T > operator*(const Vector< T > &v) const
Product with vector.
Definition: Matrix.h:283
int ncol() const
Number of columns.
Definition: Matrix.h:108
Matrix()
Empty constructor.
Definition: Matrix.h:36
RED GREEN BLUE color.
Definition: Color.h:26
Symmetric Matrix of variable size.
Definition: SymMatrix.h:22
SymMatrix & fill(T x)
Filling.
Definition: SymMatrix.h:133
friend Vector< T > linSolve(const SymMatrix &A, const Vector< T > &b)
Linear system.
Definition: SymMatrix.h:482
SymMatrix(int N)
Constructor (known size).
Definition: SymMatrix.h:43
SymMatrix operator/(T x) const
Scalar division.
Definition: SymMatrix.h:447
Base::const_iterator const_iterator
Const iterator type.
Definition: SymMatrix.h:30
SymMatrix(T *t, int N, bool handleDelete=false)
Constructor (pre-allocated).
Definition: SymMatrix.h:62
SymMatrix(const Matrix< T > &A)
Conversion from Matrix.
Definition: SymMatrix.h:69
SymMatrix & operator*=(T x)
Scalar in place multiplication.
Definition: SymMatrix.h:458
friend void EVD(const SymMatrix &A, Matrix< T > &Q, Vector< T > &Lambda)
Eigen values.
Definition: SymMatrix.h:316
static SymMatrix Identity(int n)
Identity.
Definition: SymMatrix.h:153
const T & operator()(int i, int j) const
Read access.
Definition: SymMatrix.h:167
Base::iterator iterator
Iterator type.
Definition: SymMatrix.h:28
SymMatrix()
Empty constructor.
Definition: SymMatrix.h:36
SymMatrix operator-() const
Opposite.
Definition: SymMatrix.h:417
SymMatrix & operator/=(T x)
Scalar in place division.
Definition: SymMatrix.h:469
int ncol() const
Number of columns.
Definition: SymMatrix.h:124
Matrix< T > getSubMat(int i0, int m, int j0, int n) const
Get sub matrix.
Definition: SymMatrix.h:198
void read(std::istream &in)
Binary read.
Definition: SymMatrix.h:518
SymMatrix clone() const
Cloning.
Definition: SymMatrix.h:91
Vector< T > operator*(const Vector< T > &v) const
Product with vector.
Definition: SymMatrix.h:225
SymMatrix & operator=(const SymMatrix &A)
Assignment.
Definition: SymMatrix.h:84
~SymMatrix()
Destructor.
Definition: SymMatrix.h:76
friend SymMatrix inverse(const SymMatrix &A)
Inverse.
Definition: SymMatrix.h:241
void write(std::ostream &out) const
Binary write.
Definition: SymMatrix.h:506
SymMatrix(const SymMatrix &A)
Copy constructor.
Definition: SymMatrix.h:50
friend std::ostream & operator<<(std::ostream &out, const SymMatrix &A)
ASCII write.
Definition: SymMatrix.h:533
friend SymMatrix posDefInverse(const SymMatrix &A)
Inverse (variant).
Definition: SymMatrix.h:267
SymMatrix operator+(const SymMatrix &B) const
Addition.
Definition: SymMatrix.h:372
SymMatrix & operator-=(const SymMatrix &B)
In place Substraction.
Definition: SymMatrix.h:406
friend std::istream & operator>>(std::istream &in, SymMatrix &A)
ASCII read.
Definition: SymMatrix.h:553
Vector< T > getDiagonal() const
Get diagonal.
Definition: SymMatrix.h:212
static SymMatrix Zero(int n)
Zero matrix.
Definition: SymMatrix.h:141
void setSize(int n)
Change sizes.
Definition: SymMatrix.h:107
T & operator()(int i, int j)
Write access.
Definition: SymMatrix.h:182
int nrow() const
Number of rows.
Definition: SymMatrix.h:117
friend T det(const SymMatrix &A)
Determinant.
Definition: SymMatrix.h:293
SymMatrix & operator+=(const SymMatrix &B)
In place Addition.
Definition: SymMatrix.h:383
Vector of variable size.
Definition: Vector.h:24
Vector & fill(T x)
Filling.
Definition: Vector.h:93
Vector clone() const
Cloning.
Definition: Vector.h:80
Imagine++ namespace.