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) ;
87 Matrix operator()(
size_t i_start,
size_t i_end,
size_t j_start,
size_t j_end)
const;
88 Matrix submat(
size_t istart,
size_t isize,
size_t jstart,
size_t jsize)
const;
89 SymMatrix submat(
size_t istart,
size_t iend)
const;
90 Vector getlin(
size_t i)
const;
91 void setlin(
size_t i,
const Vector& v);
93 void solveLin(
Vector * B,
int nbvect);
96 const SymMatrix& operator=(
const double d);
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++)
282 C.data()[i]+=B.
data()[i];
294 for (
size_t i=0;i<
nlin()*(
nlin()+1)/2;i++)
295 C.data()[i]-=B.
data()[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);
SymMatrix(const char *fname)
SymMatrix(size_t M, size_t N)
void operator+=(const SymMatrix &B)
SymMatrix(const SymMatrix &S, const DeepCopy)
utils::RCPtr< LinOpValue > value
SymMatrix operator*(const SymMatrix &B) const
Vector getlin(size_t i) const
SymMatrix operator/(double x) const
Vect3 operator*(const double &d, const Vect3 &v)
double operator()(size_t i, size_t j) const
SymMatrix operator-(const SymMatrix &B) const
void operator-=(const SymMatrix &B)
void setlin(size_t i, const Vector &v)
OPENMEEGMATHS_EXPORT BLAS_INT sizet_to_int(const size_t &num)
Vector solveLin(const Vector &B) const
void save(const std::string &s) const
void load(const std::string &s)
SymMatrix posdefinverse() const
void reference_data(const double *array)
SymMatrix operator+(const SymMatrix &B) const
SymMatrix inverse() const