OpenMEEG
Loading...
Searching...
No Matches
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
19namespace 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
68 #if defined(SWIGPYTHON) || defined(WIN32)
69
70 // This method is only needed for swig wrapping. But windows insists on having it
71 // in the library despite the fact it is inlined.
72
74
75 std::shared_ptr<double[]> get_shared_data_ptr() { return value; }
76 #endif
77
80
81 double operator()(const Index i,const Index j) const {
82 om_assert(i<nlin() && j<ncol());
83 return value[i+nlin()*j];
84 }
85
88
89 double& operator()(const Index i,const Index j) {
90 om_assert(i<nlin() && j<ncol());
91 return value[i+nlin()*j];
92 }
93
94 Matrix submat(const Index istart,const Index isize,const Index jstart,const Index jsize) const;
95 void insertmat(const Index istart,const Index jstart,const Matrix& B);
96
97 Vector getcol(const Index j) const;
98 void setcol(const Index j,const Vector& v);
99
100 Vector getlin(const Index i) const;
101 void setlin(const Index i,const Vector& v);
102
103 const Matrix& set(const double d);
104
105 Matrix operator*(const Matrix& B) const;
106 Matrix operator*(const SymMatrix& B) const;
108
109 Matrix operator+(const Matrix& B) const {
110 Matrix C(*this,DEEP_COPY);
111 C += B;
112 return C;
113 }
114
115 Matrix operator-(const Matrix& B) const {
116 Matrix C(*this,DEEP_COPY);
117 C -= B;
118 return C;
119 }
120
121 Matrix operator*(double x) const;
122 Matrix operator/(double x) const;
123 inline void operator+=(const Matrix& B);
124 inline void operator-=(const Matrix& B);
125 void operator*=(double x);
126 void operator/=(double x);
127
128 Vector operator*(const Vector& v) const;
129 Vector tmult(const Vector& v) const;
130 Matrix tmult(const Matrix& m) const;
131 Matrix multt(const Matrix& m) const;
132 Matrix tmultt(const Matrix& m) const;
133
135 Matrix inverse() const;
136 Matrix pinverse(const double reltol=0.0) const;
137 void svd(Matrix& U,SparseMatrix& S,Matrix& V,const bool complete=true) const;
138
141
142 double frobenius_norm() const;
143 double dot(const Matrix& B) const;
144
146
147 void save(const char* filename) const;
148
150
151 void load(const char* filename);
152
153 void save(const std::string& s) const { save(s.c_str()); }
154 void load(const std::string& s) { load(s.c_str()); }
155
157
158 void info() const;
159
160 friend class SparseMatrix;
161 friend class SymMatrix;
162 };
163
164 inline std::ostream& operator<<(std::ostream& os,const Matrix& M) {
165 for (unsigned i=0; i<M.nlin(); ++i) {
166 for (unsigned j=0; j<M.ncol(); ++j)
167 os << M(i,j) << ' ';
168 os << std::endl;
169 }
170 return os;
171 }
172
173 inline double Matrix::frobenius_norm() const {
174 const size_t sz = size();
175 if (sz==0)
176 return 0.0;
177
178 #ifdef HAVE_LAPACK
179 double work;
180 const BLAS_INT M = sizet_to_int(nlin());
181 const BLAS_INT N = sizet_to_int(ncol());
182 return DLANGE('F',M,N,data(),M,&work);
183 #else
184 double d = 0.0;
185 for (size_t i=0; i<sz; i++)
186 d += data()[i]*data()[i];
187 return sqrt(d);
188 #endif
189 }
190
191 inline Vector Matrix::operator*(const Vector& v) const {
192 om_assert(ncol()==v.nlin());
193 Vector res(nlin());
194 #ifdef HAVE_BLAS
195 const BLAS_INT M = sizet_to_int(nlin());
196 const BLAS_INT N = sizet_to_int(ncol());
197 DGEMV(CblasNoTrans,M,N,1.0,data(),M,v.data(),1,0.0,res.data(),1);
198 #else
199 res.set(0.0);
200 for (Index j=0; j<ncol(); ++j)
201 for (Index i=0; i<nlin(); ++i)
202 res(i) += (*this)(i,j)*v(j);
203 #endif
204
205 return res;
206 }
207
208 inline Matrix Matrix::submat(const Index istart,const Index isize,const Index jstart,const Index jsize) const {
209 om_assert(istart+isize<=nlin() && jstart+jsize<=ncol());
210
211 Matrix res(isize,jsize);
212
213 for (Index j=0; j<jsize; ++j) {
214 #ifdef HAVE_BLAS
215 const BLAS_INT sz = sizet_to_int(isize);
216 BLAS(dcopy,DCOPY)(sz,data()+istart+(jstart+j)*nlin(),1,res.data()+j*isize,1);
217 #else
218 for (Index i=0; i<isize; ++i)
219 res(i,j) = (*this)(istart+i,jstart+j);
220 #endif
221 }
222 return res;
223 }
224
225 inline void Matrix::insertmat(const Index istart,const Index jstart,const Matrix& B) {
226 om_assert(istart+B.nlin()<=nlin() && jstart+B.ncol()<=ncol() );
227
228 for (Index j=0; j<B.ncol(); ++j)
229 for (Index i=0; i<B.nlin(); ++i)
230 (*this)(istart+i,jstart+j) = B(i,j);
231 }
232
233 inline Vector Matrix::getcol(const Index j) const {
234 om_assert(j<ncol( ));
235 Vector res(nlin());
236 #ifdef HAVE_BLAS
237 const BLAS_INT M = sizet_to_int(nlin());
238 BLAS(dcopy,DCOPY)(M,data()+nlin()*j,1,res.data(),1);
239 #else
240 for (Index i=0; i<nlin(); ++i)
241 res(i) = (*this)(i,j);
242 #endif
243 return res;
244 }
245
246 inline Vector Matrix::getlin(const Index i) const {
247 om_assert(i<nlin());
248 Vector res(ncol());
249 #ifdef HAVE_BLAS
250 const BLAS_INT M = sizet_to_int(nlin());
251 const BLAS_INT N = sizet_to_int(ncol());
252 BLAS(dcopy,DCOPY)(N,data()+i,M,res.data(),1);
253 #else
254 for (Index j=0; j<ncol(); ++j)
255 res(j) = (*this)(i,j);
256 #endif
257 return res;
258 }
259
260 inline void Matrix::setcol(const Index j,const Vector& v) {
261 om_assert(v.size()==nlin() && j<ncol());
262 #ifdef HAVE_BLAS
263 const BLAS_INT M = sizet_to_int(nlin());
264 BLAS(dcopy,DCOPY)(M,v.data(),1,data()+nlin()*j,1);
265 #else
266 for (Index i=0; i<nlin(); ++i)
267 (*this)(i,j) = v(i);
268 #endif
269 }
270
271 inline void Matrix::setlin(const Index i,const Vector& v) {
272 om_assert(v.size()==ncol());
273 om_assert(i<nlin());
274 #ifdef HAVE_BLAS
275 const BLAS_INT M = sizet_to_int(nlin());
276 const BLAS_INT N = sizet_to_int(ncol());
277 BLAS(dcopy,DCOPY)(N,v.data(),1,data()+i,M);
278 #else
279 for (Index j=0; j<ncol(); ++j)
280 (*this)(i,j) = v(j);
281 #endif
282 }
283
284 inline Vector Matrix::tmult(const Vector& v) const {
285 om_assert(nlin()==v.nlin());
286 Vector res(ncol());
287 #ifdef HAVE_BLAS
288 const BLAS_INT M = sizet_to_int(nlin());
289 const BLAS_INT N = sizet_to_int(ncol());
290 DGEMV(CblasTrans,M,N,1.0,data(),sizet_to_int(nlin()),v.data(),1,0.0,res.data(),1);
291 #else
292 for (Index i=0; i<ncol(); ++i) {
293 res(i) = 0;
294 for (Index j=0; j<nlin(); ++j)
295 res(i) += (*this)(j,i)*v(j);
296 }
297 #endif
298
299 return res;
300 }
301
302 inline Matrix Matrix::inverse() const {
303 om_assert(nlin()==ncol());
304 #ifdef HAVE_LAPACK
305 Matrix invA(*this,DEEP_COPY);
306 // LU
307 #if defined(CLAPACK_INTERFACE)
308 const BLAS_INT M = sizet_to_int(invA.nlin());
309 const BLAS_INT N = sizet_to_int(ncol());
310 BLAS_INT* pivots = new BLAS_INT[N];
311 DGETRF(M,N,invA.data(),M,pivots);
312 DGETRI(N,invA.data(),N,pivots);
313 delete[] pivots;
314 #else
315 int Info = 0;
316 BLAS_INT M = sizet_to_int(invA.nlin());
317 BLAS_INT N = sizet_to_int(ncol());
318 int* pivots = new int[N];
319 DGETRF(M,N,invA.data(),M,pivots,Info);
320 const Dimension sz = invA.ncol()*64;
321 double* work=new double[sz];
322 DGETRI(N,invA.data(),N,pivots,work,sz,Info);
323 delete[] pivots;
324 delete[] work;
325 om_assert(Info==0);
326 #endif
327 return invA;
328 #else
329 throw OpenMEEG::maths::LinearAlgebraError("Inverse not implemented, requires LAPACK");
330 #endif
331 }
332
333 inline Matrix Matrix::operator*(const Matrix& B) const {
334 om_assert(ncol()==B.nlin());
335 Matrix C(nlin(),B.ncol());
336 #ifdef HAVE_BLAS
337 const BLAS_INT M = sizet_to_int(nlin());
338 const BLAS_INT N = sizet_to_int(ncol());
339 const BLAS_INT L = sizet_to_int(B.ncol());
340 DGEMM(CblasNoTrans,CblasNoTrans,M,L,N,1.0,data(),M,B.data(),N,0.0,C.data(),M);
341 #else
342 for (Index i=0; i<C.nlin(); ++i)
343 for (Index j=0; j<C.ncol(); ++j) {
344 C(i,j) = 0.0;
345 for (Index k=0; k<ncol(); ++k)
346 C(i,j) += (*this)(i,k)*B(k,j);
347 }
348 #endif
349 return C;
350 }
351
352 inline Matrix Matrix::tmult(const Matrix& B) const {
353 om_assert(nlin()==B.nlin());
354 Matrix C(ncol(),B.ncol());
355 #ifdef HAVE_BLAS
356 const BLAS_INT M = sizet_to_int(nlin());
357 const BLAS_INT N = sizet_to_int(ncol());
358 const BLAS_INT L = sizet_to_int(B.ncol());
359 DGEMM(CblasTrans,CblasNoTrans,N,L,M,1.0,data(),M,B.data(),M,0.0,C.data(),N);
360 #else
361 for (Index i=0; i<C.nlin(); ++i)
362 for (Index j=0; j<C.ncol(); ++j) {
363 C(i,j) = 0.0;
364 for (Index k=0; k<nlin(); ++k)
365 C(i,j) += (*this)(k,i)*B(k,j);
366 }
367 #endif
368 return C;
369 }
370
371 inline Matrix Matrix::multt(const Matrix& B) const {
372 om_assert(ncol()==B.ncol());
373 Matrix C(nlin(),B.nlin());
374 #ifdef HAVE_BLAS
375 const BLAS_INT M = sizet_to_int(nlin());
376 const BLAS_INT N = sizet_to_int(ncol());
377 const BLAS_INT L = sizet_to_int(B.nlin());
378 DGEMM(CblasNoTrans,CblasTrans,M,L,N,1.0,data(),M,B.data(),L,0.0,C.data(),M);
379 #else
380 for (Index j=0; j<C.ncol(); ++j)
381 for (Index i=0; i<C.nlin(); ++i) {
382 C(i,j) = 0.0;
383 for (Index k=0; k<ncol(); ++k)
384 C(i,j) += (*this)(i,k)*B(j,k);
385 }
386 #endif
387 return C;
388 }
389
390 inline Matrix Matrix::tmultt(const Matrix& B) const {
391 om_assert(nlin()==B.ncol());
392 Matrix C(ncol(),B.nlin());
393 #ifdef HAVE_BLAS
394 const BLAS_INT M = sizet_to_int(nlin());
395 const BLAS_INT N = sizet_to_int(ncol());
396 const BLAS_INT L = sizet_to_int(B.nlin());
397 DGEMM(CblasTrans,CblasTrans,L,N,M,1.0,data(),M,B.data(),N,0.0,C.data(),L);
398 #else
399 for (Index i=0; i<C.nlin(); ++i)
400 for (Index j=0; j<C.ncol(); ++j) {
401 C(i,j) = 0.0;
402 for (Index k=0; k<nlin(); ++k)
403 C(i,j) += (*this)(k,i)*B(j,k);
404 }
405 #endif
406 return C;
407 }
408
409 inline Matrix Matrix::operator*(const SymMatrix& B) const {
410 om_assert(ncol()==B.nlin());
411 Matrix C(nlin(),B.ncol());
412
413 // Workaround an MKL bug
414 //#ifdef HAVE_BLAS
415 #if defined(HAVE_BLAS) && !defined(USE_MKL)
416 Matrix D(B);
417 const BLAS_INT m = sizet_to_int(nlin());
418 const BLAS_INT n = sizet_to_int(B.ncol());
419 DSYMM(CblasRight,CblasUpper,m,n,1.0,D.data(),n,data(),m,0.0,C.data(),m);
420 #else
421 for (Index j=0; j<B.ncol(); ++j)
422 for (Index i=0; i<nlin(); ++i) {
423 double sum = 0.0;
424 for (size_t k=0; k<ncol(); ++k)
425 sum += (*this)(i,k)*B(k,j);
426 C(i,j) = sum;
427 }
428 #endif
429 return C;
430 }
431
432 inline void Matrix::operator+=(const Matrix& B) {
433 om_assert(nlin()==B.nlin());
434 om_assert(ncol()==B.ncol());
435 #ifdef HAVE_BLAS
436 const BLAS_INT sz = sizet_to_int(size());
437 BLAS(daxpy,DAXPY)(sz,1.0,B.data(),1,data(),1);
438 #else
439 const size_t sz = size();
440 for (size_t i=0; i<sz; ++i)
441 data()[i] += B.data()[i];
442 #endif
443 }
444
445 inline void Matrix::operator-=(const Matrix& B) {
446 om_assert(nlin()==B.nlin());
447 om_assert(ncol()==B.ncol());
448 #ifdef HAVE_BLAS
449 const BLAS_INT sz = sizet_to_int(size());
450 BLAS(daxpy,DAXPY)(sz,-1.0,B.data(),1,data(),1);
451 #else
452 const size_t sz = size();
453 for (size_t i=0; i<sz; ++i)
454 data()[i] -= B.data()[i];
455 #endif
456 }
457
458 inline double Matrix::dot(const Matrix& B) const {
459 om_assert(nlin()==B.nlin());
460 om_assert(ncol()==B.ncol());
461 #ifdef HAVE_BLAS
462 const BLAS_INT sz = sizet_to_int(size());
463 return BLAS(ddot,DDOT)(sz,data(),1,B.data(),1);
464 #else
465 const sz = size();
466 double s = 0.0;
467 for (size_t i=0; i<sz; i++)
468 s += data()[i]*B.data()[i];
469 return s;
470 #endif
471 }
472}
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:284
Matrix(const char *fname)
Definition matrix.h:40
void operator/=(double x)
double * data() const
Get Matrix data.
Definition matrix.h:66
double & operator()(const Index i, const Index j)
Get Matrix value.
Definition matrix.h:89
void save(const std::string &s) const
Definition matrix.h:153
Vector getcol(const Index j) const
Definition matrix.h:233
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:115
double operator()(const Index i, const Index j) const
Get Matrix value.
Definition matrix.h:81
Vector getlin(const Index i) const
Definition matrix.h:246
Matrix multt(const Matrix &m) const
Definition matrix.h:371
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:208
double frobenius_norm() const
Get Matrix Frobenius norm.
Definition matrix.h:173
void setcol(const Index j, const Vector &v)
Definition matrix.h:260
Matrix operator*(const SparseMatrix &B) const
void setlin(const Index i, const Vector &v)
Definition matrix.h:271
double dot(const Matrix &B) const
Definition matrix.h:458
Matrix transpose() const
void operator-=(const Matrix &B)
Definition matrix.h:445
Matrix operator*(const Matrix &B) const
Definition matrix.h:333
Matrix operator+(const Matrix &B) const
Definition matrix.h:109
void load(const char *filename)
Load Matrix from file (Format set using file name extension)
void operator+=(const Matrix &B)
Definition matrix.h:432
void load(const std::string &s)
Definition matrix.h:154
Matrix(const SparseMatrix &A)
void operator*=(double x)
void insertmat(const Index istart, const Index jstart, const Matrix &B)
Definition matrix.h:225
void reference_data(const double *vals)
Definition matrix.h:51
Matrix pinverse(const double reltol=0.0) const
Matrix inverse() const
Definition matrix.h:302
Matrix tmultt(const Matrix &m) const
Definition matrix.h:390
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
unsigned Index
Definition linop.h:33
BLAS_INT sizet_to_int(const unsigned &num)
Definition linop.h:26
std::ostream & operator<<(std::ostream &os, const Conductivity< REP > &m)
bool empty() const
Definition linop.h:96