Replace Meschach by Eigen (#2470)

* As part of the legacy code modernization, we are replacing old in-source
copy of ~20-year-old Meschach library with a modern alternative - Eigen!

* However, be aware that this change will inevitably impact results, particularly
if we expect the floating point/spike times identical results for simulations.

* In #2491, we did a detailed analysis of the root cause. The conclusion is that the
floating point differences can be safely ignored.

Co-authored-by: Michael Hines <[email protected]>

.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.git

.sanitizers/undefined.supp

@@ -1,6 +1,9 @@
 function:constructor(Object*)
 function:ctng.cpp
 implicit-integer-sign-change:double vector[2] Eigen::internal::pabs<double vector[2]>(double vector[2] const&)
+implicit-integer-sign-change:double vector[2] Eigen::internal::pnegate<double vector[2]>(double vector[2] const&)
+implicit-integer-sign-change:Eigen::internal::Packet1cd Eigen::internal::pconj<Eigen::internal::Packet1cd>(Eigen::internal::Packet1cd const&)
+implicit-integer-sign-change:Eigen::internal::Packet1cd Eigen::internal::pmul<Eigen::internal::Packet1cd>(Eigen::internal::Packet1cd const&, Eigen::internal::Packet1cd const&)
 implicit-integer-sign-change:fgets_unlimited_nltrans(HocStr*, _IO_FILE*, int)
 implicit-integer-sign-change:invlfiresha.c
 implicit-integer-sign-change:ivRegexp::Search(char const*, int, int, int)

CMakeLists.txt

@@ -243,6 +243,7 @@ if(NOT EXISTS "${CODING_CONV_CMAKE}/3rdparty.cmake")
 endif()
 include("${CODING_CONV_CMAKE}/3rdparty.cmake")
 cpp_cc_git_submodule(Random123)
+cpp_cc_git_submodule(eigen)
 
 # =================================================================================================
 # Enable sanitizer support if the NRN_SANITIZERS variable is set. Comes befores PythonHelper.cmake.

cmake/NeuronFileLists.cmake

@@ -316,76 +316,6 @@ nrn_create_file_list(
   sundialsmath.c)
 set(NRN_SUNDIALS_SRC_FILES ${SUNDIALS_CVODES} ${SUNDIALS_IDA} ${SUNDIALS_SHARED})
 
-# meschach matrix sources
-set(MESCH_FILES_LIST
-    arnoldi.c
-    bdfactor.c
-    bkpfacto.c
-    chfactor.c
-    arnoldi.c
-    bdfactor.c
-    bkpfacto.c
-    chfactor.c
-    conjgrad.c
-    copy.c
-    dmacheps.c
-    err.c
-    extras.c
-    fft.c
-    givens.c
-    hessen.c
-    hsehldr.c
-    init.c
-    iter0.c
-    iternsym.c
-    itersym.c
-    ivecop.c
-    lanczos.c
-    lufactor.c
-    machine.c
-    matlab.c
-    matop.c
-    matrixio.c
-    meminfo.c
-    memory.c
-    memstat.c
-    mfunc.c
-    norm.c
-    otherio.c
-    pxop.c
-    qrfactor.c
-    schur.c
-    solve.c
-    sparse.c
-    sparseio.c
-    spbkp.c
-    spchfctr.c
-    splufctr.c
-    sprow.c
-    spswap.c
-    submat.c
-    svd.c
-    symmeig.c
-    update.c
-    vecop.c
-    version.c
-    zcopy.c
-    zfunc.c
-    zgivens.c
-    zhessen.c
-    zhsehldr.c
-    zlufctr.c
-    zmachine.c
-    zmatio.c
-    zmatlab.c
-    zmatop.c
-    zmemory.c
-    znorm.c
-    zqrfctr.c
-    zschur.c
-    zsolve.c
-    zvecop.c)
-
 set(SPARSE_FILES_LIST bksub.cpp getelm.cpp lineq.cpp prmat.cpp subrows.cpp)
 
 # scopmath sources
@@ -530,7 +460,6 @@ nrn_create_file_list(NRN_PARALLEL_SRC_FILES ${PROJECT_SOURCE_DIR}/src/nrniv
                      nvector_nrnparallel_ld.cpp)
 nrn_create_file_list(NRN_PARALLEL_SRC_FILES ${PROJECT_SOURCE_DIR}/src/sundials/shared
                      nvector_parallel.c)
-nrn_create_file_list(NRN_MESCH_SRC_FILES ${PROJECT_SOURCE_DIR}/src/mesch ${MESCH_FILES_LIST})
 nrn_create_file_list(NRN_SPARSE_SRC_FILES ${PROJECT_SOURCE_DIR}/src/sparse ${SPARSE_FILES_LIST})
 nrn_create_file_list(NRN_SCOPMATH_SRC_FILES ${PROJECT_SOURCE_DIR}/src/scopmath
                      ${SCOPMATH_FILES_LIST})

docs/python/programming/math/matrix.rst

@@ -500,10 +500,6 @@ Matrix
 
 
 
-    .. warning::
-        Implemented only for full and sparse matrices. 
-
-
 ----
 
 
@@ -660,9 +656,6 @@ Matrix
             print()
             m.solv(b,1).printf("%8.3f", 475, 535) 
 
-    .. warning::
-        Implemented only for full and sparse matrices. 
-
 
 ----
 
@@ -790,9 +783,6 @@ Matrix
             m.printf()
 
 
-    .. warning::
-        Implemented only for full and sparse matrices. 
-
 
 ----
 
@@ -813,9 +803,6 @@ Matrix
 
         Otherwise fill the matrix row with a constant. 
 
-    .. warning::
-        Implemented only for full matrices and sparse. 
-
 
 ----
 
@@ -836,9 +823,6 @@ Matrix
 
         Otherwise fill the matrix column with a constant. 
 
-    .. warning::
-        Implemented only for full matrices. 
-
 
 ----
 
@@ -883,9 +867,6 @@ Matrix
             m.printf()
 
 
-    .. warning::
-        Implemented only for full and sparse matrices. 
-
 
 ----
 
@@ -901,9 +882,6 @@ Matrix
     Description:
         Fills the matrix with 0. 
 
-    .. warning::
-        Implemented only for full matrices. 
-
 
 ----
 
@@ -929,10 +907,6 @@ Matrix
             m.printf() 
 
 
-    .. warning::
-        Implemented only for full matrices. 
-
-
 ----
 
 

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";