68 explicit Matrix(
const Matrix& A,
const size_t M):
LinOp(A.nlin(),M,FULL,2),value(A.value) { }
73 Matrix(
const char* fname):
LinOp(0,0,FULL,2),value() { this->load(fname); }
89 bool empty()
const {
return value->empty(); }
95 size_t size()
const {
return nlin()*ncol(); };
101 double*
data()
const {
return value->data; }
107 inline double operator()(
size_t i,
size_t j)
const ;
113 inline double& operator()(
size_t i,
size_t j) ;
115 Matrix submat(
size_t istart,
size_t isize,
size_t jstart,
size_t jsize)
const;
116 void insertmat(
size_t istart,
size_t jstart,
const Matrix& B);
117 Vector getcol(
size_t j)
const;
118 void setcol(
size_t j,
const Vector& v);
119 Vector getlin(
size_t i)
const;
120 void setlin(
size_t i,
const Vector& v);
122 const Matrix& set(
const double d);
130 Matrix operator/(
double x)
const;
131 void operator+=(
const Matrix& B);
132 void operator-=(
const Matrix& B);
133 void operator*=(
double x);
134 void operator/=(
double x);
147 Matrix pinverse(
double reltol=0)
const;
154 double frobenius_norm()
const;
155 double dot(
const Matrix& B)
const;
160 void save(
const char *filename)
const;
165 void load(
const char *filename);
167 void save(
const std::string& s)
const {
save(s.c_str()); }
168 void load(
const std::string& s) {
load(s.c_str()); }
207 DGEMV(CblasNoTrans,
sizet_to_int(
nlin()),
sizet_to_int(
ncol()),1.0,
data(),
sizet_to_int(
nlin()),v.
data(),1,0.,y.data(),1);
209 for (
size_t i=0;i<
nlin();i++) {
211 for (
size_t j=0;j<
ncol();j++)
212 y(i) += (*this)(i,j)*v(j);
220 om_assert (istart+isize<=
nlin() && jstart+jsize<=
ncol());
225 for (
size_t j=0; j<jsize; j++) {
227 BLAS(dcopy,DCOPY)(sz,
data()+istart+(jstart+j)*
nlin(),1,a.
data()+j*isize,1);
229 for (
size_t i=0; i<isize; i++)
230 a(i,j) = (*this)(istart+i,jstart+j);
238 for (
size_t j=0; j<B.
ncol(); j++) {
239 for (
size_t i=0; i<B.
nlin(); i++) {
240 (*this)(istart+i,jstart+j)=B(i,j);
289 DGEMV(CblasTrans,
sizet_to_int(
nlin()),
sizet_to_int(
ncol()),1.,
data(),
sizet_to_int(
nlin()),v.
data(),1,0.,y.data(),1);
291 for (
size_t i=0;i<
ncol();i++) {
293 for (
size_t j=0;j<
nlin();j++)
294 y(i)+=(*this)(j,i)*v(j);
306 #if defined(CLAPACK_INTERFACE)
307 #if defined(__APPLE__) && defined(USE_VECLIB) // Apple Veclib Framework (Handles 32 and 64 Bits)
308 __CLPK_integer *pivots =
new __CLPK_integer[
ncol()];
309 __CLPK_integer Info = 0;
310 __CLPK_integer nlin_local = invA.
nlin();
311 __CLPK_integer nlin_local2 = invA.
nlin();
312 __CLPK_integer ncol_local = invA.
ncol();
313 __CLPK_integer sz = invA.
ncol()*64;
314 DGETRF(nlin_local,ncol_local,invA.
data(),nlin_local2,pivots,Info);
315 double *work=
new double[sz];
316 DGETRI(ncol_local,invA.
data(),ncol_local,pivots,work,sz,Info);
330 const unsigned sz = invA.
ncol()*64;
331 double *work=
new double[sz];
339 std::cerr <<
"!!!!! Inverse not implemented !!!!!" << std::endl;
349 DGEMM(CblasNoTrans,CblasNoTrans,
355 for (
size_t i=0;i<C.nlin();i++)
356 for (
size_t j=0;j<C.ncol();j++) {
358 for (
size_t k=0;k<p;k++)
359 C(i,j)+=(*this)(i,k)*B(k,j);
370 DGEMM(CblasTrans,CblasNoTrans,
376 for (
size_t i=0;i<C.nlin();i++)
377 for (
size_t j=0;j<C.ncol();j++) {
379 for (
size_t k=0;k<p;k++)
380 C(i,j)+=(*this)(k,i)*B(k,j);
391 DGEMM(CblasNoTrans,CblasTrans,
397 for (
size_t i=0;i<C.nlin();i++)
398 for (
size_t j=0;j<C.ncol();j++) {
400 for (
size_t k=0;k<p;k++)
401 C(i,j)+=(*this)(i,k)*B(j,k);
412 DGEMM(CblasTrans,CblasTrans,
418 for (
size_t i=0;i<C.nlin();i++)
419 for (
size_t j=0;j<C.ncol();j++) {
421 for (
size_t k=0;k<p;k++)
422 C(i,j)+=(*this)(k,i)*B(j,k);
437 DSYMM(CblasRight,CblasUpper,n,m,1.,D.data(),m,
data(),n,0.,C.data(),l);
439 for (
size_t j=0;j<B.
ncol();j++)
440 for (
size_t i=0;i<
ncol();i++) {
442 for (
size_t k=0;k<
ncol();k++)
443 C(i,j)+=(*this)(i,k)*B(k,j);
456 for (
size_t i=0;i<
nlin()*
ncol();i++)
457 C.data()[i]+=B.
data()[i];
469 for (
size_t i=0;i<
nlin()*
ncol();i++)
470 C.data()[i]-=B.
data()[i];
481 for (
size_t i=0;i<
nlin()*
ncol();i++)
492 for (
size_t i=0;i<
nlin()*
ncol();i++)
503 for (
size_t i=0;i<
nlin()*
ncol();i++)
size_t size() const
Get Matrix size.
Matrix operator*(const Matrix &B) const
Matrix(const char *fname)
utils::RCPtr< LinOpValue > value
Matrix(const Matrix &A, const DeepCopy)
Matrix(const size_t M, const size_t N)
Matrix submat(size_t istart, size_t isize, size_t jstart, size_t jsize) const
void load(const std::string &s)
double frobenius_norm() const
Get Matrix Frobenius norm.
double operator()(size_t i, size_t j) const
Get Matrix value.
Vect3 operator*(const double &d, const Vect3 &v)
double dot(const Matrix &B) const
Matrix multt(const Matrix &m) const
void setcol(size_t j, const Vector &v)
Matrix(const Matrix &A, const size_t M)
virtual size_t ncol() const
Vector getlin(size_t i) const
Matrix operator-(const Matrix &B) const
double * data() const
Get Matrix data.
OPENMEEGMATHS_EXPORT BLAS_INT sizet_to_int(const size_t &num)
Vector getcol(size_t j) const
void reference_data(const double *vals)
void save(const std::string &s) const
Matrix operator+(const Matrix &B) const
Vector tmult(const Vector &v) const
bool empty() const
Test if Matrix is empty.
void operator+=(const Matrix &B)
void insertmat(size_t istart, size_t jstart, const Matrix &B)
void operator-=(const Matrix &B)
Matrix tmultt(const Matrix &m) const
void setlin(size_t i, const Vector &v)