63 SymMatrix(
const char* fname):
LinOp(0,0,SYMMETRIC,2),value() { this->load(fname); }
71 size_t size()
const {
return nlin()*(nlin()+1)/2; };
74 size_t ncol()
const {
return nlin(); }
75 size_t&
ncol() {
return nlin(); }
80 bool empty()
const {
return value->empty(); }
82 double*
data()
const {
return value->data; }
84 inline double operator()(
size_t i,
size_t j)
const;
85 inline double& operator()(
size_t i,
size_t j) ;
88 Matrix submat(
size_t istart,
size_t isize,
size_t jstart,
size_t jsize)
const;
90 Vector getlin(
size_t i)
const;
91 void setlin(
size_t i,
const Vector& v);
93 void solveLin(
Vector * B,
int nbvect);
107 void operator *=(
double x);
108 void operator /=(
double x) { (*this)*=(1/x); }
116 void save(
const char *filename)
const;
117 void load(
const char *filename);
119 void save(
const std::string& s)
const {
save(s.c_str()); }
120 void load(
const std::string& s) {
load(s.c_str()); }
128 return data()[i+j*(j+1)/2];
130 return data()[j+i*(i+1)/2];
136 return data()[i+j*(j+1)/2];
138 return data()[j+i*(i+1)/2];
157 std::cout <<
"solveLin not defined" << std::endl;
173 for(
int i = 0; i < nbvect; i++)
179 std::cout <<
"solveLin not defined" << std::endl;
188 for (
size_t i=0;i<
nlin()*(
nlin()+1)/2;i++)
198 for (
size_t i=0;i<
nlin()*(
nlin()+1)/2;i++)
213 std::cerr <<
"Positive definite inverse not defined" << std::endl;
228 std::cout <<
"Big problem in det (DSPTRF)" << std::endl;
229 for (
size_t i = 0; i<
nlin(); i++){
230 if (pivots[i] >= 0) {
233 if (i <
nlin()-1 && pivots[i] == pivots[i+1]) {
234 d *= (invA(i,i)*invA(i+1,i+1)-invA(i,i+1)*invA(i+1,i));
237 std::cout <<
"Big problem in det" << std::endl;
243 std::cerr <<
"Determinant not defined without LAPACK" << std::endl;
281 for (
size_t i=0;i<
nlin()*(
nlin()+1)/2;i++)
294 for (
size_t i=0;i<
nlin()*(
nlin()+1)/2;i++)
308 double *work =
new double[this->
nlin() * 64];
316 std::cerr <<
"!!!!! Inverse not implemented !!!!!" << std::endl;
328 double *work =
new double[this->
nlin() * 64];
336 std::cerr <<
"!!!!! Inverse not implemented !!!!!" << std::endl;
347 for (
size_t i=0;i<
nlin();i++) {
349 for (
size_t j=0;j<
nlin();j++)
350 y(i)+=(*this)(i,j)*v(j);
365 for (
size_t j = 0; j <
ncol(); ++j) this->
operator()(i,j) = v(j);
void save(const std::string &s) const
SymMatrix(size_t M, size_t N)
SymMatrix posdefinverse() const
Matrix operator()(size_t i_start, size_t i_end, size_t j_start, size_t j_end) const
Matrix operator*(const Matrix &B) const
SymMatrix(const char *fname)
utils::RCPtr< LinOpValue > value
Vector solveLin(const Vector &B) const
Matrix solveLin(Matrix &B) const
Vector getlin(size_t i) const
double operator()(size_t i, size_t j) const
void operator+=(const SymMatrix &B)
const SymMatrix & operator=(const double d)
SymMatrix operator*(double x) const
SymMatrix operator+(const SymMatrix &B) const
SymMatrix operator*(const SymMatrix &B) const
void save(const char *filename) const
SymMatrix inverse() const
SymMatrix submat(size_t istart, size_t iend) const
void reference_data(const double *array)
SymMatrix(const SymMatrix &S, const DeepCopy)
void setlin(size_t i, const Vector &v)
void load(const char *filename)
SymMatrix(const Matrix &A)
SymMatrix(const Vector &v)
void load(const std::string &s)
void operator-=(const SymMatrix &B)
SymMatrix operator-(const SymMatrix &B) const
SymMatrix operator/(double x) const
Matrix submat(size_t istart, size_t isize, size_t jstart, size_t jsize) const
OPENMEEGMATHS_EXPORT BLAS_INT sizet_to_int(const size_t &num)
Vect3 operator*(const double &d, const Vect3 &v)