40 Matrix(
const char* fname):
LinOp(0,0,FULL,2),value() { load(fname); }
61 size_t size()
const {
return static_cast<std::streamoff
>(nlin())*
static_cast<std::streamoff
>(ncol()); };
66 double*
data()
const {
return value.get(); }
68 #if defined(SWIGPYTHON) || defined(WIN32)
75 std::shared_ptr<double[]> get_shared_data_ptr() {
return value; }
82 om_assert(i<nlin() && j<ncol());
83 return value[i+nlin()*j];
90 om_assert(i<nlin() && j<ncol());
91 return value[i+nlin()*j];
123 inline void operator+=(
const Matrix& B);
124 inline void operator-=(
const Matrix& B);
142 double frobenius_norm()
const;
143 double dot(
const Matrix& B)
const;
147 void save(
const char* filename)
const;
151 void load(
const char* filename);
153 void save(
const std::string& s)
const {
save(s.c_str()); }
154 void load(
const std::string& s) {
load(s.c_str()); }
165 for (
unsigned i=0; i<M.
nlin(); ++i) {
166 for (
unsigned j=0; j<M.
ncol(); ++j)
174 const size_t sz =
size();
182 return DLANGE(
'F',M,N,
data(),M,&work);
185 for (
size_t i=0; i<sz; i++)
197 DGEMV(CblasNoTrans,M,N,1.0,
data(),M,v.
data(),1,0.0,res.
data(),1);
202 res(i) += (*this)(i,j)*v(j);
209 om_assert(istart+isize<=
nlin() && jstart+jsize<=
ncol());
213 for (
Index j=0; j<jsize; ++j) {
216 BLAS(dcopy,DCOPY)(sz,
data()+istart+(jstart+j)*
nlin(),1,res.
data()+j*isize,1);
218 for (
Index i=0; i<isize; ++i)
219 res(i,j) = (*this)(istart+i,jstart+j);
230 (*
this)(istart+i,jstart+j) = B(i,j);
234 om_assert(j<
ncol( ));
241 res(i) = (*this)(i,j);
252 BLAS(dcopy,DCOPY)(N,
data()+i,M,res.
data(),1);
255 res(j) = (*this)(i,j);
277 BLAS(dcopy,DCOPY)(N,v.
data(),1,
data()+i,M);
290 DGEMV(CblasTrans,M,N,1.0,
data(),
sizet_to_int(
nlin()),v.
data(),1,0.0,res.
data(),1);
295 res(i) += (*this)(j,i)*v(j);
307 #if defined(CLAPACK_INTERFACE)
310 BLAS_INT* pivots =
new BLAS_INT[N];
311 DGETRF(M,N,invA.
data(),M,pivots);
312 DGETRI(N,invA.
data(),N,pivots);
318 int* pivots =
new int[N];
319 DGETRF(M,N,invA.
data(),M,pivots,Info);
321 double* work=
new double[sz];
322 DGETRI(N,invA.
data(),N,pivots,work,sz,Info);
329 throw OpenMEEG::maths::LinearAlgebraError(
"Inverse not implemented, requires LAPACK");
340 DGEMM(CblasNoTrans,CblasNoTrans,M,L,N,1.0,
data(),M,B.
data(),N,0.0,C.
data(),M);
346 C(i,j) += (*this)(i,k)*B(k,j);
359 DGEMM(CblasTrans,CblasNoTrans,N,L,M,1.0,
data(),M,B.
data(),M,0.0,C.
data(),N);
365 C(i,j) += (*this)(k,i)*B(k,j);
378 DGEMM(CblasNoTrans,CblasTrans,M,L,N,1.0,
data(),M,B.
data(),L,0.0,C.
data(),M);
384 C(i,j) += (*this)(i,k)*B(j,k);
397 DGEMM(CblasTrans,CblasTrans,L,N,M,1.0,
data(),M,B.
data(),N,0.0,C.
data(),L);
403 C(i,j) += (*this)(k,i)*B(j,k);
415 #if defined(HAVE_BLAS) && !defined(USE_MKL)
419 DSYMM(CblasRight,CblasUpper,m,n,1.0,D.
data(),n,
data(),m,0.0,C.
data(),m);
424 for (
size_t k=0; k<
ncol(); ++k)
425 sum += (*
this)(i,k)*B(k,j);
437 BLAS(daxpy,DAXPY)(sz,1.0,B.
data(),1,
data(),1);
439 const size_t sz =
size();
440 for (
size_t i=0; i<sz; ++i)
450 BLAS(daxpy,DAXPY)(sz,-1.0,B.
data(),1,
data(),1);
452 const size_t sz =
size();
453 for (
size_t i=0; i<sz; ++i)
463 return BLAS(ddot,DDOT)(sz,
data(),1,B.
data(),1);
467 for (
size_t i=0; i<sz; i++)
virtual Dimension ncol() const
Matrix class Matrix class.
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)
Matrix(const SymMatrix &A)
size_t size() const
Get Matrix size.
Matrix(const Matrix &A, const DeepCopy)
const Matrix & set(const double d)
Matrix operator/(double x) const
Vector tmult(const Vector &v) const
Matrix(const char *fname)
void operator/=(double x)
double * data() const
Get Matrix data.
double & operator()(const Index i, const Index j)
Get Matrix value.
void save(const std::string &s) const
Vector getcol(const Index j) const
Matrix(const Vector &v, const Dimension M, const Dimension N)
Matrix(const Dimension M, const Dimension N)
Matrix operator-(const Matrix &B) const
double operator()(const Index i, const Index j) const
Get Matrix value.
Vector getlin(const Index i) const
Matrix multt(const Matrix &m) const
bool empty() const
Test if Matrix is empty.
Matrix(const Matrix &A, const Dimension M)
void info() const
Print info on Matrix.
Matrix submat(const Index istart, const Index isize, const Index jstart, const Index jsize) const
double frobenius_norm() const
Get Matrix Frobenius norm.
void setcol(const Index j, const Vector &v)
Matrix operator*(const SparseMatrix &B) const
void setlin(const Index i, const Vector &v)
double dot(const Matrix &B) const
void operator-=(const Matrix &B)
Matrix operator*(const Matrix &B) const
Matrix operator+(const Matrix &B) const
void load(const char *filename)
Load Matrix from file (Format set using file name extension)
void operator+=(const Matrix &B)
void load(const std::string &s)
Matrix(const SparseMatrix &A)
void operator*=(double x)
void insertmat(const Index istart, const Index jstart, const Matrix &B)
void reference_data(const double *vals)
Matrix pinverse(const double reltol=0.0) const
Matrix tmultt(const Matrix &m) const
Matrix operator*(double x) const
Matrix(const std::string &fname)
Vect3 operator*(const double d, const Vect3 &V)
BLAS_INT sizet_to_int(const unsigned &num)
std::ostream & operator<<(std::ostream &os, const Conductivity< REP > &m)