OpenMEEG
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
symmatrix.h
Go to the documentation of this file.
1 /*
2 Project Name : OpenMEEG
3 
4 © INRIA and ENPC (contributors: Geoffray ADDE, Maureen CLERC, Alexandre
5 GRAMFORT, Renaud KERIVEN, Jan KYBIC, Perrine LANDREAU, Théodore PAPADOPOULO,
6 Emmanuel OLIVI
7 Maureen.Clerc.AT.inria.fr, keriven.AT.certis.enpc.fr,
8 kybic.AT.fel.cvut.cz, papadop.AT.inria.fr)
9 
10 The OpenMEEG software is a C++ package for solving the forward/inverse
11 problems of electroencephalography and magnetoencephalography.
12 
13 This software is governed by the CeCILL-B license under French law and
14 abiding by the rules of distribution of free software. You can use,
15 modify and/ or redistribute the software under the terms of the CeCILL-B
16 license as circulated by CEA, CNRS and INRIA at the following URL
17 "http://www.cecill.info".
18 
19 As a counterpart to the access to the source code and rights to copy,
20 modify and redistribute granted by the license, users are provided only
21 with a limited warranty and the software's authors, the holders of the
22 economic rights, and the successive licensors have only limited
23 liability.
24 
25 In this respect, the user's attention is drawn to the risks associated
26 with loading, using, modifying and/or developing or reproducing the
27 software by the user in light of its specific status of free software,
28 that may mean that it is complicated to manipulate, and that also
29 therefore means that it is reserved for developers and experienced
30 professionals having in-depth computer knowledge. Users are therefore
31 encouraged to load and test the software's suitability as regards their
32 requirements in conditions enabling the security of their systems and/or
33 data to be ensured and, more generally, to use and operate it in the
34 same conditions as regards security.
35 
36 The fact that you are presently reading this means that you have had
37 knowledge of the CeCILL-B license and that you accept its terms.
38 */
39 
40 #pragma once
41 
42 #include <iostream>
43 #include <cstdlib>
44 #include <string>
45 
46 #include <vector.h>
47 #include <linop.h>
48 
49 namespace OpenMEEG {
50 
51  class Matrix;
52 
53  class OPENMEEGMATHS_EXPORT SymMatrix : public LinOp {
54 
55  friend class Vector;
56 
57  utils::RCPtr<LinOpValue> value;
58 
59  public:
60 
61  SymMatrix(): LinOp(0,0,SYMMETRIC,2),value() {}
62 
63  SymMatrix(const char* fname): LinOp(0,0,SYMMETRIC,2),value() { this->load(fname); }
64  SymMatrix(size_t N): LinOp(N,N,SYMMETRIC,2),value(new LinOpValue(size())) { }
65  SymMatrix(size_t M,size_t N): LinOp(N,N,SYMMETRIC,2),value(new LinOpValue(size())) { om_assert(N==M); }
66  SymMatrix(const SymMatrix& S,const DeepCopy): LinOp(S.nlin(),S.nlin(),SYMMETRIC,2),value(new LinOpValue(S.size(),S.data())) { }
67 
68  explicit SymMatrix(const Vector& v);
69  explicit SymMatrix(const Matrix& A);
70 
71  size_t size() const { return nlin()*(nlin()+1)/2; };
72  void info() const ;
73 
74  size_t ncol() const { return nlin(); } // SymMatrix only need num_lines
75  size_t& ncol() { return nlin(); }
76 
77  void alloc_data() { value = new LinOpValue(size()); }
78  void reference_data(const double* array) { value = new LinOpValue(size(),array); }
79 
80  bool empty() const { return value->empty(); }
81  void set(double x) ;
82  double* data() const { return value->data; }
83 
84  inline double operator()(size_t i,size_t j) const;
85  inline double& operator()(size_t i,size_t j) ;
86 
87  Matrix operator()(size_t i_start, size_t i_end, size_t j_start, size_t j_end) const;
88  Matrix submat(size_t istart, size_t isize, size_t jstart, size_t jsize) const;
89  SymMatrix submat(size_t istart, size_t iend) const;
90  Vector getlin(size_t i) const;
91  void setlin(size_t i, const Vector& v);
92  Vector solveLin(const Vector &B) const;
93  void solveLin(Vector * B, int nbvect);
94  Matrix solveLin(Matrix& B) const;
95 
96  const SymMatrix& operator=(const double d);
97 
98  SymMatrix operator+(const SymMatrix& B) const;
99  SymMatrix operator-(const SymMatrix& B) const;
100  SymMatrix operator*(const SymMatrix& B) const;
101  Matrix operator*(const Matrix& B) const;
102  Vector operator*(const Vector& v) const;
103  SymMatrix operator*(double x) const;
104  SymMatrix operator/(double x) const {return (*this)*(1/x);}
105  void operator +=(const SymMatrix& B);
106  void operator -=(const SymMatrix& B);
107  void operator *=(double x);
108  void operator /=(double x) { (*this)*=(1/x); }
109 
110  SymMatrix inverse() const;
111  void invert();
112  SymMatrix posdefinverse() const;
113  double det();
114  // void eigen(Matrix & Z, Vector & D );
115 
116  void save(const char *filename) const;
117  void load(const char *filename);
118 
119  void save(const std::string& s) const { save(s.c_str()); }
120  void load(const std::string& s) { load(s.c_str()); }
121 
122  friend class Matrix;
123  };
124 
125  inline double SymMatrix::operator()(size_t i,size_t j) const {
126  om_assert(i<nlin() && j<nlin());
127  if(i<=j)
128  return data()[i+j*(j+1)/2];
129  else
130  return data()[j+i*(i+1)/2];
131  }
132 
133  inline double& SymMatrix::operator()(size_t i,size_t j) {
134  om_assert(i<nlin() && j<nlin());
135  if(i<=j)
136  return data()[i+j*(j+1)/2];
137  else
138  return data()[j+i*(i+1)/2];
139  }
140 
141  //returns the solution of (this)*X = B
142  inline Vector SymMatrix::solveLin(const Vector &B) const {
143  SymMatrix invA(*this,DEEP_COPY);
144  Vector X(B,DEEP_COPY);
145 
146  #ifdef HAVE_LAPACK
147  // Bunch Kaufman Factorization
148  BLAS_INT *pivots=new BLAS_INT[nlin()];
149  int Info = 0;
150  DSPTRF('U',sizet_to_int(invA.nlin()),invA.data(),pivots,Info);
151  // Inverse
152  DSPTRS('U',sizet_to_int(invA.nlin()),1,invA.data(),pivots,X.data(),sizet_to_int(invA.nlin()),Info);
153 
154  om_assert(Info==0);
155  delete[] pivots;
156  #else
157  std::cout << "solveLin not defined" << std::endl;
158  #endif
159  return X;
160  }
161 
162  // stores in B the solution of (this)*X = B, where B is a set of nbvect vector
163  inline void SymMatrix::solveLin(Vector * B, int nbvect) {
164  SymMatrix invA(*this,DEEP_COPY);
165 
166  #ifdef HAVE_LAPACK
167  // Bunch Kaufman Factorization
168  BLAS_INT *pivots=new BLAS_INT[nlin()];
169  int Info = 0;
170  //char *uplo="U";
171  DSPTRF('U',sizet_to_int(invA.nlin()),invA.data(),pivots,Info);
172  // Inverse
173  for(int i = 0; i < nbvect; i++)
174  DSPTRS('U',sizet_to_int(invA.nlin()),1,invA.data(),pivots,B[i].data(),sizet_to_int(invA.nlin()),Info);
175 
176  om_assert(Info==0);
177  delete[] pivots;
178  #else
179  std::cout << "solveLin not defined" << std::endl;
180  #endif
181  }
182 
183  inline void SymMatrix::operator -=(const SymMatrix &B) {
184  om_assert(nlin()==B.nlin());
185  #ifdef HAVE_BLAS
186  BLAS(daxpy,DAXPY)(sizet_to_int(nlin()*(nlin()+1)/2), -1.0, B.data(), 1, data() , 1);
187  #else
188  for (size_t i=0;i<nlin()*(nlin()+1)/2;i++)
189  data()[i]+=B.data()[i];
190  #endif
191  }
192 
193  inline void SymMatrix::operator +=(const SymMatrix &B) {
194  om_assert(nlin()==B.nlin());
195  #ifdef HAVE_BLAS
196  BLAS(daxpy,DAXPY)(sizet_to_int(nlin()*(nlin()+1)/2), 1.0, B.data(), 1, data() , 1);
197  #else
198  for (size_t i=0;i<nlin()*(nlin()+1)/2;i++)
199  data()[i]+=B.data()[i];
200  #endif
201  }
202 
204  // supposes (*this) is definite positive
205  SymMatrix invA(*this,DEEP_COPY);
206  #ifdef HAVE_LAPACK
207  // U'U factorization then inverse
208  int Info = 0;
209  DPPTRF('U', sizet_to_int(nlin()),invA.data(),Info);
210  DPPTRI('U', sizet_to_int(nlin()),invA.data(),Info);
211  om_assert(Info==0);
212  #else
213  std::cerr << "Positive definite inverse not defined" << std::endl;
214  #endif
215  return invA;
216  }
217 
218  inline double SymMatrix::det() {
219  SymMatrix invA(*this,DEEP_COPY);
220  double d = 1.0;
221  #ifdef HAVE_LAPACK
222  // Bunch Kaufmqn
223  BLAS_INT *pivots=new BLAS_INT[nlin()];
224  int Info = 0;
225  // TUDUtTt
226  DSPTRF('U', sizet_to_int(invA.nlin()), invA.data(), pivots,Info);
227  if (Info<0)
228  std::cout << "Big problem in det (DSPTRF)" << std::endl;
229  for (size_t i = 0; i< nlin(); i++){
230  if (pivots[i] >= 0) {
231  d *= invA(i,i);
232  } else { // pivots[i] < 0
233  if (i < nlin()-1 && pivots[i] == pivots[i+1]) {
234  d *= (invA(i,i)*invA(i+1,i+1)-invA(i,i+1)*invA(i+1,i));
235  i++;
236  } else {
237  std::cout << "Big problem in det" << std::endl;
238  }
239  }
240  }
241  delete[] pivots;
242  #else
243  std::cerr << "Determinant not defined without LAPACK" << std::endl;
244  exit(1);
245  #endif
246  return(d);
247  }
248 
249  // inline void SymMatrix::eigen(Matrix & Z, Vector & D ){
250  // // performs the complete eigen-decomposition.
251  // // (*this) = Z.D.Z'
252  // // -> eigenvector are columns of the Matrix Z.
253  // // (*this).Z[:,i] = D[i].Z[:,i]
254  // #ifdef HAVE_LAPACK
255  // SymMatrix symtemp(*this,DEEP_COPY);
256  // D = Vector(nlin());
257  // Z = Matrix(nlin(),nlin());
258  //
259  // int info;
260  // double lworkd;
261  // int lwork;
262  // int liwork;
263  //
264  // DSPEVD('V','U',sizet_to_int(nlin()),symtemp.data(),D.data(),Z.data(),sizet_to_int(nlin()),&lworkd,-1,&liwork,-1,info);
265  // lwork = (int) lworkd;
266  // double * work = new double[lwork];
267  // BLAS_INT *iwork = new BLAS_INT[liwork];
268  // DSPEVD('V','U',sizet_to_int(nlin()),symtemp.data(),D.data(),Z.data(),sizet_to_int(nlin()),work,lwork,iwork,liwork,info);
269  //
270  // delete[] work;
271  // delete[] iwork;
272  // #endif
273  // }
274 
275  inline SymMatrix SymMatrix::operator +(const SymMatrix &B) const {
276  om_assert(nlin()==B.nlin());
277  SymMatrix C(*this,DEEP_COPY);
278  #ifdef HAVE_BLAS
279  BLAS(daxpy,DAXPY)(sizet_to_int(nlin()*(nlin()+1)/2), 1.0, B.data(), 1, C.data() , 1);
280  #else
281  for (size_t i=0;i<nlin()*(nlin()+1)/2;i++)
282  C.data()[i]+=B.data()[i];
283  #endif
284  return C;
285  }
286 
288  {
289  om_assert(nlin()==B.nlin());
290  SymMatrix C(*this,DEEP_COPY);
291  #ifdef HAVE_BLAS
292  BLAS(daxpy,DAXPY)(sizet_to_int(nlin()*(nlin()+1)/2), -1.0, B.data(), 1, C.data() , 1);
293  #else
294  for (size_t i=0;i<nlin()*(nlin()+1)/2;i++)
295  C.data()[i]-=B.data()[i];
296  #endif
297  return C;
298  }
299 
300  inline SymMatrix SymMatrix::inverse() const {
301  #ifdef HAVE_LAPACK
302  SymMatrix invA(*this, DEEP_COPY);
303  // LU
304  BLAS_INT *pivots = new BLAS_INT[nlin()];
305  int Info = 0;
306  DSPTRF('U', sizet_to_int(nlin()), invA.data(), pivots, Info);
307  // Inverse
308  double *work = new double[this->nlin() * 64];
309  DSPTRI('U', sizet_to_int(nlin()), invA.data(), pivots, work, Info);
310  om_assert(Info==0);
311 
312  delete[] pivots;
313  delete[] work;
314  return invA;
315  #else
316  std::cerr << "!!!!! Inverse not implemented !!!!!" << std::endl;
317  exit(1);
318  #endif
319  }
320 
321  inline void SymMatrix::invert() {
322  #ifdef HAVE_LAPACK
323  // LU
324  BLAS_INT *pivots = new BLAS_INT[nlin()];
325  int Info = 0;
326  DSPTRF('U', sizet_to_int(nlin()), data(), pivots, Info);
327  // Inverse
328  double *work = new double[this->nlin() * 64];
329  DSPTRI('U', sizet_to_int(nlin()), data(), pivots, work, Info);
330 
331  om_assert(Info==0);
332  delete[] pivots;
333  delete[] work;
334  return;
335  #else
336  std::cerr << "!!!!! Inverse not implemented !!!!!" << std::endl;
337  exit(1);
338  #endif
339  }
340 
341  inline Vector SymMatrix::operator *(const Vector &v) const {
342  om_assert(nlin()==v.size());
343  Vector y(nlin());
344  #ifdef HAVE_BLAS
345  DSPMV(CblasUpper,sizet_to_int(nlin()),1.,data(),v.data(),1,0.,y.data(),1);
346  #else
347  for (size_t i=0;i<nlin();i++) {
348  y(i)=0;
349  for (size_t j=0;j<nlin();j++)
350  y(i)+=(*this)(i,j)*v(j);
351  }
352  #endif
353  return y;
354  }
355 
356  inline Vector SymMatrix::getlin(size_t i) const {
357  om_assert(i<nlin());
358  Vector v(ncol());
359  for ( size_t j = 0; j < ncol(); ++j) v.data()[j] = this->operator()(i,j);
360  return v;
361  }
362 
363  inline void SymMatrix::setlin(size_t i,const Vector& v) {
364  om_assert(v.size()==nlin() && i<nlin());
365  for ( size_t j = 0; j < ncol(); ++j) this->operator()(i,j) = v(j);
366  }
367 }
SymMatrix operator/(double x) const
Definition: symmatrix.h:104
SymMatrix operator-(const SymMatrix &B) const
Definition: symmatrix.h:287
double * data() const
Definition: vector.h:77
SymMatrix(const char *fname)
Definition: symmatrix.h:63
SymMatrix(size_t M, size_t N)
Definition: symmatrix.h:65
Vector getlin(size_t i) const
Definition: symmatrix.h:356
void operator+=(const SymMatrix &B)
Definition: symmatrix.h:193
bool empty() const
Definition: symmatrix.h:80
size_t size() const
Definition: vector.h:73
size_t size() const
Definition: symmatrix.h:71
SymMatrix(const SymMatrix &S, const DeepCopy)
Definition: symmatrix.h:66
utils::RCPtr< LinOpValue > value
Definition: symmatrix.h:57
size_t ncol() const
Definition: symmatrix.h:74
void save(const std::string &s) const
Definition: symmatrix.h:119
SymMatrix operator+(const SymMatrix &B) const
Definition: symmatrix.h:275
Vect3 operator*(const double &d, const Vect3 &v)
Definition: vect3.h:151
void operator-=(const SymMatrix &B)
Definition: symmatrix.h:183
void setlin(size_t i, const Vector &v)
Definition: symmatrix.h:363
size_t & ncol()
Definition: symmatrix.h:75
SymMatrix operator*(const SymMatrix &B) const
OPENMEEGMATHS_EXPORT BLAS_INT sizet_to_int(const size_t &num)
Definition: linop.h:56
SymMatrix(size_t N)
Definition: symmatrix.h:64
size_t nlin() const
Definition: linop.h:85
double * data() const
Definition: symmatrix.h:82
Vector solveLin(const Vector &B) const
Definition: symmatrix.h:142
void load(const std::string &s)
Definition: symmatrix.h:120
int BLAS_INT
SymMatrix posdefinverse() const
Definition: symmatrix.h:203
DeepCopy
Definition: linop.h:126
void reference_data(const double *array)
Definition: symmatrix.h:78
SymMatrix inverse() const
Definition: symmatrix.h:300
Matrix class.
Definition: matrix.h:61
double operator()(size_t i, size_t j) const
Definition: symmatrix.h:125