Replace Meschach by eigen

.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

@@ -245,6 +245,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. 
-
-
 ----
 
 

share/lib/hoc/import3d/import3d_gui.hoc

@@ -1227,14 +1227,21 @@ proc contour2centroid() {local i, j, imax, imin, ok  localobj mean, pts, d, max,
 	tobj = m.symmeig(m)
 	// major axis is the one with largest eigenvalue
 	major = m.getcol(tobj.max_ind)
+	// For plotting and round-off error consistency
+	// between meschach and eigen use the heuristic that
+	// major axis points in positive direction.
+	// i.e. if greatest value is negative, reverse orientation.
+	if (major.x[major.c.abs().max_ind()] < 0) {
+		major.mul(-1)
+	}
 	// minor is normal and in xy plane
 	minor = m.getcol(3-tobj.min_ind-tobj.max_ind)
 	minor.x[2] = 0
 	minor.div(minor.mag)
-if (g != nil) {
-g.beginline(4, 3) g.line(mean.x[0], mean.x[1])
-g.line(mean.x[0] + 20*major.x[0], mean.x[1] + 20*major.x[1]) g.flush
-}
+	if (g != nil) {
+		g.beginline(4, 3) g.line(mean.x[0], mean.x[1])
+		g.line(mean.x[0] + 20*major.x[0], mean.x[1] + 20*major.x[1]) g.flush
+	}
 	d = new Vector(pts.ncol)
 	rad = new Vector(pts.ncol)
 	for i=0, pts.ncol-1 {

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