OpenMEEG
matrix.h
Go to the documentation of this file.
1 // Project Name: OpenMEEG (http://openmeeg.github.io)
2 // © INRIA and ENPC under the French open source license CeCILL-B.
3 // See full copyright notice in the file LICENSE.txt
4 // If you make a copy of this file, you must either:
5 // - provide also LICENSE.txt and modify this header to refer to it.
6 // - replace this header by the LICENSE.txt content.
7 
8 #pragma once
9 
10 #include <OpenMEEGMathsConfig.h>
11 #include <iostream>
12 #include <cstdlib>
13 #include <string>
14 
15 #include <linop.h>
16 #include <MathsIO.H>
17 #include <symmatrix.h>
18 
19 namespace OpenMEEG {
20 
21  class SparseMatrix;
22  class SymMatrix;
23  class Vector;
24 
27 
28  class OPENMEEGMATHS_EXPORT Matrix: public LinOp {
29  protected:
30 
31  friend class Vector;
32 
34 
35  explicit Matrix(const Matrix& A,const Dimension M): LinOp(A.nlin(),M,FULL,2),value(A.value) { }
36 
37  public:
38 
39  Matrix(): LinOp(0,0,FULL,2),value() { }
40  Matrix(const char* fname): LinOp(0,0,FULL,2),value() { load(fname); }
41  Matrix(const std::string& fname): Matrix(fname.c_str()) { }
42  Matrix(const Dimension M,const Dimension N): LinOp(M,N,FULL,2),value(N*M) { }
43  Matrix(const Matrix& A,const DeepCopy): LinOp(A.nlin(),A.ncol(),FULL,2),value(A.size(),A.data()) { }
44 
45  explicit Matrix(const SymMatrix& A);
46  explicit Matrix(const SparseMatrix& A);
47 
48  Matrix(const Vector& v,const Dimension M,const Dimension N);
49 
50  void alloc_data() { value = LinOpValue(size()); }
51  void reference_data(const double* vals) { value = LinOpValue(size(),vals); }
52 
55 
56  bool empty() const { return value.empty(); }
57 
60 
61  size_t size() const { return static_cast<std::streamoff>(nlin())*static_cast<std::streamoff>(ncol()); };
62 
65 
66  double* data() const { return value.get(); }
67 
70 
71  double operator()(const Index i,const Index j) const {
72  om_assert(i<nlin() && j<ncol());
73  return value[i+nlin()*j];
74  }
75 
78 
79  double& operator()(const Index i,const Index j) {
80  om_assert(i<nlin() && j<ncol());
81  return value[i+nlin()*j];
82  }
83 
84  Matrix submat(const Index istart,const Index isize,const Index jstart,const Index jsize) const;
85  void insertmat(const Index istart,const Index jstart,const Matrix& B);
86 
87  Vector getcol(const Index j) const;
88  void setcol(const Index j,const Vector& v);
89 
90  Vector getlin(const Index i) const;
91  void setlin(const Index i,const Vector& v);
92 
93  const Matrix& set(const double d);
94 
95  Matrix operator*(const Matrix& B) const;
96  Matrix operator*(const SymMatrix& B) const;
97  Matrix operator*(const SparseMatrix& B) const;
98 
99  Matrix operator+(const Matrix& B) const {
100  Matrix C(*this,DEEP_COPY);
101  C += B;
102  return C;
103  }
104 
105  Matrix operator-(const Matrix& B) const {
106  Matrix C(*this,DEEP_COPY);
107  C -= B;
108  return C;
109  }
110 
111  Matrix operator*(double x) const;
112  Matrix operator/(double x) const;
113  inline void operator+=(const Matrix& B);
114  inline void operator-=(const Matrix& B);
115  void operator*=(double x);
116  void operator/=(double x);
117 
118  Vector operator*(const Vector& v) const;
119  Vector tmult(const Vector& v) const;
120  Matrix tmult(const Matrix& m) const;
121  Matrix multt(const Matrix& m) const;
122  Matrix tmultt(const Matrix& m) const;
123 
124  Matrix transpose() const;
125  Matrix inverse() const;
126  Matrix pinverse(const double reltol=0.0) const;
127  void svd(Matrix& U,SparseMatrix& S,Matrix& V,const bool complete=true) const;
128 
131 
132  double frobenius_norm() const;
133  double dot(const Matrix& B) const;
134 
136 
137  void save(const char* filename) const;
138 
140 
141  void load(const char* filename);
142 
143  void save(const std::string& s) const { save(s.c_str()); }
144  void load(const std::string& s) { load(s.c_str()); }
145 
147 
148  void info() const;
149 
150  friend class SparseMatrix;
151  friend class SymMatrix;
152  };
153 
154  inline std::ostream& operator<<(std::ostream& os,const Matrix& M) {
155  for (unsigned i=0; i<M.nlin(); ++i) {
156  for (unsigned j=0; j<M.ncol(); ++j)
157  os << M(i,j) << ' ';
158  os << std::endl;
159  }
160  return os;
161  }
162 
163  inline double Matrix::frobenius_norm() const {
164  const size_t sz = size();
165  if (sz==0)
166  return 0.0;
167 
168  #ifdef HAVE_LAPACK
169  double work;
170  const BLAS_INT M = sizet_to_int(nlin());
171  const BLAS_INT N = sizet_to_int(ncol());
172  return DLANGE('F',M,N,data(),M,&work);
173  #else
174  double d = 0.0;
175  for (size_t i=0; i<sz; i++)
176  d += data()[i]*data()[i];
177  return sqrt(d);
178  #endif
179  }
180 
181  inline Vector Matrix::operator*(const Vector& v) const {
182  om_assert(ncol()==v.nlin());
183  Vector res(nlin());
184  #ifdef HAVE_BLAS
185  const BLAS_INT M = sizet_to_int(nlin());
186  const BLAS_INT N = sizet_to_int(ncol());
187  DGEMV(CblasNoTrans,M,N,1.0,data(),M,v.data(),1,0.0,res.data(),1);
188  #else
189  res.set(0.0);
190  for (Index j=0; j<ncol(); ++j)
191  for (Index i=0; i<nlin(); ++i)
192  res(i) += (*this)(i,j)*v(j);
193  #endif
194 
195  return res;
196  }
197 
198  inline Matrix Matrix::submat(const Index istart,const Index isize,const Index jstart,const Index jsize) const {
199  om_assert(istart+isize<=nlin() && jstart+jsize<=ncol());
200 
201  Matrix res(isize,jsize);
202 
203  for (Index j=0; j<jsize; ++j) {
204  #ifdef HAVE_BLAS
205  const BLAS_INT sz = sizet_to_int(isize);
206  BLAS(dcopy,DCOPY)(sz,data()+istart+(jstart+j)*nlin(),1,res.data()+j*isize,1);
207  #else
208  for (Index i=0; i<isize; ++i)
209  res(i,j) = (*this)(istart+i,jstart+j);
210  #endif
211  }
212  return res;
213  }
214 
215  inline void Matrix::insertmat(const Index istart,const Index jstart,const Matrix& B) {
216  om_assert(istart+B.nlin()<=nlin() && jstart+B.ncol()<=ncol() );
217 
218  for (Index j=0; j<B.ncol(); ++j)
219  for (Index i=0; i<B.nlin(); ++i)
220  (*this)(istart+i,jstart+j) = B(i,j);
221  }
222 
223  inline Vector Matrix::getcol(const Index j) const {
224  om_assert(j<ncol( ));
225  Vector res(nlin());
226  #ifdef HAVE_BLAS
227  const BLAS_INT M = sizet_to_int(nlin());
228  BLAS(dcopy,DCOPY)(M,data()+nlin()*j,1,res.data(),1);
229  #else
230  for (Index i=0; i<nlin(); ++i)
231  res(i) = (*this)(i,j);
232  #endif
233  return res;
234  }
235 
236  inline Vector Matrix::getlin(const Index i) const {
237  om_assert(i<nlin());
238  Vector res(ncol());
239  #ifdef HAVE_BLAS
240  const BLAS_INT M = sizet_to_int(nlin());
241  const BLAS_INT N = sizet_to_int(ncol());
242  BLAS(dcopy,DCOPY)(N,data()+i,M,res.data(),1);
243  #else
244  for (Index j=0; j<ncol(); ++j)
245  res(j) = (*this)(i,j);
246  #endif
247  return res;
248  }
249 
250  inline void Matrix::setcol(const Index j,const Vector& v) {
251  om_assert(v.size()==nlin() && j<ncol());
252  #ifdef HAVE_BLAS
253  const BLAS_INT M = sizet_to_int(nlin());
254  BLAS(dcopy,DCOPY)(M,v.data(),1,data()+nlin()*j,1);
255  #else
256  for (Index i=0; i<nlin(); ++i)
257  (*this)(i,j) = v(i);
258  #endif
259  }
260 
261  inline void Matrix::setlin(const Index i,const Vector& v) {
262  om_assert(v.size()==ncol());
263  om_assert(i<nlin());
264  #ifdef HAVE_BLAS
265  const BLAS_INT M = sizet_to_int(nlin());
266  const BLAS_INT N = sizet_to_int(ncol());
267  BLAS(dcopy,DCOPY)(N,v.data(),1,data()+i,M);
268  #else
269  for (Index j=0; j<ncol(); ++j)
270  (*this)(i,j) = v(j);
271  #endif
272  }
273 
274  inline Vector Matrix::tmult(const Vector& v) const {
275  om_assert(nlin()==v.nlin());
276  Vector res(ncol());
277  #ifdef HAVE_BLAS
278  const BLAS_INT M = sizet_to_int(nlin());
279  const BLAS_INT N = sizet_to_int(ncol());
280  DGEMV(CblasTrans,M,N,1.0,data(),sizet_to_int(nlin()),v.data(),1,0.0,res.data(),1);
281  #else
282  for (Index i=0; i<ncol(); ++i) {
283  res(i) = 0;
284  for (Index j=0; j<nlin(); ++j)
285  res(i) += (*this)(j,i)*v(j);
286  }
287  #endif
288 
289  return res;
290  }
291 
292  inline Matrix Matrix::inverse() const {
293  om_assert(nlin()==ncol());
294  #ifdef HAVE_LAPACK
295  Matrix invA(*this,DEEP_COPY);
296  // LU
297  #if defined(CLAPACK_INTERFACE)
298  const BLAS_INT M = sizet_to_int(invA.nlin());
299  const BLAS_INT N = sizet_to_int(ncol());
300  BLAS_INT* pivots = new BLAS_INT[N];
301  DGETRF(M,N,invA.data(),M,pivots);
302  DGETRI(N,invA.data(),N,pivots);
303  delete[] pivots;
304  #else
305  int Info = 0;
306  BLAS_INT M = sizet_to_int(invA.nlin());
307  BLAS_INT N = sizet_to_int(ncol());
308  int* pivots = new int[N];
309  DGETRF(M,N,invA.data(),M,pivots,Info);
310  const Dimension sz = invA.ncol()*64;
311  double* work=new double[sz];
312  DGETRI(N,invA.data(),N,pivots,work,sz,Info);
313  delete[] pivots;
314  delete[] work;
315  om_assert(Info==0);
316  #endif
317  return invA;
318  #else
319  throw OpenMEEG::maths::LinearAlgebraError("Inverse not implemented, requires LAPACK");
320  #endif
321  }
322 
323  inline Matrix Matrix::operator*(const Matrix& B) const {
324  om_assert(ncol()==B.nlin());
325  Matrix C(nlin(),B.ncol());
326  #ifdef HAVE_BLAS
327  const BLAS_INT M = sizet_to_int(nlin());
328  const BLAS_INT N = sizet_to_int(ncol());
329  const BLAS_INT L = sizet_to_int(B.ncol());
330  DGEMM(CblasNoTrans,CblasNoTrans,M,L,N,1.0,data(),M,B.data(),N,0.0,C.data(),M);
331  #else
332  for (Index i=0; i<C.nlin(); ++i)
333  for (Index j=0; j<C.ncol(); ++j) {
334  C(i,j) = 0.0;
335  for (Index k=0; k<ncol(); ++k)
336  C(i,j) += (*this)(i,k)*B(k,j);
337  }
338  #endif
339  return C;
340  }
341 
342  inline Matrix Matrix::tmult(const Matrix& B) const {
343  om_assert(nlin()==B.nlin());
344  Matrix C(ncol(),B.ncol());
345  #ifdef HAVE_BLAS
346  const BLAS_INT M = sizet_to_int(nlin());
347  const BLAS_INT N = sizet_to_int(ncol());
348  const BLAS_INT L = sizet_to_int(B.ncol());
349  DGEMM(CblasTrans,CblasNoTrans,N,L,M,1.0,data(),M,B.data(),M,0.0,C.data(),N);
350  #else
351  for (Index i=0; i<C.nlin(); ++i)
352  for (Index j=0; j<C.ncol(); ++j) {
353  C(i,j) = 0.0;
354  for (Index k=0; k<nlin(); ++k)
355  C(i,j) += (*this)(k,i)*B(k,j);
356  }
357  #endif
358  return C;
359  }
360 
361  inline Matrix Matrix::multt(const Matrix& B) const {
362  om_assert(ncol()==B.ncol());
363  Matrix C(nlin(),B.nlin());
364  #ifdef HAVE_BLAS
365  const BLAS_INT M = sizet_to_int(nlin());
366  const BLAS_INT N = sizet_to_int(ncol());
367  const BLAS_INT L = sizet_to_int(B.nlin());
368  DGEMM(CblasNoTrans,CblasTrans,M,L,N,1.0,data(),M,B.data(),L,0.0,C.data(),M);
369  #else
370  for (Index j=0; j<C.ncol(); ++j)
371  for (Index i=0; i<C.nlin(); ++i) {
372  C(i,j) = 0.0;
373  for (Index k=0; k<ncol(); ++k)
374  C(i,j) += (*this)(i,k)*B(j,k);
375  }
376  #endif
377  return C;
378  }
379 
380  inline Matrix Matrix::tmultt(const Matrix& B) const {
381  om_assert(nlin()==B.ncol());
382  Matrix C(ncol(),B.nlin());
383  #ifdef HAVE_BLAS
384  const BLAS_INT M = sizet_to_int(nlin());
385  const BLAS_INT N = sizet_to_int(ncol());
386  const BLAS_INT L = sizet_to_int(B.nlin());
387  DGEMM(CblasTrans,CblasTrans,L,N,M,1.0,data(),M,B.data(),N,0.0,C.data(),L);
388  #else
389  for (Index i=0; i<C.nlin(); ++i)
390  for (Index j=0; j<C.ncol(); ++j) {
391  C(i,j) = 0.0;
392  for (Index k=0; k<nlin(); ++k)
393  C(i,j) += (*this)(k,i)*B(j,k);
394  }
395  #endif
396  return C;
397  }
398 
399  inline Matrix Matrix::operator*(const SymMatrix& B) const {
400  om_assert(ncol()==B.nlin());
401  Matrix C(nlin(),B.ncol());
402 
403  // Workaround an MKL bug
404  //#ifdef HAVE_BLAS
405  #if defined(HAVE_BLAS) && !defined(USE_MKL)
406  Matrix D(B);
407  const BLAS_INT m = sizet_to_int(nlin());
408  const BLAS_INT n = sizet_to_int(B.ncol());
409  DSYMM(CblasRight,CblasUpper,m,n,1.0,D.data(),n,data(),m,0.0,C.data(),m);
410  #else
411  for (Index j=0; j<B.ncol(); ++j)
412  for (Index i=0; i<nlin(); ++i) {
413  double sum = 0.0;
414  for (size_t k=0; k<ncol(); ++k)
415  sum += (*this)(i,k)*B(k,j);
416  C(i,j) = sum;
417  }
418  #endif
419  return C;
420  }
421 
422  inline void Matrix::operator+=(const Matrix& B) {
423  om_assert(nlin()==B.nlin());
424  om_assert(ncol()==B.ncol());
425  #ifdef HAVE_BLAS
426  const BLAS_INT sz = sizet_to_int(size());
427  BLAS(daxpy,DAXPY)(sz,1.0,B.data(),1,data(),1);
428  #else
429  const size_t sz = size();
430  for (size_t i=0; i<sz; ++i)
431  data()[i] += B.data()[i];
432  #endif
433  }
434 
435  inline void Matrix::operator-=(const Matrix& B) {
436  om_assert(nlin()==B.nlin());
437  om_assert(ncol()==B.ncol());
438  #ifdef HAVE_BLAS
439  const BLAS_INT sz = sizet_to_int(size());
440  BLAS(daxpy,DAXPY)(sz,-1.0,B.data(),1,data(),1);
441  #else
442  const size_t sz = size();
443  for (size_t i=0; i<sz; ++i)
444  data()[i] -= B.data()[i];
445  #endif
446  }
447 
448  inline double Matrix::dot(const Matrix& B) const {
449  om_assert(nlin()==B.nlin());
450  om_assert(ncol()==B.ncol());
451  #ifdef HAVE_BLAS
452  const BLAS_INT sz = sizet_to_int(size());
453  return BLAS(ddot,DDOT)(sz,data(),1,B.data(),1);
454  #else
455  const sz = size();
456  double s = 0.0;
457  for (size_t i=0; i<sz; i++)
458  s += data()[i]*B.data()[i];
459  return s;
460  #endif
461  }
462 }
Dimension nlin() const
Definition: linop.h:48
virtual Dimension ncol() const
Definition: linop.h:51
Matrix class Matrix class.
Definition: matrix.h:28
void svd(Matrix &U, SparseMatrix &S, Matrix &V, const bool complete=true) const
void save(const char *filename) const
Save Matrix to file (Format set using file name extension)
void alloc_data()
Definition: matrix.h:50
Matrix(const SymMatrix &A)
size_t size() const
Get Matrix size.
Definition: matrix.h:61
Matrix(const Matrix &A, const DeepCopy)
Definition: matrix.h:43
const Matrix & set(const double d)
Matrix operator/(double x) const
Vector tmult(const Vector &v) const
Definition: matrix.h:274
Matrix(const char *fname)
Definition: matrix.h:40
double * data() const
Get Matrix data.
Definition: matrix.h:66
void operator/=(double x)
double & operator()(const Index i, const Index j)
Get Matrix value.
Definition: matrix.h:79
void save(const std::string &s) const
Definition: matrix.h:143
Vector getcol(const Index j) const
Definition: matrix.h:223
Matrix(const Vector &v, const Dimension M, const Dimension N)
Matrix(const Dimension M, const Dimension N)
Definition: matrix.h:42
Matrix operator-(const Matrix &B) const
Definition: matrix.h:105
double operator()(const Index i, const Index j) const
Get Matrix value.
Definition: matrix.h:71
Vector getlin(const Index i) const
Definition: matrix.h:236
Matrix multt(const Matrix &m) const
Definition: matrix.h:361
bool empty() const
Test if Matrix is empty.
Definition: matrix.h:56
LinOpValue value
Definition: matrix.h:33
Matrix(const Matrix &A, const Dimension M)
Definition: matrix.h:35
void info() const
Print info on Matrix.
Matrix submat(const Index istart, const Index isize, const Index jstart, const Index jsize) const
Definition: matrix.h:198
double frobenius_norm() const
Get Matrix Frobenius norm.
Definition: matrix.h:163
void setcol(const Index j, const Vector &v)
Definition: matrix.h:250
Matrix operator*(const SparseMatrix &B) const
void setlin(const Index i, const Vector &v)
Definition: matrix.h:261
double dot(const Matrix &B) const
Definition: matrix.h:448
Matrix transpose() const
void operator-=(const Matrix &B)
Definition: matrix.h:435
Matrix operator*(const Matrix &B) const
Definition: matrix.h:323
Matrix operator+(const Matrix &B) const
Definition: matrix.h:99
void load(const char *filename)
Load Matrix from file (Format set using file name extension)
void operator+=(const Matrix &B)
Definition: matrix.h:422
void load(const std::string &s)
Definition: matrix.h:144
Matrix(const SparseMatrix &A)
void operator*=(double x)
void insertmat(const Index istart, const Index jstart, const Matrix &B)
Definition: matrix.h:215
void reference_data(const double *vals)
Definition: matrix.h:51
Matrix pinverse(const double reltol=0.0) const
Matrix inverse() const
Definition: matrix.h:292
Matrix tmultt(const Matrix &m) const
Definition: matrix.h:380
Matrix operator*(double x) const
Matrix(const std::string &fname)
Definition: matrix.h:41
Dimension ncol() const
Definition: symmatrix.h:42
void set(const double x)
size_t size() const
Definition: vector.h:40
double * data() const
Definition: vector.h:44
Vect3 operator*(const double d, const Vect3 &V)
Definition: vect3.h:105
DeepCopy
Definition: linop.h:84
@ DEEP_COPY
Definition: linop.h:84
unsigned Dimension
Definition: linop.h:32
std::ostream & operator<<(std::ostream &os, const Conductivity< REP > &m)
unsigned Index
Definition: linop.h:33
BLAS_INT sizet_to_int(const unsigned &num)
Definition: linop.h:26
bool empty() const
Definition: linop.h:96