diff --git a/.gitmodules b/.gitmodules index a2538c164c..05fd71cce8 100644 --- a/.gitmodules +++ b/.gitmodules @@ -22,3 +22,6 @@ [submodule "external/CLI11"] path = external/CLI11 url = https://github.com/CLIUtils/CLI11.git +[submodule "external/eigen"] + path = external/eigen + url = https://gitlab.com/libeigen/eigen diff --git a/CMakeLists.txt b/CMakeLists.txt index e479b995aa..725298fdf0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -245,6 +245,8 @@ if(NOT EXISTS "${CODING_CONV_CMAKE}/3rdparty.cmake") endif() include("${CODING_CONV_CMAKE}/3rdparty.cmake") cpp_cc_git_submodule(Random123) +option(EIGEN_BUILD_TESTING "" OFF) +cpp_cc_git_submodule(eigen BUILD PACKAGE Eigen3::Eigen REQUIRED) # ================================================================================================= # Enable sanitizer support if the NRN_SANITIZERS variable is set. Comes befores PythonHelper.cmake. diff --git a/external/eigen b/external/eigen new file mode 160000 index 0000000000..328b5f9085 --- /dev/null +++ b/external/eigen @@ -0,0 +1 @@ +Subproject commit 328b5f90858f93344ebc1484df8cadfd2a5da6dd diff --git a/src/ivoc/matrix.cpp b/src/ivoc/matrix.cpp index ee65d49fdd..824a651e89 100644 --- a/src/ivoc/matrix.cpp +++ b/src/ivoc/matrix.cpp @@ -16,17 +16,6 @@ extern int hoc_return_type_code; extern double hoc_scan(FILE*); extern Object** hoc_temp_objptr(Object*); -#if 0 - extern void install_matrix_method(const char* name, double (*)(...)); - extern void* matrix_arg(int); - extern double* matrix_pelm(void*, int i, int j); - extern int matrix_nrow(void*); - extern int matrix_ncol(void*); - extern int matrix_type(void*); // FULL 1, SPARSE 2, BAND 3 - extern MAT* matrix_full(void*); // hoc_execerror if void* not right type - extern SPMAT* matrix_sparse(void*); -#endif - static void check_domain(int i, int j) { if (i > j || i < 0) { auto const tmp = "index=" + std::to_string(i) + " max_index=" + std::to_string(j) + "\n"; diff --git a/src/ivoc/ocmatrix.cpp b/src/ivoc/ocmatrix.cpp index 123de08771..82f4acb9af 100644 --- a/src/ivoc/ocmatrix.cpp +++ b/src/ivoc/ocmatrix.cpp @@ -4,47 +4,26 @@ #define v_elem(v, i) (*(vector_vec(v) + i)) -#include "ivocvect.h" -#include "oc2iv.h" - -#undef error +#include +#include +#include -extern "C" { -#undef OUT /* /usr/x86_64-w64-mingw32/sys-root/mingw/include/windef.h */ -#include "matrix.h" //meschach -#include "matrix2.h" -#include "sparse.h" -#include "sparse2.h" -extern MAT* m_get(int, int); -} // extern "C" - -int nrn_matrix_dim(void*, int); +#include "ivocvect.h" +#include "oc_ansi.h" #include "ocmatrix.h" -using std::vector; int nrn_matrix_dim(void* vm, int d) { OcMatrix* m = (OcMatrix*) vm; return d ? m->ncol() : m->nrow(); } -static void Vect2VEC(Vect* v1, VEC& v2) { -#ifdef WIN32 - v2.ve = vector_vec(v1); - v2.dim = vector_capacity(v1); - v2.max_dim = vector_buffer_size(v1); -#else - v2.ve = v1->data(); - v2.dim = v1->size(); - v2.max_dim = v1->buffer_size(); -#endif +static Eigen::Map Vect2VEC(Vect* v1) { + return Eigen::Map(v1->data(), v1->size()); } -OcMatrix::OcMatrix(int type) { - obj_ = nullptr; - type_ = type; -} -OcMatrix::~OcMatrix() {} +OcMatrix::OcMatrix(int type) + : type_(type) {} OcMatrix* OcMatrix::instance(int nrow, int ncol, int type) { switch (type) { @@ -60,12 +39,12 @@ void OcMatrix::unimp() { hoc_execerror("Matrix method not implemented for this type matrix", 0); } -void OcMatrix::nonzeros(vector& m, vector& n) { +void OcMatrix::nonzeros(std::vector& m, std::vector& n) { m.clear(); n.clear(); for (int i = 0; i < nrow(); i++) { for (int j = 0; j < ncol(); j++) { - if (getval(i, j)) { + if (getval(i, j) != 0) { m.push_back(i); n.push_back(j); } @@ -73,21 +52,6 @@ void OcMatrix::nonzeros(vector& m, vector& n) { } } -void OcSparseMatrix::nonzeros(vector& m, vector& n) { - m.clear(); - n.clear(); - for (int i = 0; i < m_->m; i++) { - SPROW* const r = m_->row + i; - row_elt* r_elt = r->elt; - for (int k = 0; k < r->len; k++) { - int j = r_elt[k].col; - m.push_back(i); - n.push_back(j); - } - } -} - - OcFullMatrix* OcMatrix::full() { if (type_ != MFULL) { // could clone one maybe hoc_execerror("Matrix is not a FULL matrix (type 1)", 0); @@ -96,243 +60,174 @@ OcFullMatrix* OcMatrix::full() { } OcFullMatrix::OcFullMatrix(int nrow, int ncol) - : OcMatrix(MFULL) { - lu_factor_ = nullptr; - lu_pivot_ = nullptr; - m_ = m_get(nrow, ncol); -} -OcFullMatrix::~OcFullMatrix() { - if (lu_factor_) { - M_FREE(lu_factor_); - PX_FREE(lu_pivot_); - } - M_FREE(m_); + : OcMatrix(MFULL) + , m_(nrow, ncol) { + // The constructor of Eigen::Matrix does not initialize values + m_.setZero(); } + double* OcFullMatrix::mep(int i, int j) { - return &m_->me[i][j]; + return &m_(i, j); } double OcFullMatrix::getval(int i, int j) { - return m_->me[i][j]; + return m_(i, j); } int OcFullMatrix::nrow() { - return m_->m; + return m_.rows(); } int OcFullMatrix::ncol() { - return m_->n; + return m_.cols(); } void OcFullMatrix::resize(int i, int j) { - m_resize(m_, i, j); + m_.conservativeResize(i, j); } void OcFullMatrix::mulv(Vect* vin, Vect* vout) { - VEC v1, v2; - Vect2VEC(vin, v1); - Vect2VEC(vout, v2); - mv_mlt(m_, &v1, &v2); + auto v1 = Vect2VEC(vin); + auto v2 = Vect2VEC(vout); + v2 = m_ * v1; } void OcFullMatrix::mulm(Matrix* in, Matrix* out) { - m_mlt(m_, in->full()->m_, out->full()->m_); + out->full()->m_ = m_ * in->full()->m_; } void OcFullMatrix::muls(double s, Matrix* out) { - sm_mlt(s, m_, out->full()->m_); + out->full()->m_ = s * m_; } void OcFullMatrix::add(Matrix* in, Matrix* out) { - m_add(m_, in->full()->m_, out->full()->m_); + out->full()->m_ = m_ + in->full()->m_; } void OcFullMatrix::copy(Matrix* out) { - m_copy(m_, out->full()->m_); + out->full()->m_ = m_; } void OcFullMatrix::bcopy(Matrix* out, int i0, int j0, int n0, int m0, int i1, int j1) { - m_move(m_, i0, j0, n0, m0, out->full()->m_, i1, j1); + out->full()->m_.block(i1, j1, n0, m0) = m_.block(i0, j0, n0, m0); } void OcFullMatrix::transpose(Matrix* out) { - m_transp(m_, out->full()->m_); + if (out->full()->m_ == m_) { + m_.transposeInPlace(); + } else { + out->full()->m_ = m_.transpose(); + } } +// As only symmetric matrix are accepted, eigenvalues are not complex void OcFullMatrix::symmeigen(Matrix* mout, Vect* vout) { - VEC v1; - Vect2VEC(vout, v1); - symmeig(m_, mout->full()->m_, &v1); + auto v1 = Vect2VEC(vout); + Eigen::EigenSolver es(m_); + v1 = es.eigenvalues().real(); + mout->full()->m_ = es.eigenvectors().real(); } void OcFullMatrix::svd1(Matrix* u, Matrix* v, Vect* d) { - VEC v1; - Vect2VEC(d, v1); - svd(m_, u ? u->full()->m_ : nullptr, v ? v->full()->m_ : nullptr, &v1); + auto v1 = Vect2VEC(d); + Eigen::JacobiSVD svd(m_); + v1 = svd.singularValues(); + if (u) { + u->full()->m_ = svd.matrixU().transpose(); + } + if (v) { + v->full()->m_ = svd.matrixV().transpose(); + } } void OcFullMatrix::getrow(int k, Vect* out) { - VEC v1; - Vect2VEC(out, v1); - get_row(m_, k, &v1); + auto v1 = Vect2VEC(out); + v1 = m_.row(k); } void OcFullMatrix::getcol(int k, Vect* out) { - VEC v1; - Vect2VEC(out, v1); - get_col(m_, k, &v1); + auto v1 = Vect2VEC(out); + v1 = m_.col(k); } void OcFullMatrix::getdiag(int k, Vect* out) { - int i, j, row, col; - row = nrow(); - col = ncol(); + auto vout = m_.diagonal(k); if (k >= 0) { - for (i = 0, j = k; i < row && j < col; ++i, ++j) { -#ifdef WIN32 - v_elem(out, i) = m_entry(m_, i, j); -#else - out->elem(i) = m_entry(m_, i, j); -#endif + for (int i = 0, j = k; i < nrow() && j < ncol(); ++i, ++j) { + out->elem(i) = vout(i); } } else { - // Yes for negative diagonal we set the vector from the middle - // The output vector should ALWAYS be the size of biggest diagonal - for (i = -k, j = 0; i < row && j < col; ++i, ++j) { -#ifdef WIN32 - v_elem(out, i) = m_entry(m_, i, j); -#else - out->elem(i) = m_entry(m_, i, j); -#endif + for (int i = -k, j = 0; i < nrow() && j < ncol(); ++i, ++j) { + out->elem(i) = vout(j); } } } void OcFullMatrix::setrow(int k, Vect* in) { - VEC v1; - Vect2VEC(in, v1); - set_row(m_, k, &v1); + auto v1 = Vect2VEC(in); + m_.row(k) = v1; } void OcFullMatrix::setcol(int k, Vect* in) { - VEC v1; - Vect2VEC(in, v1); - set_col(m_, k, &v1); + auto v1 = Vect2VEC(in); + m_.col(k) = v1; } void OcFullMatrix::setdiag(int k, Vect* in) { - int i, j, row, col; - row = nrow(); - col = ncol(); + auto out = m_.diagonal(k); if (k >= 0) { - for (i = 0, j = k; i < row && j < col; ++i, ++j) { -#ifdef WIN32 - m_set_val(m_, i, j, v_elem(in, i)); -#else - m_set_val(m_, i, j, in->elem(i)); -#endif + for (int i = 0, j = k; i < nrow() && j < ncol(); ++i, ++j) { + out(i) = in->elem(i); } } else { - // Yes for negative diagonal we set the vector from the middle - // The input vector should ALWAYS be the size of biggest diagonal - for (i = -k, j = 0; i < row && j < col; ++i, ++j) { -#ifdef WIN32 - m_set_val(m_, i, j, v_elem(in, i)); -#else - m_set_val(m_, i, j, in->elem(i)); -#endif + for (int i = -k, j = 0; i < nrow() && j < ncol(); ++i, ++j) { + out(j) = in->elem(i); } } + m_.diagonal(k) = out; } void OcFullMatrix::setrow(int k, double in) { - int i, col = ncol(); - for (i = 0; i < col; ++i) { - m_set_val(m_, k, i, in); - } + m_.row(k).fill(in); } void OcFullMatrix::setcol(int k, double in) { - int i, row = nrow(); - for (i = 0; i < row; ++i) { - m_set_val(m_, i, k, in); - } + m_.col(k).fill(in); } void OcFullMatrix::setdiag(int k, double in) { - int i, j, row, col; - row = nrow(); - col = ncol(); - if (k >= 0) { - for (i = 0, j = k; i < row && j < col; ++i, ++j) { - m_set_val(m_, i, j, in); - } - } else { - for (i = -k, j = 0; i < row && j < col; ++i, ++j) { - m_set_val(m_, i, j, in); - } - } + m_.diagonal(k).fill(in); } void OcFullMatrix::zero() { - m_zero(m_); + m_.setZero(); } void OcFullMatrix::ident() { - m_ident(m_); + m_.setIdentity(); } void OcFullMatrix::exp(Matrix* out) { - m_exp(m_, 0., out->full()->m_); + out->full()->m_ = m_.exp(); } void OcFullMatrix::pow(int i, Matrix* out) { - m_pow(m_, i, out->full()->m_); + out->full()->m_ = m_.pow(i).eval(); } void OcFullMatrix::inverse(Matrix* out) { - m_inverse(m_, out->full()->m_); + out->full()->m_ = m_.inverse(); } void OcFullMatrix::solv(Vect* in, Vect* out, bool use_lu) { - bool call_lufac = true; - if (!lu_factor_) { - lu_factor_ = m_get(nrow(), nrow()); - lu_pivot_ = px_get(nrow()); - } else if (use_lu && lu_factor_->m == nrow()) { - call_lufac = false; + if (!lu_ || !use_lu || lu_->rows() != m_.rows()) { + lu_ = std::make_unique>(m_); } - VEC v1, v2; - Vect2VEC(in, v1); - Vect2VEC(out, v2); - if (call_lufac) { - m_resize(lu_factor_, nrow(), nrow()); - m_copy(m_, lu_factor_); - px_resize(lu_pivot_, nrow()); - LUfactor(lu_factor_, lu_pivot_); - } - LUsolve(lu_factor_, lu_pivot_, &v1, &v2); + auto v1 = Vect2VEC(in); + auto v2 = Vect2VEC(out); + v2 = lu_->solve(v1); } double OcFullMatrix::det(int* e) { - int n = nrow(); - MAT* lu = m_get(n, n); - PERM* piv = px_get(n); - m_copy(m_, lu); - LUfactor(lu, piv); - double m = 1.0; *e = 0; - for (int i = 0; i < n; ++i) { - m *= lu->me[i][i]; - if (m == 0.0) { - break; - } - while (std::abs(m) >= 1e12) { - m *= 1e-12; - *e += 12; - } - while (std::abs(m) < 1e-12) { - m *= 1e12; - *e -= 12; - } - } + double m = m_.determinant(); if (m) { while (std::abs(m) >= 10.0) { m *= 0.1; @@ -343,213 +238,144 @@ double OcFullMatrix::det(int* e) { *e -= 1; } } - m *= double(px_sign(piv)); - M_FREE(lu); - PX_FREE(piv); return m; } //-------------------------- OcSparseMatrix::OcSparseMatrix(int nrow, int ncol) - : OcMatrix(MSPARSE) { - /* sp_get -- get sparse matrix - -- len is number of elements available for each row without - allocating further memory */ - - int len = 4; - m_ = sp_get(nrow, ncol, len); - lu_factor_ = nullptr; - lu_pivot_ = nullptr; -} -OcSparseMatrix::~OcSparseMatrix() { - if (lu_factor_) { - SP_FREE(lu_factor_); - PX_FREE(lu_pivot_); - } - SP_FREE(m_); -} - -// returns pointer to sparse element. nullptr if it does not exist. -double* OcSparseMatrix::pelm(int i, int j) { - SPROW* r = m_->row + i; - int idx = sprow_idx(r, j); - if (idx >= 0) { - return &r->elt[idx].val; - } else { - return nullptr; - } -} + : OcMatrix(MSPARSE) + , m_(nrow, ncol) {} double* OcSparseMatrix::mep(int i, int j) { - SPROW* r = m_->row + i; - int idx = sprow_idx(r, j); - if (idx >= 0) { - return &r->elt[idx].val; - } - // does not exist so create it with a value of 0 - sp_set_val(m_, i, j, 0.); - // and try again - idx = sprow_idx(r, j); - return &r->elt[idx].val; + return &m_.coeffRef(i, j); } void OcSparseMatrix::zero() { - sp_zero(m_); + for (int k = 0; k < m_.outerSize(); ++k) { + for (decltype(m_)::InnerIterator it(m_, k); it; ++it) { + it.valueRef() = 0.; + } + } } double OcSparseMatrix::getval(int i, int j) { - return sp_get_val(m_, i, j); + return m_.coeff(i, j); } + int OcSparseMatrix::nrow() { - return m_->m; + return m_.rows(); } + int OcSparseMatrix::ncol() { - return m_->n; + return m_.cols(); } + void OcSparseMatrix::mulv(Vect* vin, Vect* vout) { - VEC v1, v2; - Vect2VEC(vin, v1); - Vect2VEC(vout, v2); - sp_mv_mlt(m_, &v1, &v2); + auto v1 = Vect2VEC(vin); + auto v2 = Vect2VEC(vout); + v2 = m_ * v1; } void OcSparseMatrix::solv(Vect* in, Vect* out, bool use_lu) { - bool call_lufac = true; - if (!lu_factor_) { - lu_factor_ = sp_get(nrow(), nrow(), 4); - lu_pivot_ = px_get(nrow()); - } else if (use_lu && lu_factor_->m == nrow()) { - call_lufac = false; + if (!lu_ || !use_lu || lu_->rows() != m_.rows()) { + m_.makeCompressed(); + lu_ = std::make_unique>(m_); + lu_->analyzePattern(m_); + lu_->factorize(m_); } - VEC v1, v2; - Vect2VEC(in, v1); - Vect2VEC(out, v2); - if (call_lufac) { - sp_resize(lu_factor_, nrow(), nrow()); - sp_copy2(m_, lu_factor_); - px_resize(lu_pivot_, nrow()); - spLUfactor(lu_factor_, lu_pivot_, .9); - } - spLUsolve(lu_factor_, lu_pivot_, &v1, &v2); + auto v1 = Vect2VEC(in); + auto v2 = Vect2VEC(out); + v2 = lu_->solve(v1); } void OcSparseMatrix::setrow(int k, Vect* in) { - VEC v1; - Vect2VEC(in, v1); - int i, n = ncol(); - double* p; - for (i = 0; i < n; ++i) { - if ((p = pelm(k, i)) != nullptr) { -#ifdef WIN32 - *p = v_elem(in, i); - } else if (v_elem(in, i)) { - sp_set_val(m_, k, i, v_elem(in, i)); -#else - *p = in->elem(i); - } else if (in->elem(i)) { - sp_set_val(m_, k, i, in->elem(i)); -#endif - } + int col = m_.cols(); + for (int i = 0; i < col; ++i) { + m_.coeffRef(k, i) = in->elem(i); } } void OcSparseMatrix::setcol(int k, Vect* in) { - VEC v1; - Vect2VEC(in, v1); - int i, n = nrow(); - double* p; - for (i = 0; i < n; ++i) { - if ((p = pelm(i, k)) != nullptr) { -#ifdef WIN32 - *p = v_elem(in, i); - } else if (v_elem(in, i)) { - sp_set_val(m_, i, k, v_elem(in, i)); -#else - *p = in->elem(i); - } else if (in->elem(i)) { - sp_set_val(m_, i, k, in->elem(i)); -#endif - } + int row = m_.rows(); + for (int i = 0; i < row; ++i) { + m_.coeffRef(i, k) = in->elem(i); } } void OcSparseMatrix::setdiag(int k, Vect* in) { - int i, j, row, col; - row = nrow(); - col = ncol(); - double* p; + int row = m_.rows(); + int col = m_.cols(); if (k >= 0) { - for (i = 0, j = k; i < row && j < col; ++i, ++j) { - if ((p = pelm(i, j)) != nullptr) { -#ifdef WIN32 - *p = v_elem(in, i); - } else if (v_elem(in, i)) { - sp_set_val(m_, i, j, v_elem(in, i)); -#else - *p = in->elem(i); - } else if (in->elem(i)) { - sp_set_val(m_, i, j, in->elem(i)); -#endif - } + for (int i = 0, j = k; i < row && j < col; ++i, ++j) { + m_.coeffRef(i, j) = in->elem(i); } } else { - for (i = -k, j = 0; i < row && j < col; ++i, ++j) { - // Yes for negative diagonal we set the vector from the middle - // The input vector should ALWAYS be the size of biggest diagonal - if ((p = pelm(i, j)) != nullptr) { -#ifdef WIN32 - *p = v_elem(in, i); - } else if (v_elem(in, i)) { - sp_set_val(m_, i, j, v_elem(in, i)); -#else - *p = in->elem(i); - } else if (in->elem(i)) { - sp_set_val(m_, i, j, in->elem(i)); -#endif - } + for (int i = -k, j = 0; i < row && j < col; ++i, ++j) { + m_.coeffRef(i, j) = in->elem(i); } } } void OcSparseMatrix::setrow(int k, double in) { - int i, col = ncol(); - for (i = 0; i < col; ++i) { - sp_set_val(m_, k, i, in); + int col = m_.cols(); + for (int i = 0; i < col; ++i) { + m_.coeffRef(k, i) = in; } } void OcSparseMatrix::setcol(int k, double in) { - int i, row = nrow(); - for (i = 0; i < row; ++i) { - sp_set_val(m_, i, k, in); + int row = m_.rows(); + for (int i = 0; i < row; ++i) { + m_.coeffRef(i, k) = in; } } void OcSparseMatrix::ident(void) { - setdiag(0, 1); + m_.setIdentity(); } void OcSparseMatrix::setdiag(int k, double in) { - int i, j, row, col; - row = nrow(); - col = ncol(); + int row = m_.rows(); + int col = m_.cols(); if (k >= 0) { - for (i = 0, j = k; i < row && j < col; ++i, ++j) { - sp_set_val(m_, i, j, in); + for (int i = 0, j = k; i < row && j < col; ++i, ++j) { + m_.coeffRef(i, j) = in; } } else { - for (i = -k, j = 0; i < row && j < col; ++i, ++j) { - sp_set_val(m_, i, j, in); + for (int i = -k, j = 0; i < row && j < col; ++i, ++j) { + m_.coeffRef(i, j) = in; } } } int OcSparseMatrix::sprowlen(int i) { - return m_->row[i].len; + 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) { - *j = m_->row[i].elt[jindx].col; - return m_->row[i].elt[jindx].val; + int acc = 0; + for (decltype(m_)::InnerIterator it(m_, i); it; ++it) { + if (acc == jindx) { + *j = it.col(); + return it.value(); + } + acc += 1; + } + return 0; +} + +void OcSparseMatrix::nonzeros(std::vector& m, std::vector& n) { + m.clear(); + n.clear(); + 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()); + } + } } diff --git a/src/ivoc/ocmatrix.h b/src/ivoc/ocmatrix.h index b1545aa04e..9273ae09ec 100644 --- a/src/ivoc/ocmatrix.h +++ b/src/ivoc/ocmatrix.h @@ -1,26 +1,25 @@ #ifndef ocmatrix_h #define ocmatrix_h -#ifndef MATRIXH -#define MAT void -#define SPMAT void -#define PERM void -#endif - +#include #include -using std::vector; + +#include +#include +#include struct Object; class IvocVect; +class OcMatrix; +using Matrix = OcMatrix; class OcFullMatrix; -#define Vect IvocVect -#define Matrix OcMatrix +using Vect = IvocVect; class OcMatrix { public: enum { MFULL = 1, MSPARSE, MBAND }; static OcMatrix* instance(int nrow, int ncol, int type = MFULL); - virtual ~OcMatrix(); + virtual ~OcMatrix() = default; virtual double* mep(int i, int j) { unimp(); @@ -46,7 +45,7 @@ class OcMatrix { unimp(); } - virtual void nonzeros(vector& m, vector& n); + virtual void nonzeros(std::vector& m, std::vector& n); OcFullMatrix* full(); @@ -145,94 +144,86 @@ class OcMatrix { OcMatrix(int type); public: - Object* obj_; + Object* obj_{}; private: - int type_; + int type_{}; }; extern Matrix* matrix_arg(int); -class OcFullMatrix: public OcMatrix { // type 1 +class OcFullMatrix final: public OcMatrix { // type 1 public: OcFullMatrix(int, int); - virtual ~OcFullMatrix(); - - virtual double* mep(int, int); - virtual double getval(int i, int j); - virtual int nrow(); - virtual int ncol(); - virtual void resize(int, int); - - virtual void mulv(Vect* in, Vect* out); - virtual void mulm(Matrix* in, Matrix* out); - virtual void muls(double, Matrix* out); - virtual void add(Matrix*, Matrix* out); - virtual void getrow(int, Vect* out); - virtual void getcol(int, Vect* out); - virtual void getdiag(int, Vect* out); - virtual void setrow(int, Vect* in); - virtual void setcol(int, Vect* in); - virtual void setdiag(int, Vect* in); - virtual void setrow(int, double in); - virtual void setcol(int, double in); - virtual void setdiag(int, double in); - virtual void zero(); - virtual void ident(); - virtual void exp(Matrix* out); - virtual void pow(int, Matrix* out); - virtual void inverse(Matrix* out); - virtual void solv(Vect* vin, Vect* vout, bool use_lu); - virtual void copy(Matrix* out); - virtual void bcopy(Matrix* mout, int i0, int j0, int n0, int m0, int i1, int j1); - virtual void transpose(Matrix* out); - virtual void symmeigen(Matrix* mout, Vect* vout); - virtual void svd1(Matrix* u, Matrix* v, Vect* d); - virtual double det(int* exponent); + ~OcFullMatrix() override = default; + + double* mep(int, int) override; + double getval(int i, int j) override; + int nrow() override; + int ncol() override; + void resize(int, int) override; + + void mulv(Vect* in, Vect* out) override; + void mulm(Matrix* in, Matrix* out) override; + void muls(double, Matrix* out) override; + void add(Matrix*, Matrix* out) override; + void getrow(int, Vect* out) override; + void getcol(int, Vect* out) override; + void getdiag(int, Vect* out) override; + void setrow(int, Vect* in) override; + void setcol(int, Vect* in) override; + void setdiag(int, Vect* in) override; + void setrow(int, double in) override; + void setcol(int, double in) override; + void setdiag(int, double in) override; + void zero() override; + void ident() override; + void exp(Matrix* out) override; + void pow(int, Matrix* out) override; + void inverse(Matrix* out) override; + void solv(Vect* vin, Vect* vout, bool use_lu) override; + void copy(Matrix* out) override; + void bcopy(Matrix* mout, int i0, int j0, int n0, int m0, int i1, int j1) override; + void transpose(Matrix* out) override; + void symmeigen(Matrix* mout, Vect* vout) override; + void svd1(Matrix* u, Matrix* v, Vect* d) override; + double det(int* exponent) override; private: - MAT* m_; - MAT* lu_factor_; - PERM* lu_pivot_; + Eigen::Matrix m_{}; + std::unique_ptr> lu_{}; }; -class OcSparseMatrix: public OcMatrix { // type 2 +class OcSparseMatrix final: public OcMatrix { // type 2 public: OcSparseMatrix(int, int); - virtual ~OcSparseMatrix(); + ~OcSparseMatrix() override = default; - virtual double* mep(int, int); - virtual double* pelm(int, int); // nullptr if element does not exist - virtual int nrow(); - virtual int ncol(); - virtual double getval(int, int); - virtual void ident(void); - virtual void mulv(Vect* in, Vect* out); - virtual void solv(Vect* vin, Vect* vout, bool use_lu); + double* mep(int, int) override; + int nrow() override; + int ncol() override; + double getval(int, int) override; + void ident(void) override; + void mulv(Vect* in, Vect* out) override; + void solv(Vect* vin, Vect* vout, bool use_lu) override; - virtual void setrow(int, Vect* in); - virtual void setcol(int, Vect* in); - virtual void setdiag(int, Vect* in); - virtual void setrow(int, double in); - virtual void setcol(int, double in); - virtual void setdiag(int, double in); + void setrow(int, Vect* in) override; + void setcol(int, Vect* in) override; + void setdiag(int, Vect* in) override; + void setrow(int, double in) override; + void setcol(int, double in) override; + void setdiag(int, double in) override; - virtual void nonzeros(vector& m, vector& n); + void nonzeros(std::vector& m, std::vector& n) override; - virtual int sprowlen(int); // how many elements in row - virtual double spgetrowval(int i, int jindx, int* j); + int sprowlen(int) override; // how many elements in row + double spgetrowval(int i, int jindx, int* j) override; - virtual void zero(); + void zero() override; private: - SPMAT* m_; - SPMAT* lu_factor_; - PERM* lu_pivot_; + Eigen::SparseMatrix m_{}; + std::unique_ptr> lu_{}; }; -#ifndef MATRIXH -#undef MAT -#undef SPMAT -#endif - #endif diff --git a/src/nrniv/CMakeLists.txt b/src/nrniv/CMakeLists.txt index c296eaedae..2ef6b54a20 100644 --- a/src/nrniv/CMakeLists.txt +++ b/src/nrniv/CMakeLists.txt @@ -411,6 +411,7 @@ include_directories(${NRN_INCLUDE_DIRS}) # ============================================================================= add_library(nrniv_lib ${NRN_LIBRARY_TYPE} ${NRN_NRNIV_LIB_SRC_FILES}) target_link_libraries(nrniv_lib nrngnu) +target_link_libraries(nrniv_lib Eigen3::Eigen) cpp_cc_configure_sanitizers(TARGET nrniv_lib) # Source-directory .cpp needs to find generated .hpp. target_include_directories(nrniv_lib PUBLIC "${NRN_OC_GEN}") diff --git a/test/unit_tests/matrix.cpp b/test/unit_tests/matrix.cpp index fc91b654ca..500e748820 100644 --- a/test/unit_tests/matrix.cpp +++ b/test/unit_tests/matrix.cpp @@ -166,21 +166,21 @@ SCENARIO("A Matrix", "[neuron_ivoc][OcMatrix]") { m.pow(2, &m); REQUIRE(compareMatrix(m, {{42., 72., 72.}, {42., 114., 114.}, {42., 72., 72.}})); } - { // Mescach computing of det has an error - // int e{}; - // double det = m.det(&e); - // REQUIRE(det == 0.); - // REQUIRE(e == 0); + { + int e{}; + double det = m.det(&e); + REQUIRE(det == 0.); + REQUIRE(e == 0); } *m.mep(2, 0) = 1; *m.mep(2, 2) = 2; { - // Mescach computing of det has an error - // int e{}; - // double det = m.det(&e); - // REQUIRE(det == -1.2348_a); - // REQUIRE(e == 5); - } { + int e{}; + double det = m.det(&e); + REQUIRE(det == -1.2348_a); + REQUIRE(e == 5); + } + { OcFullMatrix n(4, 3); m.inverse(&n); n.resize(3, 3); // ??? @@ -202,7 +202,7 @@ SCENARIO("A Matrix", "[neuron_ivoc][OcMatrix]") { { IvocVect v(3); m.getdiag(-2, &v); - REQUIRE(v.vec()[2] == Catch::Detail::Approx({1.})); + REQUIRE(v.vec()[2] == Catch::Detail::Approx(1.)); v.vec() = {1., 0., 0.}; m.setdiag(2, &v); REQUIRE(compareMatrix(m, {{42., 72., 1.}, {72., 114., 114.}, {1., 114., 2.}})); @@ -324,8 +324,6 @@ SCENARIO("A Matrix", "[neuron_ivoc][OcMatrix]") { REQUIRE(*pmep == 1); pmep = m.mep(1, 0); REQUIRE(*pmep == 0); - pmep = m.pelm(0, 1); - REQUIRE(pmep == nullptr); } { int col{};