OpenMEEG
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
matrix.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 <OpenMEEGMathsConfig.h>
43 #include <iostream>
44 #include <cstdlib>
45 #include <string>
46 
47 #include <linop.h>
48 #include <MathsIO.H>
49 #include <symmatrix.h>
50 
51 namespace OpenMEEG {
52 
53  class SparseMatrix;
54  class SymMatrix;
55  class Vector;
56 
61  class OPENMEEGMATHS_EXPORT Matrix: public LinOp {
62  protected:
63 
64  friend class Vector;
65 
66  utils::RCPtr<LinOpValue> value;
67 
68  explicit Matrix(const Matrix& A,const size_t M): LinOp(A.nlin(),M,FULL,2),value(A.value) { }
69 
70  public:
71 
72  Matrix(): LinOp(0,0,FULL,2),value() { }
73  Matrix(const char* fname): LinOp(0,0,FULL,2),value() { this->load(fname); }
74  Matrix(const size_t M,const size_t N): LinOp(M,N,FULL,2),value(new LinOpValue(N*M)) { }
75  Matrix(const Matrix& A,const DeepCopy): LinOp(A.nlin(),A.ncol(),FULL,2),value(new LinOpValue(A.size(),A.data())) { }
76 
77  explicit Matrix(const SymMatrix& A);
78  explicit Matrix(const SparseMatrix& A);
79 
80  Matrix(const Vector& v,const size_t M,const size_t N);
81 
82  void alloc_data() { value = new LinOpValue(size()); }
83  void reference_data(const double* vals) { value = new LinOpValue(size(),vals); }
84 
89  bool empty() const { return value->empty(); }
90 
95  size_t size() const { return nlin()*ncol(); };
96 
101  double* data() const { return value->data; }
102 
107  inline double operator()(size_t i,size_t j) const ;
108 
113  inline double& operator()(size_t i,size_t j) ;
114 
115  Matrix submat(size_t istart, size_t isize, size_t jstart, size_t jsize) const;
116  void insertmat(size_t istart, size_t jstart, const Matrix& B);
117  Vector getcol(size_t j) const;
118  void setcol(size_t j, const Vector& v);
119  Vector getlin(size_t i) const;
120  void setlin(size_t i, const Vector& v);
121 
122  const Matrix& set(const double d);
123 
124  Matrix operator*(const Matrix& B) const;
125  Matrix operator*(const SymMatrix& B) const;
126  Matrix operator*(const SparseMatrix& B) const;
127  Matrix operator+(const Matrix& B) const;
128  Matrix operator-(const Matrix& B) const;
129  Matrix operator*(double x) const;
130  Matrix operator/(double x) const;
131  void operator+=(const Matrix& B);
132  void operator-=(const Matrix& B);
133  void operator*=(double x);
134  void operator/=(double x);
135 
136  Vector operator*(const Vector& v) const;
137  Vector tmult(const Vector &v) const;
138  Matrix tmult(const Matrix &m) const;
139  Matrix multt(const Matrix &m) const;
140  Matrix tmultt(const Matrix &m) const;
141 
142  Vector mean() const;
143  Vector tmean() const;
144 
145  Matrix transpose() const;
146  Matrix inverse() const;
147  Matrix pinverse(double reltol=0) const;
148  void svd(Matrix &U, SparseMatrix &S, Matrix &V, bool complete=true) const;
149 
154  double frobenius_norm() const;
155  double dot(const Matrix& B) const;
156 
160  void save(const char *filename) const;
161 
165  void load(const char *filename);
166 
167  void save(const std::string& s) const { save(s.c_str()); }
168  void load(const std::string& s) { load(s.c_str()); }
169 
173  void info() const;
174 
175  friend class SparseMatrix;
176  friend class SymMatrix;
177  };
178 
179  inline double Matrix::operator()(size_t i,size_t j) const {
180  om_assert(i<nlin() && j<ncol());
181  return value->data[i+nlin()*j];
182  }
183  inline double& Matrix::operator()(size_t i,size_t j) {
184  om_assert(i<nlin() && j<ncol());
185  return value->data[i+nlin()*j];
186  }
187 
188  inline double Matrix::frobenius_norm() const {
189  #ifdef HAVE_LAPACK
190  if ( nlin()*ncol() != 0 ) {
191  double work;
192  return DLANGE('F', sizet_to_int(nlin()), sizet_to_int(ncol()), data(), sizet_to_int(nlin()), &work);
193  } else {
194  return 0;
195  }
196  #else
197  double d = 0.;
198  for (size_t i=0; i<nlin()*ncol(); i++) d+=data()[i]*data()[i];
199  return sqrt(d);
200  #endif
201  }
202 
203  inline Vector Matrix::operator*(const Vector &v) const {
204  om_assert(ncol()==v.nlin());
205  Vector y(nlin());
206  #ifdef HAVE_BLAS
207  DGEMV(CblasNoTrans,sizet_to_int(nlin()),sizet_to_int(ncol()),1.0,data(),sizet_to_int(nlin()),v.data(),1,0.,y.data(),1);
208  #else
209  for (size_t i=0;i<nlin();i++) {
210  y(i) = 0;
211  for (size_t j=0;j<ncol();j++)
212  y(i) += (*this)(i,j)*v(j);
213  }
214  #endif
215 
216  return y;
217  }
218 
219  inline Matrix Matrix::submat(size_t istart, size_t isize, size_t jstart, size_t jsize) const {
220  om_assert (istart+isize<=nlin() && jstart+jsize<=ncol());
221 
222  Matrix a(isize,jsize);
223  const int sz = sizet_to_int(isize);
224 
225  for (size_t j=0; j<jsize; j++) {
226  #ifdef HAVE_BLAS
227  BLAS(dcopy,DCOPY)(sz,data()+istart+(jstart+j)*nlin(),1,a.data()+j*isize,1);
228  #else
229  for (size_t i=0; i<isize; i++)
230  a(i,j) = (*this)(istart+i,jstart+j);
231  #endif
232  }
233  return a;
234  }
235 
236  inline void Matrix::insertmat(size_t istart, size_t jstart, const Matrix& B) {
237  om_assert (istart+B.nlin()<=nlin() && jstart+B.ncol()<=ncol() );
238  for (size_t j=0; j<B.ncol(); j++) {
239  for (size_t i=0; i<B.nlin(); i++) {
240  (*this)(istart+i,jstart+j)=B(i,j);
241  }
242  }
243  }
244 
245  inline Vector Matrix::getcol(size_t j) const {
246  om_assert(j<ncol());
247  Vector v(nlin());
248  #ifdef HAVE_BLAS
249  BLAS(dcopy,DCOPY)(sizet_to_int(nlin()),data()+nlin()*j,1,v.data(),1);
250  #else
251  for (size_t i=0;i<nlin();i++) v.data()[i]=data()[i+nlin()*j];
252  #endif
253  return v;
254  }
255 
256  inline Vector Matrix::getlin(size_t i) const {
257  om_assert(i<nlin());
258  Vector v(ncol());
259  #ifdef HAVE_BLAS
260  BLAS(dcopy,DCOPY)(sizet_to_int(ncol()),data()+i,sizet_to_int(nlin()),v.data(),1);
261  #else
262  for (size_t j=0;j<ncol();j++) v.data()[j]=data()[i+nlin()*j];
263  #endif
264  return v;
265  }
266 
267  inline void Matrix::setcol(size_t j,const Vector& v) {
268  om_assert(v.size()==nlin() && j<ncol());
269  #ifdef HAVE_BLAS
270  BLAS(dcopy,DCOPY)(sizet_to_int(nlin()),v.data(),1,data()+nlin()*j,1);
271  #else
272  for (size_t i=0;i<nlin();i++) data()[i+nlin()*j]=v.data()[i];
273  #endif
274  }
275 
276  inline void Matrix::setlin(size_t i,const Vector& v) {
277  om_assert(v.size()==ncol() && i<nlin());
278  #ifdef HAVE_BLAS
279  BLAS(dcopy,DCOPY)(sizet_to_int(ncol()),v.data(),1,data()+i,sizet_to_int(nlin()));
280  #else
281  for (size_t j=0;j<ncol();j++) data()[i+nlin()*j]=v.data()[j];
282  #endif
283  }
284 
285  inline Vector Matrix::tmult(const Vector &v) const {
286  om_assert(nlin()==v.nlin());
287  Vector y(ncol());
288  #ifdef HAVE_BLAS
289  DGEMV(CblasTrans,sizet_to_int(nlin()),sizet_to_int(ncol()),1.,data(),sizet_to_int(nlin()),v.data(),1,0.,y.data(),1);
290  #else
291  for (size_t i=0;i<ncol();i++) {
292  y(i)=0;
293  for (size_t j=0;j<nlin();j++)
294  y(i)+=(*this)(j,i)*v(j);
295  }
296  #endif
297 
298  return y;
299  }
300 
301  inline Matrix Matrix::inverse() const {
302  #ifdef HAVE_LAPACK
303  om_assert(nlin()==ncol());
304  Matrix invA(*this,DEEP_COPY);
305  // LU
306  #if defined(CLAPACK_INTERFACE)
307  #if defined(__APPLE__) && defined(USE_VECLIB) // Apple Veclib Framework (Handles 32 and 64 Bits)
308  __CLPK_integer *pivots = new __CLPK_integer[ncol()];
309  __CLPK_integer Info = 0;
310  __CLPK_integer nlin_local = invA.nlin();
311  __CLPK_integer nlin_local2 = invA.nlin();
312  __CLPK_integer ncol_local = invA.ncol();
313  __CLPK_integer sz = invA.ncol()*64;
314  DGETRF(nlin_local,ncol_local,invA.data(),nlin_local2,pivots,Info);
315  double *work=new double[sz];
316  DGETRI(ncol_local,invA.data(),ncol_local,pivots,work,sz,Info);
317  delete[] pivots;
318  delete[] work;
319  om_assert(Info==0);
320  #else
321  BLAS_INT *pivots=new BLAS_INT[sizet_to_int(ncol())];
322  DGETRF(sizet_to_int(invA.nlin()),sizet_to_int(invA.ncol()),invA.data(),sizet_to_int(invA.nlin()),pivots);
323  DGETRI(sizet_to_int(invA.ncol()),invA.data(),sizet_to_int(invA.ncol()),pivots);
324  delete[] pivots;
325  #endif
326  #else
327  int Info = 0;
328  int *pivots=new int[sizet_to_int(ncol())];
329  DGETRF(sizet_to_int(invA.nlin()),sizet_to_int(invA.ncol()),invA.data(),sizet_to_int(invA.nlin()),pivots,Info);
330  const unsigned sz = invA.ncol()*64;
331  double *work=new double[sz];
332  DGETRI(sizet_to_int(invA.ncol()),invA.data(),sizet_to_int(invA.ncol()),pivots,work,sz,Info);
333  delete[] pivots;
334  delete[] work;
335  om_assert(Info==0);
336  #endif
337  return invA;
338  #else
339  std::cerr << "!!!!! Inverse not implemented !!!!!" << std::endl;
340  exit(1);
341  #endif
342  }
343 
344  inline Matrix Matrix::operator *(const Matrix &B) const {
345  om_assert(ncol()==B.nlin());
346  size_t p=ncol();
347  Matrix C(nlin(),B.ncol());
348  #ifdef HAVE_BLAS
349  DGEMM(CblasNoTrans,CblasNoTrans,
350  sizet_to_int(C.nlin()),sizet_to_int(C.ncol()),sizet_to_int(p),
351  1.,data(),sizet_to_int(nlin()),
352  B.data(),sizet_to_int(B.nlin()),
353  0.,C.data(),sizet_to_int(C.nlin()));
354  #else
355  for (size_t i=0;i<C.nlin();i++)
356  for (size_t j=0;j<C.ncol();j++) {
357  C(i,j)=0;
358  for (size_t k=0;k<p;k++)
359  C(i,j)+=(*this)(i,k)*B(k,j);
360  }
361  #endif
362  return C;
363  }
364 
365  inline Matrix Matrix::tmult(const Matrix &B) const {
366  om_assert(nlin()==B.nlin());
367  size_t p=nlin();
368  Matrix C(ncol(),B.ncol());
369  #ifdef HAVE_BLAS
370  DGEMM(CblasTrans,CblasNoTrans,
371  sizet_to_int(C.nlin()),sizet_to_int(C.ncol()),sizet_to_int(p),
372  1.,data(),sizet_to_int(nlin()),
373  B.data(),sizet_to_int(B.nlin()),
374  0.,C.data(),sizet_to_int(C.nlin()));
375  #else
376  for (size_t i=0;i<C.nlin();i++)
377  for (size_t j=0;j<C.ncol();j++) {
378  C(i,j)=0;
379  for (size_t k=0;k<p;k++)
380  C(i,j)+=(*this)(k,i)*B(k,j);
381  }
382  #endif
383  return C;
384  }
385 
386  inline Matrix Matrix::multt(const Matrix &B) const {
387  om_assert(ncol()==B.ncol());
388  size_t p=ncol();
389  Matrix C(nlin(),B.nlin());
390  #ifdef HAVE_BLAS
391  DGEMM(CblasNoTrans,CblasTrans,
392  sizet_to_int(C.nlin()),sizet_to_int(C.ncol()),sizet_to_int(p),
393  1.,data(),sizet_to_int(nlin()),
394  B.data(),sizet_to_int(B.nlin()),
395  0.,C.data(),sizet_to_int(C.nlin()));
396  #else
397  for (size_t i=0;i<C.nlin();i++)
398  for (size_t j=0;j<C.ncol();j++) {
399  C(i,j)=0;
400  for (size_t k=0;k<p;k++)
401  C(i,j)+=(*this)(i,k)*B(j,k);
402  }
403  #endif
404  return C;
405  }
406 
407  inline Matrix Matrix::tmultt(const Matrix &B) const {
408  om_assert(nlin()==B.ncol());
409  size_t p=nlin();
410  Matrix C(ncol(),B.nlin());
411  #ifdef HAVE_BLAS
412  DGEMM(CblasTrans,CblasTrans,
413  sizet_to_int(C.nlin()),sizet_to_int(C.ncol()),sizet_to_int(p),
414  1.,data(),sizet_to_int(nlin()),
415  B.data(),sizet_to_int(B.nlin()),
416  0.,C.data(),sizet_to_int(C.nlin()));
417  #else
418  for (size_t i=0;i<C.nlin();i++)
419  for (size_t j=0;j<C.ncol();j++) {
420  C(i,j)=0;
421  for (size_t k=0;k<p;k++)
422  C(i,j)+=(*this)(k,i)*B(j,k);
423  }
424  #endif
425  return C;
426  }
427 
428  inline Matrix Matrix::operator*(const SymMatrix &B) const {
429  om_assert(ncol()==B.ncol());
430  Matrix C(nlin(),B.ncol());
431 
432  #ifdef HAVE_BLAS
433  Matrix D(B);
434  const int n = sizet_to_int(nlin());
435  const int m = sizet_to_int(D.ncol());
436  const int l = sizet_to_int(C.nlin());
437  DSYMM(CblasRight,CblasUpper,n,m,1.,D.data(),m,data(),n,0.,C.data(),l);
438  #else
439  for (size_t j=0;j<B.ncol();j++)
440  for (size_t i=0;i<ncol();i++) {
441  C(i,j)=0;
442  for (size_t k=0;k<ncol();k++)
443  C(i,j)+=(*this)(i,k)*B(k,j);
444  }
445  #endif
446  return C;
447  }
448 
449  inline Matrix Matrix::operator+(const Matrix &B) const {
450  om_assert(ncol()==B.ncol());
451  om_assert(nlin()==B.nlin());
452  Matrix C(*this,DEEP_COPY);
453  #ifdef HAVE_BLAS
454  BLAS(daxpy,DAXPY)(sizet_to_int(nlin()*ncol()), 1.0, B.data(), 1, C.data() , 1);
455  #else
456  for (size_t i=0;i<nlin()*ncol();i++)
457  C.data()[i]+=B.data()[i];
458  #endif
459  return C;
460  }
461 
462  inline Matrix Matrix::operator-(const Matrix &B) const {
463  om_assert(ncol()==B.ncol());
464  om_assert(nlin()==B.nlin());
465  Matrix C(*this,DEEP_COPY);
466  #ifdef HAVE_BLAS
467  BLAS(daxpy,DAXPY)(sizet_to_int(nlin()*ncol()), -1.0, B.data(), 1, C.data() , 1);
468  #else
469  for (size_t i=0;i<nlin()*ncol();i++)
470  C.data()[i]-=B.data()[i];
471  #endif
472  return C;
473  }
474 
475  inline void Matrix::operator+=(const Matrix &B) {
476  om_assert(ncol()==B.ncol());
477  om_assert(nlin()==B.nlin());
478  #ifdef HAVE_BLAS
479  BLAS(daxpy,DAXPY)(sizet_to_int(nlin()*ncol()), 1.0, B.data(), 1, data() , 1);
480  #else
481  for (size_t i=0;i<nlin()*ncol();i++)
482  data()[i]+=B.data()[i];
483  #endif
484  }
485 
486  inline void Matrix::operator-=(const Matrix &B) {
487  om_assert(ncol()==B.ncol());
488  om_assert(nlin()==B.nlin());
489  #ifdef HAVE_BLAS
490  BLAS(daxpy,DAXPY)(sizet_to_int(nlin()*ncol()), -1.0, B.data(), 1, data() , 1);
491  #else
492  for (size_t i=0;i<nlin()*ncol();i++)
493  data()[i]-=B.data()[i];
494  #endif
495  }
496 
497  inline double Matrix::dot(const Matrix& b) const {
498  om_assert(nlin()==b.nlin()&&ncol()==b.ncol());
499  #ifdef HAVE_BLAS
500  return BLAS(ddot,DDOT)(sizet_to_int(nlin()*ncol()),data(),1,b.data(),1);
501  #else
502  double s=0;
503  for (size_t i=0;i<nlin()*ncol();i++)
504  s+=data()[i]*b.data()[i];
505  return s;
506  #endif
507  }
508 }
size_t size() const
Get Matrix size.
Definition: matrix.h:95
void alloc_data()
Definition: matrix.h:82
double * data() const
Definition: vector.h:77
Matrix operator*(const Matrix &B) const
Definition: matrix.h:344
Matrix(const char *fname)
Definition: matrix.h:73
utils::RCPtr< LinOpValue > value
Definition: matrix.h:66
Matrix(const Matrix &A, const DeepCopy)
Definition: matrix.h:75
size_t size() const
Definition: vector.h:73
Matrix(const size_t M, const size_t N)
Definition: matrix.h:74
Matrix submat(size_t istart, size_t isize, size_t jstart, size_t jsize) const
Definition: matrix.h:219
void load(const std::string &s)
Definition: matrix.h:168
double frobenius_norm() const
Get Matrix Frobenius norm.
Definition: matrix.h:188
size_t ncol() const
Definition: symmatrix.h:74
double operator()(size_t i, size_t j) const
Get Matrix value.
Definition: matrix.h:179
Matrix inverse() const
Definition: matrix.h:301
Vect3 operator*(const double &d, const Vect3 &v)
Definition: vect3.h:151
double dot(const Matrix &B) const
Definition: matrix.h:497
Matrix multt(const Matrix &m) const
Definition: matrix.h:386
void setcol(size_t j, const Vector &v)
Definition: matrix.h:267
Matrix(const Matrix &A, const size_t M)
Definition: matrix.h:68
virtual size_t ncol() const
Definition: linop.h:80
Vector getlin(size_t i) const
Definition: matrix.h:256
Matrix operator-(const Matrix &B) const
Definition: matrix.h:462
double * data() const
Get Matrix data.
Definition: matrix.h:101
OPENMEEGMATHS_EXPORT BLAS_INT sizet_to_int(const size_t &num)
Definition: linop.h:56
Vector getcol(size_t j) const
Definition: matrix.h:245
void reference_data(const double *vals)
Definition: matrix.h:83
void save(const std::string &s) const
Definition: matrix.h:167
size_t nlin() const
Definition: linop.h:77
Matrix operator+(const Matrix &B) const
Definition: matrix.h:449
Vector tmult(const Vector &v) const
Definition: matrix.h:285
int BLAS_INT
bool empty() const
Test if Matrix is empty.
Definition: matrix.h:89
void operator+=(const Matrix &B)
Definition: matrix.h:475
DeepCopy
Definition: linop.h:113
void insertmat(size_t istart, size_t jstart, const Matrix &B)
Definition: matrix.h:236
Matrix class.
Definition: matrix.h:61
void operator-=(const Matrix &B)
Definition: matrix.h:486
Matrix tmultt(const Matrix &m) const
Definition: matrix.h:407
void setlin(size_t i, const Vector &v)
Definition: matrix.h:276