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
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
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
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: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
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