Skip to content

Commit

Permalink
Constify OcMatrix (#3148)
Browse files Browse the repository at this point in the history 
Minify usage of mep as it get a non stable pointer on the matrix itself.
  • Loading branch information
@alkino
alkino authored Oct 30, 2024
1 parent cd7ac7d commit d52f6ae
Show file tree
Hide file tree
Showing 8 changed files with 162 additions and 153 deletions.
17 changes: 7 additions & 10 deletions src/ivoc/matrix.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,13 +64,10 @@ static double m_ncol(void* v) {

static double m_setval(void* v) {
Matrix* m = (Matrix*) v;
int i, j;
double val, *pval;
i = (int) chkarg(1, 0, m->nrow() - 1);
j = (int) chkarg(2, 0, m->ncol() - 1);
val = *getarg(3);
pval = m->mep(i, j);
*pval = val;
int i = (int) chkarg(1, 0, m->nrow() - 1);
int j = (int) chkarg(2, 0, m->ncol() - 1);
double val = *getarg(3);
m->coeff(i, j) = val;
return val;
}

Expand Down Expand Up @@ -171,7 +168,7 @@ static double m_scanf(void* v) {
m->resize(nrow, ncol);
for (i = 0; i < nrow; ++i)
for (j = 0; j < ncol; ++j) {
*(m->mep(i, j)) = hoc_scan(f);
m->coeff(i, j) = hoc_scan(f);
}
return 0.;
}
Expand Down Expand Up @@ -602,7 +599,7 @@ static Object** m_set(void* v) {
int k;
for (k = 0, i = 0; i < nrow; ++i) {
for (j = 0; j < ncol; ++j) {
*(m->mep(i, j)) = *getarg(++k);
m->coeff(i, j) = *getarg(++k);
}
}
return temp_objvar(m);
Expand Down Expand Up @@ -641,7 +638,7 @@ static Object** m_from_vector(void* v) {
double* ve = vector_vec(vout);
for (j = 0; j < ncol; ++j)
for (i = 0; i < nrow; ++i) {
*(m->mep(i, j)) = ve[k++];
m->coeff(i, j) = ve[k++];
}
return temp_objvar(m);
}
Expand Down
83 changes: 41 additions & 42 deletions src/ivoc/ocmatrix.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,21 +34,20 @@ OcMatrix* OcMatrix::instance(int nrow, int ncol, int type) {
}
}

void OcMatrix::unimp() {
hoc_execerror("Matrix method not implemented for this type matrix", 0);
void OcMatrix::unimp() const {
hoc_execerror("Matrix method not implemented for this type matrix", nullptr);
}

void OcMatrix::nonzeros(std::vector<int>& m, std::vector<int>& n) {
m.clear();
n.clear();
std::vector<std::pair<int, int>> OcMatrix::nonzeros() const {
std::vector<std::pair<int, int>> nzs;
for (int i = 0; i < nrow(); i++) {
for (int j = 0; j < ncol(); j++) {
if (getval(i, j) != 0) {
m.push_back(i);
n.push_back(j);
if (getval(i, j) != 0.) {
nzs.emplace_back(i, j);
}
}
}
return nzs;
}

OcFullMatrix* OcMatrix::full() {
Expand All @@ -65,16 +64,16 @@ OcFullMatrix::OcFullMatrix(int nrow, int ncol)
m_.setZero();
}

double* OcFullMatrix::mep(int i, int j) {
return &m_(i, j);
double& OcFullMatrix::coeff(int i, int j) {
return m_(i, j);
}
double OcFullMatrix::getval(int i, int j) {
double OcFullMatrix::getval(int i, int j) const {
return m_(i, j);
}
int OcFullMatrix::nrow() {
int OcFullMatrix::nrow() const {
return m_.rows();
}
int OcFullMatrix::ncol() {
int OcFullMatrix::ncol() const {
return m_.cols();
}

Expand All @@ -84,29 +83,29 @@ void OcFullMatrix::resize(int i, int j) {
m_.conservativeResizeLike(v);
}

void OcFullMatrix::mulv(Vect* vin, Vect* vout) {
void OcFullMatrix::mulv(Vect* vin, Vect* vout) const {
auto v1 = Vect2VEC(vin);
auto v2 = Vect2VEC(vout);
v2 = m_ * v1;
}

void OcFullMatrix::mulm(Matrix* in, Matrix* out) {
void OcFullMatrix::mulm(Matrix* in, Matrix* out) const {
out->full()->m_ = m_ * in->full()->m_;
}

void OcFullMatrix::muls(double s, Matrix* out) {
void OcFullMatrix::muls(double s, Matrix* out) const {
out->full()->m_ = s * m_;
}

void OcFullMatrix::add(Matrix* in, Matrix* out) {
void OcFullMatrix::add(Matrix* in, Matrix* out) const {
out->full()->m_ = m_ + in->full()->m_;
}

void OcFullMatrix::copy(Matrix* out) {
void OcFullMatrix::copy(Matrix* out) const {
out->full()->m_ = m_;
}

void OcFullMatrix::bcopy(Matrix* out, int i0, int j0, int n0, int m0, int i1, int j1) {
void OcFullMatrix::bcopy(Matrix* out, int i0, int j0, int n0, int m0, int i1, int j1) const {
out->full()->m_.block(i1, j1, n0, m0) = m_.block(i0, j0, n0, m0);
}

Expand All @@ -119,14 +118,14 @@ void OcFullMatrix::transpose(Matrix* out) {
}

// As only symmetric matrix are accepted, eigenvalues are not complex
void OcFullMatrix::symmeigen(Matrix* mout, Vect* vout) {
void OcFullMatrix::symmeigen(Matrix* mout, Vect* vout) const {
auto v1 = Vect2VEC(vout);
Eigen::EigenSolver<Eigen::MatrixXd> es(m_);
v1 = es.eigenvalues().real();
mout->full()->m_ = es.eigenvectors().real();
}

void OcFullMatrix::svd1(Matrix* u, Matrix* v, Vect* d) {
void OcFullMatrix::svd1(Matrix* u, Matrix* v, Vect* d) const {
auto v1 = Vect2VEC(d);
Eigen::JacobiSVD<Eigen::MatrixXd> svd(m_, Eigen::ComputeFullU | Eigen::ComputeFullV);
v1 = svd.singularValues();
Expand All @@ -138,17 +137,17 @@ void OcFullMatrix::svd1(Matrix* u, Matrix* v, Vect* d) {
}
}

void OcFullMatrix::getrow(int k, Vect* out) {
void OcFullMatrix::getrow(int k, Vect* out) const {
auto v1 = Vect2VEC(out);
v1 = m_.row(k);
}

void OcFullMatrix::getcol(int k, Vect* out) {
void OcFullMatrix::getcol(int k, Vect* out) const {
auto v1 = Vect2VEC(out);
v1 = m_.col(k);
}

void OcFullMatrix::getdiag(int k, Vect* out) {
void OcFullMatrix::getdiag(int k, Vect* out) const {
auto vout = m_.diagonal(k);
if (k >= 0) {
for (int i = 0, j = k; i < nrow() && j < ncol(); ++i, ++j) {
Expand Down Expand Up @@ -205,15 +204,15 @@ void OcFullMatrix::ident() {
m_.setIdentity();
}

void OcFullMatrix::exp(Matrix* out) {
void OcFullMatrix::exp(Matrix* out) const {
out->full()->m_ = m_.exp();
}

void OcFullMatrix::pow(int i, Matrix* out) {
void OcFullMatrix::pow(int i, Matrix* out) const {
out->full()->m_ = m_.pow(i).eval();
}

void OcFullMatrix::inverse(Matrix* out) {
void OcFullMatrix::inverse(Matrix* out) const {
out->full()->m_ = m_.inverse();
}

Expand All @@ -226,7 +225,7 @@ void OcFullMatrix::solv(Vect* in, Vect* out, bool use_lu) {
v2 = lu_->solve(v1);
}

double OcFullMatrix::det(int* e) {
double OcFullMatrix::det(int* e) const {
*e = 0;
double m = m_.determinant();
if (m) {
Expand All @@ -248,8 +247,8 @@ OcSparseMatrix::OcSparseMatrix(int nrow, int ncol)
: OcMatrix(MSPARSE)
, m_(nrow, ncol) {}

double* OcSparseMatrix::mep(int i, int j) {
return &m_.coeffRef(i, j);
double& OcSparseMatrix::coeff(int i, int j) {
return m_.coeffRef(i, j);
}

void OcSparseMatrix::zero() {
Expand All @@ -260,19 +259,19 @@ void OcSparseMatrix::zero() {
}
}

double OcSparseMatrix::getval(int i, int j) {
double OcSparseMatrix::getval(int i, int j) const {
return m_.coeff(i, j);
}

int OcSparseMatrix::nrow() {
int OcSparseMatrix::nrow() const {
return m_.rows();
}

int OcSparseMatrix::ncol() {
int OcSparseMatrix::ncol() const {
return m_.cols();
}

void OcSparseMatrix::mulv(Vect* vin, Vect* vout) {
void OcSparseMatrix::mulv(Vect* vin, Vect* vout) const {
auto v1 = Vect2VEC(vin);
auto v2 = Vect2VEC(vout);
v2 = m_ * v1;
Expand Down Expand Up @@ -330,7 +329,7 @@ void OcSparseMatrix::setcol(int k, double in) {
}
}

void OcSparseMatrix::ident(void) {
void OcSparseMatrix::ident() {
m_.setIdentity();
}

Expand All @@ -348,15 +347,15 @@ void OcSparseMatrix::setdiag(int k, double in) {
}
}

int OcSparseMatrix::sprowlen(int i) {
int OcSparseMatrix::sprowlen(int i) const {
int acc = 0;
for (decltype(m_)::InnerIterator it(m_, i); it; ++it) {
acc += 1;
}
return acc;
}

double OcSparseMatrix::spgetrowval(int i, int jindx, int* j) {
double OcSparseMatrix::spgetrowval(int i, int jindx, int* j) const {
int acc = 0;
for (decltype(m_)::InnerIterator it(m_, i); it; ++it) {
if (acc == jindx) {
Expand All @@ -368,13 +367,13 @@ double OcSparseMatrix::spgetrowval(int i, int jindx, int* j) {
return 0;
}

void OcSparseMatrix::nonzeros(std::vector<int>& m, std::vector<int>& n) {
m.clear();
n.clear();
std::vector<std::pair<int, int>> OcSparseMatrix::nonzeros() const {
std::vector<std::pair<int, int>> nzs;
nzs.reserve(m_.nonZeros());
for (int k = 0; k < m_.outerSize(); ++k) {
for (decltype(m_)::InnerIterator it(m_, k); it; ++it) {
m.push_back(it.row());
n.push_back(it.col());
nzs.emplace_back(it.row(), it.col());
}
}
return nzs;
}
Loading

0 comments on commit d52f6ae

Please sign in to comment.