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(); }
72 om_assert(i<nlin() && j<ncol());
73 return value[i+nlin()*j];
80 om_assert(i<nlin() && j<ncol());
81 return value[i+nlin()*j];
113 inline void operator+=(
const Matrix& B);
114 inline void operator-=(
const Matrix& B);
132 double frobenius_norm()
const;
133 double dot(
const Matrix& B)
const;
137 void save(
const char* filename)
const;
141 void load(
const char* filename);
143 void save(
const std::string& s)
const {
save(s.c_str()); }
144 void load(
const std::string& s) {
load(s.c_str()); }
155 for (
unsigned i=0; i<M.
nlin(); ++i) {
156 for (
unsigned j=0; j<M.
ncol(); ++j)
164 const size_t sz =
size();
172 return DLANGE(
'F',M,N,
data(),M,&work);
175 for (
size_t i=0; i<sz; i++)
187 DGEMV(CblasNoTrans,M,N,1.0,
data(),M,v.
data(),1,0.0,res.
data(),1);
192 res(i) += (*this)(i,j)*v(j);
199 om_assert(istart+isize<=
nlin() && jstart+jsize<=
ncol());
203 for (
Index j=0; j<jsize; ++j) {
206 BLAS(dcopy,DCOPY)(sz,
data()+istart+(jstart+j)*
nlin(),1,res.
data()+j*isize,1);
208 for (
Index i=0; i<isize; ++i)
209 res(i,j) = (*this)(istart+i,jstart+j);
220 (*
this)(istart+i,jstart+j) = B(i,j);
224 om_assert(j<
ncol( ));
231 res(i) = (*this)(i,j);
242 BLAS(dcopy,DCOPY)(N,
data()+i,M,res.
data(),1);
245 res(j) = (*this)(i,j);
267 BLAS(dcopy,DCOPY)(N,v.
data(),1,
data()+i,M);
280 DGEMV(CblasTrans,M,N,1.0,
data(),
sizet_to_int(
nlin()),v.
data(),1,0.0,res.
data(),1);
285 res(i) += (*this)(j,i)*v(j);
297 #if defined(CLAPACK_INTERFACE)
300 BLAS_INT* pivots =
new BLAS_INT[N];
301 DGETRF(M,N,invA.
data(),M,pivots);
302 DGETRI(N,invA.
data(),N,pivots);
308 int* pivots =
new int[N];
309 DGETRF(M,N,invA.
data(),M,pivots,Info);
311 double* work=
new double[sz];
312 DGETRI(N,invA.
data(),N,pivots,work,sz,Info);
319 throw OpenMEEG::maths::LinearAlgebraError(
"Inverse not implemented, requires LAPACK");
330 DGEMM(CblasNoTrans,CblasNoTrans,M,L,N,1.0,
data(),M,B.
data(),N,0.0,C.
data(),M);
336 C(i,j) += (*this)(i,k)*B(k,j);
349 DGEMM(CblasTrans,CblasNoTrans,N,L,M,1.0,
data(),M,B.
data(),M,0.0,C.
data(),N);
355 C(i,j) += (*this)(k,i)*B(k,j);
368 DGEMM(CblasNoTrans,CblasTrans,M,L,N,1.0,
data(),M,B.
data(),L,0.0,C.
data(),M);
374 C(i,j) += (*this)(i,k)*B(j,k);
387 DGEMM(CblasTrans,CblasTrans,L,N,M,1.0,
data(),M,B.
data(),N,0.0,C.
data(),L);
393 C(i,j) += (*this)(k,i)*B(j,k);
405 #if defined(HAVE_BLAS) && !defined(USE_MKL)
409 DSYMM(CblasRight,CblasUpper,m,n,1.0,D.
data(),n,
data(),m,0.0,C.
data(),m);
414 for (
size_t k=0; k<
ncol(); ++k)
415 sum += (*
this)(i,k)*B(k,j);
427 BLAS(daxpy,DAXPY)(sz,1.0,B.
data(),1,
data(),1);
429 const size_t sz =
size();
430 for (
size_t i=0; i<sz; ++i)
440 BLAS(daxpy,DAXPY)(sz,-1.0,B.
data(),1,
data(),1);
442 const size_t sz =
size();
443 for (
size_t i=0; i<sz; ++i)
453 return BLAS(ddot,DDOT)(sz,
data(),1,B.
data(),1);
457 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)