From 31c9245eee45a92f39d0937985756523fc578ef4 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 9 Oct 2024 16:08:21 +0200 Subject: [PATCH 01/17] Use nb::iterator instead of Py_Iterator (#3120) --- src/nrnpython/nrnpy_p2h.cpp | 69 +++++++++---------------------------- src/utils/enumerate.h | 4 +-- 2 files changed, 18 insertions(+), 55 deletions(-) diff --git a/src/nrnpython/nrnpy_p2h.cpp b/src/nrnpython/nrnpy_p2h.cpp index f84cb29bcf..8ea32f97b2 100644 --- a/src/nrnpython/nrnpy_p2h.cpp +++ b/src/nrnpython/nrnpy_p2h.cpp @@ -14,6 +14,7 @@ #include "parse.hpp" #include + namespace nb = nanobind; static char* nrnpyerr_str(); @@ -926,8 +927,8 @@ static Object* py_alltoall_type(int size, int type) { size = 0; // calculate dest size (cannot be -1 so cannot return it) } - char* s = NULL; - int* scnt = NULL; + std::vector s{}; + std::vector scnt{}; int* sdispl = NULL; char* r = NULL; int* rcnt = NULL; @@ -937,49 +938,18 @@ static Object* py_alltoall_type(int size, int type) { // for alltoall, each rank handled identically // for scatter, root handled as list all, other ranks handled as None if (type == 1 || nrnmpi_myid == root) { // psrc is list of nhost items - - scnt = new int[np]; - for (int i = 0; i < np; ++i) { - scnt[i] = 0; - } - - PyObject* iterator = PyObject_GetIter(psrc); - PyObject* p; - - size_t bufsz = 100000; // 100k buffer to start with - if (size > 0) { // or else the positive number specified - bufsz = size; - } - if (size >= 0) { // otherwise count only - s = new char[bufsz]; - } - size_t curpos = 0; - for (size_t i = 0; (p = PyIter_Next(iterator)) != NULL; ++i) { - if (p == Py_None) { - scnt[i] = 0; - Py_DECREF(p); + scnt.reserve(np); + for (const nb::handle& p: nb::list(psrc)) { + if (p.is_none()) { + scnt.push_back(0); continue; } - auto b = pickle(p); + const std::vector b = pickle(p.ptr()); if (size >= 0) { - if (curpos + b.size() >= bufsz) { - bufsz = bufsz * 2 + b.size(); - char* s2 = new char[bufsz]; - for (size_t i = 0; i < curpos; ++i) { - s2[i] = s[i]; - } - delete[] s; - s = s2; - } - for (size_t j = 0; j < b.size(); ++j) { - s[curpos + j] = b[j]; - } + s.insert(std::end(s), std::begin(b), std::end(b)); } - curpos += b.size(); - scnt[i] = static_cast(b.size()); - Py_DECREF(p); + scnt.push_back(static_cast(b.size())); } - Py_DECREF(iterator); // scatter equivalent to alltoall NONE list for not root ranks. } else if (type == 5 && nrnmpi_myid != root) { @@ -996,26 +966,23 @@ static Object* py_alltoall_type(int size, int type) { } sdispl = mk_displ(ones); rcnt = new int[np]; - nrnmpi_int_alltoallv(scnt, ones, sdispl, rcnt, ones, sdispl); + nrnmpi_int_alltoallv(scnt.data(), ones, sdispl, rcnt, ones, sdispl); delete[] ones; delete[] sdispl; // exchange - sdispl = mk_displ(scnt); + sdispl = mk_displ(scnt.data()); rdispl = mk_displ(rcnt); if (size < 0) { pdest = PyTuple_New(2); PyTuple_SetItem(pdest, 0, Py_BuildValue("l", (long) sdispl[np])); PyTuple_SetItem(pdest, 1, Py_BuildValue("l", (long) rdispl[np])); - delete[] scnt; delete[] sdispl; delete[] rcnt; delete[] rdispl; } else { char* r = new char[rdispl[np] + 1]; // force > 0 for all None case - nrnmpi_char_alltoallv(s, scnt, sdispl, r, rcnt, rdispl); - delete[] s; - delete[] scnt; + nrnmpi_char_alltoallv(s.data(), scnt.data(), sdispl, r, rcnt, rdispl); delete[] sdispl; pdest = char2pylist(r, np, rcnt, rdispl); @@ -1029,18 +996,14 @@ static Object* py_alltoall_type(int size, int type) { // destination counts rcnt = new int[1]; - nrnmpi_int_scatter(scnt, rcnt, 1, root); + nrnmpi_int_scatter(scnt.data(), rcnt, 1, root); std::vector r(rcnt[0] + 1); // rcnt[0] can be 0 // exchange if (nrnmpi_myid == root) { - sdispl = mk_displ(scnt); + sdispl = mk_displ(scnt.data()); } - nrnmpi_char_scatterv(s, scnt, sdispl, r.data(), rcnt[0], root); - if (s) - delete[] s; - if (scnt) - delete[] scnt; + nrnmpi_char_scatterv(s.data(), scnt.data(), sdispl, r.data(), rcnt[0], root); if (sdispl) delete[] sdispl; diff --git a/src/utils/enumerate.h b/src/utils/enumerate.h index 1a0b38f2ad..e9492887bc 100644 --- a/src/utils/enumerate.h +++ b/src/utils/enumerate.h @@ -99,7 +99,7 @@ constexpr auto enumerate(T&& iterable) { ++iter; } auto operator*() const { - return std::tie(i, *iter); + return std::forward_as_tuple(i, *iter); } }; struct iterable_wrapper { @@ -129,7 +129,7 @@ constexpr auto renumerate(T&& iterable) { ++iter; } auto operator*() const { - return std::tie(i, *iter); + return std::forward_as_tuple(i, *iter); } }; struct iterable_wrapper { From c9ee32e3f082f2ac84c9e80e83a8a92d1f52ea1d Mon Sep 17 00:00:00 2001 From: JCGoran Date: Wed, 9 Oct 2024 16:39:08 +0200 Subject: [PATCH 02/17] Add NRN_PARALLEL_BUILDS to setup.py (#3122) The environment variable `NRN_PARALLEL_BUILDS` can be used to control the build parallelism. Essentially, `-j${NRN_PARALLEL_BUILDS}`. --- docs/install/python_wheels.md | 2 ++ setup.py | 7 ++++++- 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/docs/install/python_wheels.md b/docs/install/python_wheels.md index 42387bfad8..8d0b1755c9 100644 --- a/docs/install/python_wheels.md +++ b/docs/install/python_wheels.md @@ -171,6 +171,8 @@ bash packaging/python/build_wheels.bash linux 3* coreneuron ``` Where we are passing `3*` to build the wheels with `CoreNEURON` support for all python 3 versions. +You can also control the level of parallelization used for the build using the `NRN_PARALLEL_BUILDS` env variable (default: 4). + ### macOS As mentioned above, for macOS all dependencies have to be available on a system. You have to then clone NEURON repository and execute: diff --git a/setup.py b/setup.py index e15b6e5ac0..ba0922b273 100644 --- a/setup.py +++ b/setup.py @@ -247,7 +247,12 @@ def _run_cmake(self, ext): if self.cmake_defs: cmake_args += ["-D" + opt for opt in self.cmake_defs.split(",")] - build_args = ["--config", cfg, "--", "-j4"] # , 'VERBOSE=1'] + build_args = [ + "--config", + cfg, + "--", + f"-j{os.environ.get('NRN_PARALLEL_BUILDS', 4)}", + ] env = os.environ.copy() env["CXXFLAGS"] = "{} -DVERSION_INFO='{}'".format( From aecfda00df71b504a78f8640b580dc79e867e469 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Mon, 14 Oct 2024 16:12:35 +0200 Subject: [PATCH 03/17] cleaning around hoccommand_exec_help1 (#3124) --- src/nrnpython/nrnpy_p2h.cpp | 108 +++++++++++++----------------------- 1 file changed, 38 insertions(+), 70 deletions(-) diff --git a/src/nrnpython/nrnpy_p2h.cpp b/src/nrnpython/nrnpy_p2h.cpp index 8ea32f97b2..fc0a4d085b 100644 --- a/src/nrnpython/nrnpy_p2h.cpp +++ b/src/nrnpython/nrnpy_p2h.cpp @@ -18,7 +18,7 @@ namespace nb = nanobind; static char* nrnpyerr_str(); -static PyObject* nrnpy_pyCallObject(PyObject*, PyObject*); +static nb::object nrnpy_pyCallObject(nb::callable, PyObject*); static PyObject* main_module; static PyObject* main_namespace; @@ -42,15 +42,12 @@ static void p_destruct(void* v) { Member_func p_members[] = {{nullptr, nullptr}}; static void call_python_with_section(Object* pyact, Section* sec) { - PyObject* po = ((Py2Nrn*) pyact->u.this_pointer)->po_; - PyObject* r; + nb::callable po = nb::borrow(((Py2Nrn*) pyact->u.this_pointer)->po_); nanobind::gil_scoped_acquire lock{}; - PyObject* args = PyTuple_Pack(1, (PyObject*) newpysechelp(sec)); - r = nrnpy_pyCallObject(po, args); - Py_XDECREF(args); - Py_XDECREF(r); - if (!r) { + nb::tuple args = nb::make_tuple(reinterpret_cast(newpysechelp(sec))); + nb::object r = nrnpy_pyCallObject(po, args.ptr()); + if (!r.is_valid()) { char* mes = nrnpyerr_str(); if (mes) { Fprintf(stderr, "%s\n", mes); @@ -109,12 +106,12 @@ Object* nrnpy_pyobject_in_obj(PyObject* po) { return on; } -static PyObject* nrnpy_pyCallObject(PyObject* callable, PyObject* args) { +static nb::object nrnpy_pyCallObject(nb::callable callable, PyObject* args) { // When hoc calls a PythonObject method, then in case python // calls something back in hoc, the hoc interpreter must be // at the top level HocTopContextSet - PyObject* p = PyObject_CallObject(callable, args); + PyObject* p = PyObject_CallObject(callable.ptr(), args); #if 0 printf("PyObject_CallObject callable\n"); PyObject_Print(callable, stdout, 0); @@ -142,7 +139,7 @@ printf("\nreturn %p\n", p); } } **/ - return p; + return nb::steal(p); } static void py2n_component(Object* ob, Symbol* sym, int nindex, int isfunc) { @@ -199,7 +196,7 @@ static void py2n_component(Object* ob, Symbol* sym, int nindex, int isfunc) { } } // printf("PyObject_CallObject %s %p\n", sym->name, tail); - result = nrnpy_pyCallObject(tail, args); + result = nrnpy_pyCallObject(nb::borrow(tail), args).release().ptr(); Py_DECREF(args); // PyObject_Print(result, stdout, 0); // printf(" result of call\n"); @@ -330,36 +327,22 @@ static void hpoasgn(Object* o, int type) { } } -static PyObject* hoccommand_exec_help1(PyObject* po) { - PyObject* r; - // PyObject_Print(po, stdout, 0); - // printf("\n"); - if (PyTuple_Check(po)) { - PyObject* args = PyTuple_GetItem(po, 1); - if (!PyTuple_Check(args)) { - args = PyTuple_Pack(1, args); - } else { - Py_INCREF(args); +static nb::object hoccommand_exec_help1(nb::object po) { + if (nb::tuple::check_(po)) { + nb::object args = po[1]; + if (!nb::tuple::check_(args)) { + args = nb::make_tuple(args); } - // PyObject_Print(PyTuple_GetItem(po, 0), stdout, 0); - // printf("\n"); - // PyObject_Print(args, stdout, 0); - // printf("\n"); - // printf("threadstate %p\n", PyThreadState_GET()); - r = nrnpy_pyCallObject(PyTuple_GetItem(po, 0), args); - Py_DECREF(args); + return nrnpy_pyCallObject(po[0], args.ptr()); } else { - PyObject* args = PyTuple_New(0); - r = nrnpy_pyCallObject(po, args); - Py_DECREF(args); + return nrnpy_pyCallObject(nb::borrow(po), nb::tuple().ptr()); } - return r; } -static PyObject* hoccommand_exec_help(Object* ho) { +static nb::object hoccommand_exec_help(Object* ho) { PyObject* po = ((Py2Nrn*) ho->u.this_pointer)->po_; // printf("%s\n", hoc_object_name(ho)); - return hoccommand_exec_help1(po); + return hoccommand_exec_help1(nb::borrow(po)); } static double praxis_efun(Object* ho, Object* v) { @@ -370,9 +353,9 @@ static double praxis_efun(Object* ho, Object* v) { PyObject* po = Py_BuildValue("(OO)", pc, pv); Py_XDECREF(pc); Py_XDECREF(pv); - PyObject* r = hoccommand_exec_help1(po); + nb::object r = hoccommand_exec_help1(nb::borrow(po)); Py_XDECREF(po); - if (!r) { + if (!r.is_valid()) { char* mes = nrnpyerr_str(); if (mes) { Fprintf(stderr, "%s\n", mes); @@ -384,18 +367,14 @@ static double praxis_efun(Object* ho, Object* v) { } return 1e9; // SystemExit? } - PyObject* pn = PyNumber_Float(r); - double x = PyFloat_AsDouble(pn); - Py_XDECREF(pn); - Py_XDECREF(r); - return x; + return static_cast(nb::float_(r)); } static int hoccommand_exec(Object* ho) { nanobind::gil_scoped_acquire lock{}; - PyObject* r = hoccommand_exec_help(ho); - if (r == NULL) { + nb::object r = hoccommand_exec_help(ho); + if (!r.is_valid()) { char* mes = nrnpyerr_str(); if (mes) { std::string tmp{"Python Callback failed [hoccommand_exec]:\n"}; @@ -407,21 +386,18 @@ static int hoccommand_exec(Object* ho) { PyErr_Print(); } } - Py_XDECREF(r); - return (r != NULL); + return r.is_valid(); } static int hoccommand_exec_strret(Object* ho, char* buf, int size) { nanobind::gil_scoped_acquire lock{}; - PyObject* r = hoccommand_exec_help(ho); - if (r) { - PyObject* pn = PyObject_Str(r); - Py2NRNString str(pn); - Py_XDECREF(pn); + nb::object r = hoccommand_exec_help(ho); + if (r.is_valid()) { + nb::str pn(r); + Py2NRNString str(pn.ptr()); strncpy(buf, str.c_str(), size); buf[size - 1] = '\0'; - Py_XDECREF(r); } else { char* mes = nrnpyerr_str(); if (mes) { @@ -433,20 +409,16 @@ static int hoccommand_exec_strret(Object* ho, char* buf, int size) { PyErr_Print(); } } - return (r != NULL); + return r.is_valid(); } static void grphcmdtool(Object* ho, int type, double x, double y, int key) { - PyObject* po = ((Py2Nrn*) ho->u.this_pointer)->po_; - PyObject* r; + nb::callable po = nb::borrow(((Py2Nrn*) ho->u.this_pointer)->po_); nanobind::gil_scoped_acquire lock{}; - PyObject* args = PyTuple_Pack( - 4, PyInt_FromLong(type), PyFloat_FromDouble(x), PyFloat_FromDouble(y), PyInt_FromLong(key)); - r = nrnpy_pyCallObject(po, args); - Py_XDECREF(args); - Py_XDECREF(r); - if (!r) { + nb::tuple args = nb::make_tuple(type, x, y, key); + nb::object r = nrnpy_pyCallObject(po, args.ptr()); + if (!r.is_valid()) { char* mes = nrnpyerr_str(); if (mes) { Fprintf(stderr, "%s\n", mes); @@ -492,8 +464,7 @@ static Object* callable_with_args(Object* ho, int narg) { } static double func_call(Object* ho, int narg, int* err) { - PyObject* po = ((Py2Nrn*) ho->u.this_pointer)->po_; - PyObject* r; + nb::callable po = nb::borrow(((Py2Nrn*) ho->u.this_pointer)->po_); nanobind::gil_scoped_acquire lock{}; PyObject* args = PyTuple_New((Py_ssize_t) narg); @@ -512,10 +483,10 @@ static double func_call(Object* ho, int narg, int* err) { } } - r = nrnpy_pyCallObject(po, args); + nb::object r = nrnpy_pyCallObject(po, args); Py_XDECREF(args); double rval = 0.0; - if (r == NULL) { + if (!r.is_valid()) { if (!err || *err) { char* mes = nrnpyerr_str(); if (mes) { @@ -535,12 +506,9 @@ static double func_call(Object* ho, int narg, int* err) { *err = 1; } } else { - if (nrnpy_numbercheck(r)) { - PyObject* pn = PyNumber_Float(r); - rval = PyFloat_AsDouble(pn); - Py_XDECREF(pn); + if (nrnpy_numbercheck(r.ptr())) { + rval = static_cast(nb::float_(r)); } - Py_XDECREF(r); if (err) { *err = 0; } // success From 1ef583b31ee1809cd3b24061811427abd52a0919 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Fri, 18 Oct 2024 11:44:42 +0200 Subject: [PATCH 04/17] Avoid warning comming from internal of eigen (#3132) --- src/nrnpython/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/nrnpython/CMakeLists.txt b/src/nrnpython/CMakeLists.txt index 1a9696de9d..c857e3c37b 100644 --- a/src/nrnpython/CMakeLists.txt +++ b/src/nrnpython/CMakeLists.txt @@ -74,7 +74,7 @@ else() target_link_libraries(nrnpython ${NRN_DEFAULT_PYTHON_LIBRARIES}) endif() target_link_libraries(nrnpython fmt::fmt) - target_include_directories(nrnpython PUBLIC ${PROJECT_SOURCE_DIR}/${NRN_3RDPARTY_DIR}/eigen) + target_include_directories(nrnpython SYSTEM PUBLIC ${PROJECT_SOURCE_DIR}/${NRN_3RDPARTY_DIR}/eigen) target_include_directories(nrnpython PUBLIC ${PROJECT_BINARY_DIR}/src/nrniv/oc_generated) make_nanobind_target(nanobind ${NRN_DEFAULT_PYTHON_INCLUDES}) target_link_libraries(nrnpython nanobind) From 27efa9271386b24f3cd7090daa46090be5ef14e1 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Fri, 18 Oct 2024 13:56:00 +0200 Subject: [PATCH 05/17] Remove only a small part (#3133) --- .gitlab-ci.yml | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 6c755ea369..6072f41ea1 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -148,19 +148,6 @@ spack_setup: # BLUECONFIGS_BRANCH does not correspond to a Spack package called blueconfigs SPACK_SETUP_IGNORE_PACKAGE_VARIABLES: BLUECONFIGS -simulation_stack: - stage: .pre - # Take advantage of GitHub PR description parsing in the spack_setup job. - needs: [spack_setup] - trigger: - branch: $BLUECONFIGS_BRANCH - project: hpc/sim/blueconfigs - # NEURON CI status depends on the BlueConfigs CI status. - strategy: depend - variables: - GITLAB_PIPELINES_BRANCH: $GITLAB_PIPELINES_BRANCH - SPACK_ENV_FILE_URL: $SPACK_SETUP_COMMIT_MAPPING_URL - # Performance seems to be terrible when we get too many jobs on a single node. .build: extends: [.spack_build] From 33263d01436d1b89fde299e1c63f3edeca4370ac Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Fri, 18 Oct 2024 16:08:58 +0200 Subject: [PATCH 06/17] Reduce warnings for gnu directory (#3134) --- src/gnu/mcran4.cpp | 15 +-------------- 1 file changed, 1 insertion(+), 14 deletions(-) diff --git a/src/gnu/mcran4.cpp b/src/gnu/mcran4.cpp index a9553687f6..fb913a8a9e 100644 --- a/src/gnu/mcran4.cpp +++ b/src/gnu/mcran4.cpp @@ -92,19 +92,6 @@ uint32_t nrnRan4int(uint32_t* idx1, uint32_t idx2) { /*n ^= (((u >> 16) | (u << 16)) ^ 0x178b0f3c) + w * v;*/ n ^= (((u >> 16) | (u << 16)) ^ 0xe874f0c3) + w * v; return n; - - w = n ^ 0x03bcdc3c; - v = w >> 16; - w &= 0xffff; - u = (v - w) * (v + w); - m ^= (((u >> 16) | (u << 16)) ^ 0x96aa3a59) + w * v; - - w = m ^ 0x0f33d1b2; - v = w >> 16; - w &= 0xffff; - u = (v - w) * (v + w); - n ^= (((u >> 16) | (u << 16)) ^ 0xaa5835b9) + w * v; - return n; } /* @@ -128,7 +115,7 @@ uint32_t nrnRan4int(uint32_t* idx1, uint32_t idx2) { */ static const double SHIFT32 = 1.0 / 4294967296.0; /* 2^-32 */ double nrnRan4dbl(uint32_t* idx1, uint32_t idx2) { - uint32_t hi, lo, extra; + uint32_t hi; hi = (uint32_t) nrnRan4int(idx1, idx2); /*top 32 bits*/ /* // lo = (extra // low bits From cf2318c92fe72ca01e3c3554fdf9917a5c3bb2fd Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Fri, 18 Oct 2024 17:56:04 +0200 Subject: [PATCH 07/17] Remove complex part of sparse13 (#3130) --- src/sparse13/spalloc.cpp | 45 --- src/sparse13/spbuild.cpp | 213 ---------- src/sparse13/spconfig.h | 26 -- src/sparse13/spdefs.h | 248 ------------ src/sparse13/spfactor.cpp | 238 +---------- src/sparse13/spmatrix.h | 50 --- src/sparse13/spoutput.cpp | 98 +---- src/sparse13/spsolve.cpp | 328 ---------------- src/sparse13/sputils.cpp | 809 +------------------------------------- 9 files changed, 23 insertions(+), 2032 deletions(-) diff --git a/src/sparse13/spalloc.cpp b/src/sparse13/spalloc.cpp index ab85ebefd9..233093c61f 100644 --- a/src/sparse13/spalloc.cpp +++ b/src/sparse13/spalloc.cpp @@ -15,7 +15,6 @@ * spWhereSingular * spGetSize * spSetReal - * spSetComplex * spFillinCount * spElementCount * @@ -45,11 +44,6 @@ * or implied warranty. */ -#ifndef lint -static char copyright[] = "Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert"; -static char RCSid[] = "@(#)$Header$"; -#endif - /* * IMPORTS * @@ -126,12 +120,10 @@ char* spCreate(int Size, BOOLEAN Complex, int* pError) } /* Test for valid type. */ -#if NOT spCOMPLEX if (Complex) { *pError = spPANIC; return NULL; } -#endif #if NOT REAL if (NOT Complex) { *pError = spPANIC; @@ -193,9 +185,6 @@ char* spCreate(int Size, BOOLEAN Complex, int* pError) /* Take out the trash. */ Matrix->TrashCan.Real = 0.0; -#if spCOMPLEX - Matrix->TrashCan.Imag = 0.0; -#endif Matrix->TrashCan.Row = 0; Matrix->TrashCan.Col = 0; Matrix->TrashCan.NextInRow = NULL; @@ -233,24 +222,6 @@ char* spCreate(int Size, BOOLEAN Complex, int* pError) Matrix->IntToExtColMap[I] = I; } -#if TRANSLATE - /* Allocate space in memory for ExtToIntColMap vector. */ - if ((Matrix->ExtToIntColMap = ALLOC(int, SizePlusOne)) == NULL) - goto MemoryError; - - /* Allocate space in memory for ExtToIntRowMap vector. */ - if ((Matrix->ExtToIntRowMap = ALLOC(int, SizePlusOne)) == NULL) - goto MemoryError; - - /* Initialize MapExtToInt vectors. */ - for (I = 1; I <= AllocatedSize; I++) { - Matrix->ExtToIntColMap[I] = -1; - Matrix->ExtToIntRowMap[I] = -1; - } - Matrix->ExtToIntColMap[0] = 0; - Matrix->ExtToIntRowMap[0] = 0; -#endif - /* Allocate space for fill-ins and initial set of elements. */ InitializeElementBlocks(Matrix, SPACE_FOR_ELEMENTS * AllocatedSize, SPACE_FOR_FILL_INS * AllocatedSize); @@ -671,14 +642,7 @@ int spGetSize(char* eMatrix, BOOLEAN External) /* Begin `spGetSize'. */ ASSERT(IS_SPARSE(Matrix)); -#if TRANSLATE - if (External) - return Matrix->ExtSize; - else - return Matrix->Size; -#else return Matrix->Size; -#endif } /* @@ -700,15 +664,6 @@ void spSetReal(char* eMatrix) return; } -void spSetComplex(char* eMatrix) -{ - /* Begin `spSetComplex'. */ - - ASSERT(IS_SPARSE((MatrixPtr)eMatrix) AND spCOMPLEX); - ((MatrixPtr)eMatrix)->Complex = YES; - return; -} - /* * ELEMENT OR FILL-IN COUNT * diff --git a/src/sparse13/spbuild.cpp b/src/sparse13/spbuild.cpp index 3b7e9a1da2..f0a5867609 100644 --- a/src/sparse13/spbuild.cpp +++ b/src/sparse13/spbuild.cpp @@ -20,11 +20,9 @@ * * >>> Other functions contained in this file: * spcFindElementInCol - * Translate * spcCreateElement * spcLinkRows * EnlargeMatrix - * ExpandTranslationArrays */ /* @@ -43,11 +41,6 @@ * or implied warranty. */ -#ifndef lint -static char copyright[] = "Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert"; -static char RCSid[] = "@(#)$Header$"; -#endif - /* * IMPORTS * @@ -70,9 +63,7 @@ static char RCSid[] = "@(#)$Header$"; extern ElementPtr spcGetFillin(MatrixPtr Matrix); extern ElementPtr spcGetElement(MatrixPtr Matrix); -static void Translate(MatrixPtr Matrix, int* Row, int* Col); static void EnlargeMatrix(MatrixPtr Matrix, int NewSize); -static void ExpandTranslationArrays(MatrixPtr Matrix, int NewSize); /* * CLEAR MATRIX @@ -98,18 +89,6 @@ void spClear(char* eMatrix) ASSERT(IS_SPARSE(Matrix)); /* Clear matrix. */ -#if spCOMPLEX - if (Matrix->PreviousMatrixWasComplex OR Matrix->Complex) { - for (I = Matrix->Size; I > 0; I--) { - pElement = Matrix->FirstInCol[I]; - while (pElement != NULL) { - pElement->Real = 0.0; - pElement->Imag = 0.0; - pElement = pElement->NextInCol; - } - } - } else -#endif { for (I = Matrix->Size; I > 0; I--) { pElement = Matrix->FirstInCol[I]; @@ -122,9 +101,6 @@ void spClear(char* eMatrix) /* Empty the trash. */ Matrix->TrashCan.Real = 0.0; -#if spCOMPLEX - Matrix->TrashCan.Imag = 0.0; -#endif Matrix->Error = spOKAY; Matrix->Factored = NO; @@ -183,17 +159,8 @@ RealNumber* spGetElement(char* eMatrix, int Row, int Col) if ((Row == 0) OR(Col == 0)) return &Matrix->TrashCan.Real; -#if NOT TRANSLATE ASSERT(Matrix->NeedsOrdering); -#endif - -#if TRANSLATE - Translate(Matrix, &Row, &Col); - if (Matrix->Error == spNO_MEMORY) - return NULL; -#endif -#if NOT TRANSLATE #if NOT EXPANDABLE ASSERT(Row <= Matrix->Size AND Col <= Matrix->Size); #endif @@ -204,7 +171,6 @@ RealNumber* spGetElement(char* eMatrix, int Row, int Col) EnlargeMatrix(Matrix, MAX(Row, Col)); if (Matrix->Error == spNO_MEMORY) return NULL; -#endif #endif /* @@ -290,109 +256,6 @@ ElementPtr spcFindElementInCol(MatrixPtr Matrix, ElementPtr* LastAddr, int Row, return NULL; } -#if TRANSLATE - -/* - * TRANSLATE EXTERNAL INDICES TO INTERNAL - * - * Convert internal row and column numbers to internal row and column numbers. - * Also updates Ext/Int maps. - * - * - * >>> Arguments: - * Matrix (MatrixPtr) - * Pointer to the matrix. - * Row (int *) - * Upon entry Row is either a external row number of an external node - * number. Upon entry, the internal equivalent is supplied. - * Col (int *) - * Upon entry Column is either a external column number of an external node - * number. Upon entry, the internal equivalent is supplied. - * - * >>> Local variables: - * ExtCol (int) - * Temporary variable used to hold the external column or node number - * during the external to internal column number translation. - * ExtRow (int) - * Temporary variable used to hold the external row or node number during - * the external to internal row number translation. - * IntCol (int) - * Temporary variable used to hold the internal column or node number - * during the external to internal column number translation. - * IntRow (int) - * Temporary variable used to hold the internal row or node number during - * the external to internal row number translation. - */ - -static void Translate(MatrixPtr Matrix, int* Row, int* Col) -{ - int IntRow, IntCol, ExtRow, ExtCol; - - /* Begin `Translate'. */ - ExtRow = *Row; - ExtCol = *Col; - - /* Expand translation arrays if necessary. */ - if ((ExtRow > Matrix->AllocatedExtSize) OR(ExtCol > Matrix->AllocatedExtSize)) { - ExpandTranslationArrays(Matrix, MAX(ExtRow, ExtCol)); - if (Matrix->Error == spNO_MEMORY) - return; - } - - /* Set ExtSize if necessary. */ - if ((ExtRow > Matrix->ExtSize) OR(ExtCol > Matrix->ExtSize)) - Matrix->ExtSize = MAX(ExtRow, ExtCol); - - /* Translate external row or node number to internal row or node number. */ - if ((IntRow = Matrix->ExtToIntRowMap[ExtRow]) == -1) { - Matrix->ExtToIntRowMap[ExtRow] = ++Matrix->CurrentSize; - Matrix->ExtToIntColMap[ExtRow] = Matrix->CurrentSize; - IntRow = Matrix->CurrentSize; - -#if NOT EXPANDABLE - ASSERT(IntRow <= Matrix->Size); -#endif - -#if EXPANDABLE - /* Re-size Matrix if necessary. */ - if (IntRow > Matrix->Size) - EnlargeMatrix(Matrix, IntRow); - if (Matrix->Error == spNO_MEMORY) - return; -#endif - - Matrix->IntToExtRowMap[IntRow] = ExtRow; - Matrix->IntToExtColMap[IntRow] = ExtRow; - } - - /* Translate external column or node number to internal column or node number.*/ - if ((IntCol = Matrix->ExtToIntColMap[ExtCol]) == -1) { - Matrix->ExtToIntRowMap[ExtCol] = ++Matrix->CurrentSize; - Matrix->ExtToIntColMap[ExtCol] = Matrix->CurrentSize; - IntCol = Matrix->CurrentSize; - -#if NOT EXPANDABLE - ASSERT(IntCol <= Matrix->Size); -#endif - -#if EXPANDABLE - /* Re-size Matrix if necessary. */ - if (IntCol > Matrix->Size) - EnlargeMatrix(Matrix, IntCol); - if (Matrix->Error == spNO_MEMORY) - return; -#endif - - Matrix->IntToExtRowMap[IntCol] = ExtCol; - Matrix->IntToExtColMap[IntCol] = ExtCol; - } - - *Row = IntRow; - *Col = IntCol; - return; -} -#endif - #if QUAD_ELEMENT /* * ADDITION OF ADMITTANCE TO MATRIX BY INDEX @@ -649,9 +512,6 @@ ElementPtr spcCreateElement(MatrixPtr Matrix, int Row, int Col, ElementPtr* Last pElement->Row = Row; pElement->Col = Col; pElement->Real = 0.0; -#if spCOMPLEX - pElement->Imag = 0.0; -#endif #if INITIALIZE pElement->pInitInfo = NULL; #endif @@ -708,9 +568,6 @@ ElementPtr spcCreateElement(MatrixPtr Matrix, int Row, int Col, ElementPtr* Last pElement->Col = Col; #endif pElement->Real = 0.0; -#if spCOMPLEX - pElement->Imag = 0.0; -#endif #if INITIALIZE pElement->pInitInfo = NULL; #endif @@ -847,58 +704,6 @@ static void EnlargeMatrix(MatrixPtr Matrix, int NewSize) return; } -#if TRANSLATE - -/* - * EXPAND TRANSLATION ARRAYS - * - * Increases the size arrays that are used to translate external to internal - * row and column numbers. - * - * >>> Arguments: - * Matrix (MatrixPtr) - * Pointer to the matrix. - * NewSize (int) - * The new size of the translation arrays. - * - * >>> Local variables: - * OldAllocatedSize (int) - * The allocated size of the translation arrays before being expanded. - */ - -static void ExpandTranslationArrays(MatrixPtr Matrix, int NewSize) -{ - int I, OldAllocatedSize = Matrix->AllocatedExtSize; - - /* Begin `ExpandTranslationArrays'. */ - Matrix->ExtSize = NewSize; - - if (NewSize <= OldAllocatedSize) - return; - - /* Expand the translation arrays ExtToIntRowMap and ExtToIntColMap. */ - NewSize = MAX(NewSize, EXPANSION_FACTOR * OldAllocatedSize); - Matrix->AllocatedExtSize = NewSize; - - if ((REALLOC(Matrix->ExtToIntRowMap, int, NewSize + 1)) == NULL) { - Matrix->Error = spNO_MEMORY; - return; - } - if ((REALLOC(Matrix->ExtToIntColMap, int, NewSize + 1)) == NULL) { - Matrix->Error = spNO_MEMORY; - return; - } - - /* Initialize the new portion of the vectors. */ - for (I = OldAllocatedSize + 1; I <= NewSize; I++) { - Matrix->ExtToIntRowMap[I] = -1; - Matrix->ExtToIntColMap[I] = -1; - } - - return; -} -#endif - #if INITIALIZE /* * INITIALIZE MATRIX @@ -960,18 +765,6 @@ int (*pInit)(); /* Begin `spInitialize'. */ ASSERT(IS_SPARSE(Matrix)); -#if spCOMPLEX - /* Clear imaginary part of matrix if matrix is real but was complex. */ - if (Matrix->PreviousMatrixWasComplex AND NOT Matrix->Complex) { - for (J = Matrix->Size; J > 0; J--) { - pElement = Matrix->FirstInCol[J]; - while (pElement != NULL) { - pElement->Imag = 0.0; - pElement = pElement->NextInCol; - } - } - } -#endif /* spCOMPLEX */ /* Initialize the matrix. */ for (J = Matrix->Size; J > 0; J--) { @@ -980,9 +773,6 @@ int (*pInit)(); while (pElement != NULL) { if (pElement->pInitInfo == NULL) { pElement->Real = 0.0; -#if spCOMPLEX - pElement->Imag = 0.0; -#endif } else { Error = (*pInit)((RealNumber*)pElement, pElement->pInitInfo, Matrix->IntToExtRowMap[pElement->Row], Col); @@ -997,9 +787,6 @@ int (*pInit)(); /* Empty the trash. */ Matrix->TrashCan.Real = 0.0; -#if spCOMPLEX - Matrix->TrashCan.Imag = 0.0; -#endif Matrix->Error = spOKAY; Matrix->Factored = NO; diff --git a/src/sparse13/spconfig.h b/src/sparse13/spconfig.h index ae52e16216..22cd3f9dcf 100644 --- a/src/sparse13/spconfig.h +++ b/src/sparse13/spconfig.h @@ -67,9 +67,6 @@ * both real and complex systems at the same time, but there is a * slight speed and memory advantage if the routines are complied * to handle only real systems of equations. - * spCOMPLEX - * This specifies that the routines will be complied to handle - * complex systems of equations. * EXPANDABLE * Setting this compiler flag true (1) makes the matrix * expandable before it has been factored. If the matrix is @@ -140,16 +137,6 @@ * must have an allocated length of one plus the size of the * matrix. ARRAY_OFFSET must be either 0 or 1, no other offsets * are valid. - * spSEPARATED_COMPLEX_VECTORS - * This specifies the format for complex vectors. If this is set - * false then a complex vector is made up of one double sized - * array of RealNumber's in which the real and imaginary numbers - * are placed in the alternately array in the array. In other - * words, the first entry would be Complex[1].Real, then comes - * Complex[1].Imag, then Complex[1].Real, etc. If - * spSEPARATED_COMPLEX_VECTORS is set true, then each complex - * vector is represented by two arrays of RealNumbers, one with - * the real terms, the other with the imaginary. [NO] * MODIFIED_MARKOWITZ * This specifies that the modified Markowitz method of pivot * selection is to be used. The modified Markowitz method differs @@ -243,12 +230,6 @@ /* Begin options. */ #define REAL YES #define EXPANDABLE YES -#if defined(cmplx_spPrefix) -/* NEURON's nonlinz.cpp uses cmplx_spGetElement after previous use of matrix */ -#define TRANSLATE YES -#else -#define TRANSLATE NO /* instead of YES */ -#endif #define INITIALIZE NO /* instead of YES */ #define DIAGONAL_PIVOTING YES #define ARRAY_OFFSET NOT FORTRAN @@ -275,13 +256,6 @@ * with user code, so use 0 for NO and 1 for YES. */ #endif /* spINSIDE_SPARSE */ -#if defined(cmplx_spPrefix) -#define spCOMPLEX 1 -#define spSEPARATED_COMPLEX_VECTORS 1 -#else -#define spCOMPLEX 0 /* instead of 1 */ -#define spSEPARATED_COMPLEX_VECTORS 0 -#endif #ifdef spINSIDE_SPARSE /* diff --git a/src/sparse13/spdefs.h b/src/sparse13/spdefs.h index f7977c77ad..78ae57ea8b 100644 --- a/src/sparse13/spdefs.h +++ b/src/sparse13/spdefs.h @@ -41,7 +41,6 @@ #ifdef lint #undef REAL -#undef spCOMPLEX #undef EXPANDABLE #undef TRANSLATE #undef INITIALIZE @@ -60,7 +59,6 @@ #undef DEBUG #define REAL YES -#define spCOMPLEX YES #define EXPANDABLE YES #define TRANSLATE YES #define INITIALIZE YES @@ -133,225 +131,7 @@ } /* Macro function that returns the approx absolute value of a complex number. */ -#if spCOMPLEX -#define ELEMENT_MAG(ptr) (ABS((ptr)->Real) + ABS((ptr)->Imag)) -#else #define ELEMENT_MAG(ptr) ((ptr)->Real < 0.0 ? -(ptr)->Real : (ptr)->Real) -#endif - -/* Complex assignment statements. */ -#define CMPLX_ASSIGN(to, from) \ - { \ - (to).Real = (from).Real; \ - (to).Imag = (from).Imag; \ - } -#define CMPLX_CONJ_ASSIGN(to, from) \ - { \ - (to).Real = (from).Real; \ - (to).Imag = -(from).Imag; \ - } -#define CMPLX_NEGATE_ASSIGN(to, from) \ - { \ - (to).Real = -(from).Real; \ - (to).Imag = -(from).Imag; \ - } -#define CMPLX_CONJ_NEGATE_ASSIGN(to, from) \ - { \ - (to).Real = -(from).Real; \ - (to).Imag = (from).Imag; \ - } -#define CMPLX_CONJ(a) (a).Imag = -(a).Imag -#define CMPLX_NEGATE(a) \ - { \ - (a).Real = -(a).Real; \ - (a).Imag = -(a).Imag; \ - } - -/* Macro that returns the approx magnitude (L-1 norm) of a complex number. */ -#define CMPLX_1_NORM(a) (ABS((a).Real) + ABS((a).Imag)) - -/* Macro that returns the approx magnitude (L-infinity norm) of a complex. */ -#define CMPLX_INF_NORM(a) (MAX(ABS((a).Real), ABS((a).Imag))) - -/* Macro function that returns the magnitude (L-2 norm) of a complex number. */ -#define CMPLX_2_NORM(a) (sqrt((a).Real * (a).Real + (a).Imag * (a).Imag)) - -/* Macro function that performs complex addition. */ -#define CMPLX_ADD(to, from_a, from_b) \ - { \ - (to).Real = (from_a).Real + (from_b).Real; \ - (to).Imag = (from_a).Imag + (from_b).Imag; \ - } - -/* Macro function that performs complex subtraction. */ -#define CMPLX_SUBT(to, from_a, from_b) \ - { \ - (to).Real = (from_a).Real - (from_b).Real; \ - (to).Imag = (from_a).Imag - (from_b).Imag; \ - } - -/* Macro function that is equivalent to += operator for complex numbers. */ -#define CMPLX_ADD_ASSIGN(to, from) \ - { \ - (to).Real += (from).Real; \ - (to).Imag += (from).Imag; \ - } - -/* Macro function that is equivalent to -= operator for complex numbers. */ -#define CMPLX_SUBT_ASSIGN(to, from) \ - { \ - (to).Real -= (from).Real; \ - (to).Imag -= (from).Imag; \ - } - -/* Macro function that multiplies a complex number by a scalar. */ -#define SCLR_MULT(to, sclr, cmplx) \ - { \ - (to).Real = (sclr) * (cmplx).Real; \ - (to).Imag = (sclr) * (cmplx).Imag; \ - } - -/* Macro function that multiply-assigns a complex number by a scalar. */ -#define SCLR_MULT_ASSIGN(to, sclr) \ - { \ - (to).Real *= (sclr); \ - (to).Imag *= (sclr); \ - } - -/* Macro function that multiplies two complex numbers. */ -#define CMPLX_MULT(to, from_a, from_b) \ - { \ - (to).Real = (from_a).Real * (from_b).Real - (from_a).Imag * (from_b).Imag; \ - (to).Imag = (from_a).Real * (from_b).Imag + (from_a).Imag * (from_b).Real; \ - } - -/* Macro function that implements to *= from for complex numbers. */ -#define CMPLX_MULT_ASSIGN(to, from) \ - { \ - RealNumber to_real_ = (to).Real; \ - (to).Real = to_real_ * (from).Real - (to).Imag * (from).Imag; \ - (to).Imag = to_real_ * (from).Imag + (to).Imag * (from).Real; \ - } - -/* Macro function that multiplies two complex numbers, the first of which is - * conjugated. */ -#define CMPLX_CONJ_MULT(to, from_a, from_b) \ - { \ - (to).Real = (from_a).Real * (from_b).Real + (from_a).Imag * (from_b).Imag; \ - (to).Imag = (from_a).Real * (from_b).Imag - (from_a).Imag * (from_b).Real; \ - } - -/* Macro function that multiplies two complex numbers and then adds them - * to another. to = add + mult_a * mult_b */ -#define CMPLX_MULT_ADD(to, mult_a, mult_b, add) \ - { \ - (to).Real = (mult_a).Real * (mult_b).Real - (mult_a).Imag * (mult_b).Imag + (add).Real; \ - (to).Imag = (mult_a).Real * (mult_b).Imag + (mult_a).Imag * (mult_b).Real + (add).Imag; \ - } - -/* Macro function that subtracts the product of two complex numbers from - * another. to = subt - mult_a * mult_b */ -#define CMPLX_MULT_SUBT(to, mult_a, mult_b, subt) \ - { \ - (to).Real = (subt).Real - (mult_a).Real * (mult_b).Real + (mult_a).Imag * (mult_b).Imag; \ - (to).Imag = (subt).Imag - (mult_a).Real * (mult_b).Imag - (mult_a).Imag * (mult_b).Real; \ - } - -/* Macro function that multiplies two complex numbers and then adds them - * to another. to = add + mult_a* * mult_b where mult_a* represents mult_a - * conjugate. */ -#define CMPLX_CONJ_MULT_ADD(to, mult_a, mult_b, add) \ - { \ - (to).Real = (mult_a).Real * (mult_b).Real + (mult_a).Imag * (mult_b).Imag + (add).Real; \ - (to).Imag = (mult_a).Real * (mult_b).Imag - (mult_a).Imag * (mult_b).Real + (add).Imag; \ - } - -/* Macro function that multiplies two complex numbers and then adds them - * to another. to += mult_a * mult_b */ -#define CMPLX_MULT_ADD_ASSIGN(to, from_a, from_b) \ - { \ - (to).Real += (from_a).Real * (from_b).Real - (from_a).Imag * (from_b).Imag; \ - (to).Imag += (from_a).Real * (from_b).Imag + (from_a).Imag * (from_b).Real; \ - } - -/* Macro function that multiplies two complex numbers and then subtracts them - * from another. */ -#define CMPLX_MULT_SUBT_ASSIGN(to, from_a, from_b) \ - { \ - (to).Real -= (from_a).Real * (from_b).Real - (from_a).Imag * (from_b).Imag; \ - (to).Imag -= (from_a).Real * (from_b).Imag + (from_a).Imag * (from_b).Real; \ - } - -/* Macro function that multiplies two complex numbers and then adds them - * to the destination. to += from_a* * from_b where from_a* represents from_a - * conjugate. */ -#define CMPLX_CONJ_MULT_ADD_ASSIGN(to, from_a, from_b) \ - { \ - (to).Real += (from_a).Real * (from_b).Real + (from_a).Imag * (from_b).Imag; \ - (to).Imag += (from_a).Real * (from_b).Imag - (from_a).Imag * (from_b).Real; \ - } - -/* Macro function that multiplies two complex numbers and then subtracts them - * from the destination. to -= from_a* * from_b where from_a* represents from_a - * conjugate. */ -#define CMPLX_CONJ_MULT_SUBT_ASSIGN(to, from_a, from_b) \ - { \ - (to).Real -= (from_a).Real * (from_b).Real + (from_a).Imag * (from_b).Imag; \ - (to).Imag -= (from_a).Real * (from_b).Imag - (from_a).Imag * (from_b).Real; \ - } - -/* - * Macro functions that provide complex division. - */ - -/* Complex division: to = num / den */ -#define CMPLX_DIV(to, num, den) \ - { \ - RealNumber r_, s_; \ - if (((den).Real >= (den).Imag AND(den).Real > -(den).Imag) OR((den).Real < (den).Imag AND(den).Real <= -(den).Imag)) { \ - r_ = (den).Imag / (den).Real; \ - s_ = (den).Real + r_ * (den).Imag; \ - (to).Real = ((num).Real + r_ * (num).Imag) / s_; \ - (to).Imag = ((num).Imag - r_ * (num).Real) / s_; \ - } else { \ - r_ = (den).Real / (den).Imag; \ - s_ = (den).Imag + r_ * (den).Real; \ - (to).Real = (r_ * (num).Real + (num).Imag) / s_; \ - (to).Imag = (r_ * (num).Imag - (num).Real) / s_; \ - } \ - } - -/* Complex division and assignment: num /= den */ -#define CMPLX_DIV_ASSIGN(num, den) \ - { \ - RealNumber r_, s_, t_; \ - if (((den).Real >= (den).Imag AND(den).Real > -(den).Imag) OR((den).Real < (den).Imag AND(den).Real <= -(den).Imag)) { \ - r_ = (den).Imag / (den).Real; \ - s_ = (den).Real + r_ * (den).Imag; \ - t_ = ((num).Real + r_ * (num).Imag) / s_; \ - (num).Imag = ((num).Imag - r_ * (num).Real) / s_; \ - (num).Real = t_; \ - } else { \ - r_ = (den).Real / (den).Imag; \ - s_ = (den).Imag + r_ * (den).Real; \ - t_ = (r_ * (num).Real + (num).Imag) / s_; \ - (num).Imag = (r_ * (num).Imag - (num).Real) / s_; \ - (num).Real = t_; \ - } \ - } - -/* Complex reciprocation: to = 1.0 / den */ -#define CMPLX_RECIPROCAL(to, den) \ - { \ - RealNumber r_; \ - if (((den).Real >= (den).Imag AND(den).Real > -(den).Imag) OR((den).Real < (den).Imag AND(den).Real <= -(den).Imag)) { \ - r_ = (den).Imag / (den).Real; \ - (to).Imag = -r_ * ((to).Real = 1.0 / ((den).Real + r_ * (den).Imag)); \ - } else { \ - r_ = (den).Real / (den).Imag; \ - (to).Real = -r_ * ((to).Imag = -1.0 / ((den).Imag + r_ * (den).Real)); \ - } \ - } /* * ASSERT and ABORT @@ -421,26 +201,6 @@ typedef spREAL RealNumber, *RealVector; -/* - * COMPLEX NUMBER DATA STRUCTURE - * - * >>> Structure fields: - * Real (RealNumber) - * The real portion of the number. Real must be the first - * field in this structure. - * Imag (RealNumber) - * The imaginary portion of the number. This field must follow - * immediately after Real. - */ - -/* Begin `ComplexNumber'. */ - -typedef struct -{ - RealNumber Real; - RealNumber Imag; -} ComplexNumber, *ComplexVector; - /* * MATRIX ELEMENT DATA STRUCTURE * @@ -455,11 +215,6 @@ typedef struct * Real (RealNumber) * The real portion of the value of the element. Real must be the first * field in this structure. - * Imag (RealNumber) - * The imaginary portion of the value of the element. If the matrix - * routines are not compiled to handle complex matrices, then this - * field does not exist. If it exists, it must follow immediately after - * Real. * Row (int) * The row number of the element. * Col (int) @@ -488,9 +243,6 @@ typedef struct struct MatrixElement { RealNumber Real; -#if spCOMPLEX - RealNumber Imag; -#endif int Row; int Col; struct MatrixElement* NextInRow; diff --git a/src/sparse13/spfactor.cpp b/src/sparse13/spfactor.cpp index c3f960f6f4..f01c612deb 100644 --- a/src/sparse13/spfactor.cpp +++ b/src/sparse13/spfactor.cpp @@ -13,7 +13,6 @@ * spPartition * * >>> Other functions contained in this file: - * FactorComplexMatrix CreateInternalVectors * CountMarkowitz MarkowitzProducts * SearchForPivot SearchForSingleton * QuicklySearchDiagonal SearchDiagonal @@ -21,7 +20,7 @@ * FindBiggestInColExclude ExchangeRowsAndCols * spcRowExchange spcColExchange * ExchangeColElements ExchangeRowElements - * RealRowColElimination ComplexRowColElimination + * RealRowColElimination * UpdateMarkowitzNumbers CreateFillin * MatrixIsSingular ZeroPivot * WriteStatus @@ -43,11 +42,6 @@ * or implied warranty. */ -#ifndef lint -static char copyright[] = "Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert"; -static char RCSid[] = "@(#)$Header$"; -#endif - /* * IMPORTS * @@ -67,7 +61,6 @@ static char RCSid[] = "@(#)$Header$"; #include /* avoid "declared implicitly `extern' and later `static' " warnings. */ -static int FactorComplexMatrix(MatrixPtr Matrix); static void CreateInternalVectors(MatrixPtr Matrix); static void CountMarkowitz(MatrixPtr Matrix, RealVector RHS, int Step); static void MarkowitzProducts(MatrixPtr Matrix, int Step); @@ -82,7 +75,6 @@ static void ExchangeRowsAndCols(MatrixPtr Matrix, ElementPtr pPivot, int Step); static void ExchangeColElements(MatrixPtr Matrix, int Row1, ElementPtr Element1, int Row2, ElementPtr Element2, int Column); static void ExchangeRowElements(MatrixPtr Matrix, int Col1, ElementPtr Element1, int Col2, ElementPtr Element2, int Row); static void RealRowColElimination(MatrixPtr Matrix, ElementPtr pPivot); -static void ComplexRowColElimination(MatrixPtr Matrix, ElementPtr pPivot); static void UpdateMarkowitzNumbers(MatrixPtr Matrix, ElementPtr pPivot); static ElementPtr CreateFillin(MatrixPtr Matrix, int Row, int Col); static int MatrixIsSingular(MatrixPtr Matrix, int Step); @@ -215,10 +207,7 @@ int spOrderAndFactor(char* eMatrix, RealNumber* RHS, RealNumber RelThreshold, Re pPivot = Matrix->Diag[Step]; LargestInCol = FindLargestInCol(pPivot->NextInCol); if ((LargestInCol * RelThreshold < ELEMENT_MAG(pPivot))) { - if (Matrix->Complex) - ComplexRowColElimination(Matrix, pPivot); - else - RealRowColElimination(Matrix, pPivot); + RealRowColElimination(Matrix, pPivot); } else { ReorderingRequired = YES; break; /* for loop */ @@ -265,10 +254,7 @@ int spOrderAndFactor(char* eMatrix, RealNumber* RHS, RealNumber RelThreshold, Re return MatrixIsSingular(Matrix, Step); ExchangeRowsAndCols(Matrix, pPivot, Step); - if (Matrix->Complex) - ComplexRowColElimination(Matrix, pPivot); - else - RealRowColElimination(Matrix, pPivot); + RealRowColElimination(Matrix, pPivot); if (Matrix->Error >= spFATAL) return Matrix->Error; @@ -335,10 +321,6 @@ int spFactor(char* eMatrix) } if (NOT Matrix->Partitioned) spPartition(eMatrix, spDEFAULT_PARTITION); -#if spCOMPLEX - if (Matrix->Complex) - return FactorComplexMatrix(Matrix); -#endif #if REAL Size = Matrix->Size; @@ -412,117 +394,6 @@ int spFactor(char* eMatrix) #endif /* REAL */ } -#if spCOMPLEX -/* - * FACTOR COMPLEX MATRIX - * - * This routine is the companion routine to spFactor(), it - * handles complex matrices. It is otherwise identical. - * - * >>> Returned: - * The error code is returned. Possible errors are listed below. - * - * >>> Arguments: - * Matrix (char *) - * Pointer to matrix. - * - * >>> Possible errors: - * spSINGULAR - * Error is cleared in this function. - */ - -static int FactorComplexMatrix(MatrixPtr Matrix) -{ - ElementPtr pElement; - ElementPtr pColumn; - int Step, Size; - ComplexNumber Mult, Pivot; - - /* Begin `FactorComplexMatrix'. */ - ASSERT(Matrix->Complex); - - Size = Matrix->Size; - pElement = Matrix->Diag[1]; - if (ELEMENT_MAG(pElement) == 0.0) - return ZeroPivot(Matrix, 1); - /* Cmplx expr: *pPivot = 1.0 / *pPivot. */ - CMPLX_RECIPROCAL(*pElement, *pElement); - - /* Start factorization. */ - for (Step = 2; Step <= Size; Step++) { - if (Matrix->DoCmplxDirect[Step]) { /* Update column using direct addressing scatter-gather. */ - ComplexNumber* Dest; - Dest = (ComplexNumber*)Matrix->Intermediate; - - /* Scatter. */ - pElement = Matrix->FirstInCol[Step]; - while (pElement != NULL) { - Dest[pElement->Row] = *(ComplexNumber*)pElement; - pElement = pElement->NextInCol; - } - - /* Update column. */ - pColumn = Matrix->FirstInCol[Step]; - while (pColumn->Row < Step) { - pElement = Matrix->Diag[pColumn->Row]; - /* Cmplx expr: Mult = Dest[pColumn->Row] * (1.0 / *pPivot). */ - CMPLX_MULT(Mult, Dest[pColumn->Row], *pElement); - CMPLX_ASSIGN(*pColumn, Mult); - while ((pElement = pElement->NextInCol) != NULL) { /* Cmplx expr: Dest[pElement->Row] -= Mult * pElement */ - CMPLX_MULT_SUBT_ASSIGN(Dest[pElement->Row], Mult, *pElement); - } - pColumn = pColumn->NextInCol; - } - - /* Gather. */ - pElement = Matrix->Diag[Step]->NextInCol; - while (pElement != NULL) { - *(ComplexNumber*)pElement = Dest[pElement->Row]; - pElement = pElement->NextInCol; - } - - /* Check for singular matrix. */ - Pivot = Dest[Step]; - if (CMPLX_1_NORM(Pivot) == 0.0) - return ZeroPivot(Matrix, Step); - CMPLX_RECIPROCAL(*Matrix->Diag[Step], Pivot); - } else { /* Update column using direct addressing scatter-gather. */ - ComplexNumber** pDest; - pDest = (ComplexNumber**)Matrix->Intermediate; - - /* Scatter. */ - pElement = Matrix->FirstInCol[Step]; - while (pElement != NULL) { - pDest[pElement->Row] = (ComplexNumber*)pElement; - pElement = pElement->NextInCol; - } - - /* Update column. */ - pColumn = Matrix->FirstInCol[Step]; - while (pColumn->Row < Step) { - pElement = Matrix->Diag[pColumn->Row]; - /* Cmplx expr: Mult = *pDest[pColumn->Row] * (1.0 / *pPivot). */ - CMPLX_MULT(Mult, *pDest[pColumn->Row], *pElement); - CMPLX_ASSIGN(*pDest[pColumn->Row], Mult); - while ((pElement = pElement->NextInCol) != NULL) { /* Cmplx expr: *pDest[pElement->Row] -= Mult * pElement */ - CMPLX_MULT_SUBT_ASSIGN(*pDest[pElement->Row], Mult, *pElement); - } - pColumn = pColumn->NextInCol; - } - - /* Check for singular matrix. */ - pElement = Matrix->Diag[Step]; - if (ELEMENT_MAG(pElement) == 0.0) - return ZeroPivot(Matrix, Step); - CMPLX_RECIPROCAL(*pElement, *pElement); - } - } - - Matrix->Factored = YES; - return (Matrix->Error = spOKAY); -} -#endif /* spCOMPLEX */ - /* * PARTITION MATRIX * @@ -566,7 +437,7 @@ void spPartition(char* eMatrix, int Mode) ElementPtr pElement, pColumn; int Step, Size; int *Nc, *No, *Nm; - BOOLEAN *DoRealDirect, *DoCmplxDirect; + BOOLEAN *DoRealDirect; /* Begin `spPartition'. */ ASSERT(IS_SPARSE(Matrix)); @@ -574,7 +445,6 @@ void spPartition(char* eMatrix, int Mode) return; Size = Matrix->Size; DoRealDirect = Matrix->DoRealDirect; - DoCmplxDirect = Matrix->DoCmplxDirect; Matrix->Partitioned = YES; /* If partition is specified by the user, this is easy. */ @@ -584,18 +454,12 @@ void spPartition(char* eMatrix, int Mode) for (Step = 1; Step <= Size; Step++) #if REAL DoRealDirect[Step] = YES; -#endif -#if spCOMPLEX - DoCmplxDirect[Step] = YES; #endif return; } else if (Mode == spINDIRECT_PARTITION) { for (Step = 1; Step <= Size; Step++) #if REAL DoRealDirect[Step] = NO; -#endif -#if spCOMPLEX - DoCmplxDirect[Step] = NO; #endif return; } else @@ -641,9 +505,6 @@ void spPartition(char* eMatrix, int Mode) #if REAL DoRealDirect[Step] = (Nm[Step] + No[Step] > 3 * Nc[Step] - 2 * Nm[Step]); -#endif -#if spCOMPLEX - DoCmplxDirect[Step] = (Nm[Step] + No[Step] > 7 * Nc[Step] - 4 * Nm[Step]); #endif } @@ -703,25 +564,12 @@ static void CreateInternalVectors(MatrixPtr Matrix) Matrix->Error = spNO_MEMORY; } #endif -#if spCOMPLEX - if (Matrix->DoCmplxDirect == NULL) { - if ((Matrix->DoCmplxDirect = ALLOC(BOOLEAN, Size + 1)) == NULL) - Matrix->Error = spNO_MEMORY; - } -#endif /* Create Intermediate vectors for use in MatrixSolve. */ -#if spCOMPLEX - if (Matrix->Intermediate == NULL) { - if ((Matrix->Intermediate = ALLOC(RealNumber, 2 * (Size + 1))) == NULL) - Matrix->Error = spNO_MEMORY; - } -#else if (Matrix->Intermediate == NULL) { if ((Matrix->Intermediate = ALLOC(RealNumber, Size + 1)) == NULL) Matrix->Error = spNO_MEMORY; } -#endif if (Matrix->Error != spNO_MEMORY) Matrix->InternalVectorsAllocated = YES; @@ -765,17 +613,8 @@ static void CountMarkowitz(MatrixPtr Matrix, RealVector RHS, int Step) /* Correct array pointer for ARRAY_OFFSET. */ #if NOT ARRAY_OFFSET -#if spSEPARATED_COMPLEX_VECTORS OR NOT spCOMPLEX if (RHS != NULL) --RHS; -#else - if (RHS != NULL) { - if (Matrix->Complex) - RHS -= 2; - else - --RHS; - } -#endif #endif /* Generate MarkowitzRow Count for each row. */ @@ -793,19 +632,9 @@ static void CountMarkowitz(MatrixPtr Matrix, RealVector RHS, int Step) /* Include nonzero elements in the RHS vector. */ ExtRow = Matrix->IntToExtRowMap[I]; -#if spSEPARATED_COMPLEX_VECTORS OR NOT spCOMPLEX if (RHS != NULL) if (RHS[ExtRow] != 0.0) Count++; -#else - if (RHS != NULL) { - if (Matrix->Complex) { - if ((RHS[2 * ExtRow] != 0.0) OR(RHS[2 * ExtRow + 1] != 0.0)) - Count++; - } else if (RHS[I] != 0.0) - Count++; - } -#endif Matrix->MarkowitzRow[I] = Count; } @@ -2015,10 +1844,6 @@ void spcRowExchange(MatrixPtr Matrix, int Row1, int Row2) SWAP(int, Matrix->MarkowitzRow[Row1], Matrix->MarkowitzRow[Row2]); SWAP(ElementPtr, Matrix->FirstInRow[Row1], Matrix->FirstInRow[Row2]); SWAP(int, Matrix->IntToExtRowMap[Row1], Matrix->IntToExtRowMap[Row2]); -#if TRANSLATE - Matrix->ExtToIntRowMap[Matrix->IntToExtRowMap[Row1]] = Row1; - Matrix->ExtToIntRowMap[Matrix->IntToExtRowMap[Row2]] = Row2; -#endif return; } @@ -2102,10 +1927,6 @@ void spcColExchange(MatrixPtr Matrix, int Col1, int Col2) SWAP(int, Matrix->MarkowitzCol[Col1], Matrix->MarkowitzCol[Col2]); SWAP(ElementPtr, Matrix->FirstInCol[Col1], Matrix->FirstInCol[Col2]); SWAP(int, Matrix->IntToExtColMap[Col1], Matrix->IntToExtColMap[Col2]); -#if TRANSLATE - Matrix->ExtToIntColMap[Matrix->IntToExtColMap[Col1]] = Col1; - Matrix->ExtToIntColMap[Matrix->IntToExtColMap[Col2]] = Col2; -#endif return; } @@ -2454,57 +2275,6 @@ static void RealRowColElimination(MatrixPtr Matrix, ElementPtr pPivot) * spNO_MEMORY */ -static void ComplexRowColElimination(MatrixPtr Matrix, ElementPtr pPivot) -{ -#if spCOMPLEX - ElementPtr pSub; - int Row; - ElementPtr pLower, pUpper; - - /* Begin `ComplexRowColElimination'. */ - - /* Test for zero pivot. */ - if (ELEMENT_MAG(pPivot) == 0.0) { - (void)MatrixIsSingular(Matrix, pPivot->Row); - return; - } - CMPLX_RECIPROCAL(*pPivot, *pPivot); - - pUpper = pPivot->NextInRow; - while (pUpper != NULL) { - /* Calculate upper triangular element. */ - /* Cmplx expr: *pUpper = *pUpper * (1.0 / *pPivot). */ - CMPLX_MULT_ASSIGN(*pUpper, *pPivot); - - pSub = pUpper->NextInCol; - pLower = pPivot->NextInCol; - while (pLower != NULL) { - Row = pLower->Row; - - /* Find element in row that lines up with current lower triangular element. */ - while (pSub != NULL AND pSub->Row < Row) - pSub = pSub->NextInCol; - - /* Test to see if desired element was not found, if not, create fill-in. */ - if (pSub == NULL OR pSub->Row > Row) { - pSub = CreateFillin(Matrix, Row, pUpper->Col); - if (pSub == NULL) { - Matrix->Error = spNO_MEMORY; - return; - } - } - - /* Cmplx expr: pElement -= *pUpper * pLower. */ - CMPLX_MULT_SUBT_ASSIGN(*pSub, *pUpper, *pLower); - pSub = pSub->NextInCol; - pLower = pLower->NextInCol; - } - pUpper = pUpper->NextInRow; - } - return; -#endif /* spCOMPLEX */ -} - /* * UPDATE MARKOWITZ NUMBERS * diff --git a/src/sparse13/spmatrix.h b/src/sparse13/spmatrix.h index a923809ef8..ce5878c324 100644 --- a/src/sparse13/spmatrix.h +++ b/src/sparse13/spmatrix.h @@ -227,11 +227,6 @@ struct spTemplate { */ /* Begin function declarations. */ - -#if defined(__STDC__) || defined(__cplusplus) - -/* For compilers that understand function prototypes. */ - extern void spClear(char*); extern spREAL spCondition(char*, spREAL, int*); extern char* spCreate(int, int, int*); @@ -260,7 +255,6 @@ extern void spPrint(char*, int, int, int); extern spREAL spPseudoCondition(char*); extern spREAL spRoundoff(char*, spREAL); extern void spScale(char*, spREAL[], spREAL[]); -extern void spSetComplex(char*); extern void spSetReal(char*); extern void spStripFills(char*); extern void spWhereSingular(char*, int*, int*); @@ -273,48 +267,4 @@ extern void spMultTransposed(char* eMatrix, spREAL* RHS, spREAL* Solution, std:: extern void spSolve(char* eMatrix, spREAL* RHS, spREAL* Solution, std::optional iRHS = std::nullopt, std::optional iSolution = std::nullopt); extern void spSolveTransposed(char*, spREAL*, spREAL*, std::optional = std::nullopt, std::optional = std::nullopt); -#else /* NOT defined(__STDC__) */ - -/* For compilers that do not understand function prototypes. */ - -extern void spClear(); -extern spREAL spCondition(); -extern char* spCreate(); -extern void spDeleteRowAndCol(); -extern void spDestroy(); -extern void spDeterminant(); -extern int spElementCount(); -extern int spError(); -extern int spFactor(); -extern int spFileMatrix(); -extern int spFileStats(); -extern int spFileVector(); -extern int spFillinCount(); -extern int spGetAdmittance(); -extern spREAL* spGetElement(char*, int, int); -extern char* spGetInitInfo(); -extern int spGetOnes(); -extern int spGetQuad(); -extern int spGetSize(); -extern int spInitialize(); -extern void spInstallInitInfo(); -extern spREAL spLargestElement(); -extern void spMNA_Preorder(); -extern void spMultiply(); -extern void spMultTransposed(); -extern spREAL spNorm(); -extern int spOrderAndFactor(); -extern void spPartition(); -extern void spPrint(); -extern spREAL spPseudoCondition(); -extern spREAL spRoundoff(); -extern void spScale(); -extern void spSetComplex(); -extern void spSetReal(); -extern void spSolve(); -extern void spSolveTransposed(); -extern void spStripFills(); -extern void spWhereSingular(); -#endif /* defined(__STDC__) */ - #endif /* spOKAY */ diff --git a/src/sparse13/spoutput.cpp b/src/sparse13/spoutput.cpp index 534025bcd2..5c33082f3e 100644 --- a/src/sparse13/spoutput.cpp +++ b/src/sparse13/spoutput.cpp @@ -33,11 +33,6 @@ * or implied warranty. */ -#ifndef lint -static char copyright[] = "Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert"; -static char RCSid[] = "$Header$"; -#endif - /* * IMPORTS * @@ -143,11 +138,7 @@ void spPrint(char* eMatrix, int PrintReordered, int Data, int Header) Size = Matrix->Size; /* Create a packed external to internal row and column translation array. */ -#if TRANSLATE - Top = Matrix->AllocatedExtSize; -#else Top = Matrix->AllocatedSize; -#endif CALLOC(PrintOrdToIntRowMap, int, Top + 1); CALLOC(PrintOrdToIntColMap, int, Top + 1); if (PrintOrdToIntRowMap == NULL OR PrintOrdToIntColMap == NULL) { @@ -286,19 +277,6 @@ void spPrint(char* eMatrix, int PrintReordered, int Data, int Header) } printf("\n"); -#if spCOMPLEX - if (Matrix->Complex AND Data) { - printf(" "); - for (J = StartCol; J <= StopCol; J++) { - if (pImagElements[J - StartCol] != NULL) { - printf(" %8.2lgj", - (double)pImagElements[J - StartCol]->Imag); - } else - printf(" "); - } - printf("\n"); - } -#endif /* spCOMPLEX */ } /* Calculate index of first column in next group. */ @@ -440,31 +418,6 @@ int spFileMatrix(char* eMatrix, char* File, char* Label, int Reordered, int Data return 0; } -#if spCOMPLEX - if (Data AND Matrix->Complex) { - for (I = 1; I <= Size; I++) { - pElement = Matrix->FirstInCol[I]; - while (pElement != NULL) { - if (Reordered) { - Row = pElement->Row; - Col = I; - } else { - Row = Matrix->IntToExtRowMap[pElement->Row]; - Col = Matrix->IntToExtColMap[I]; - } - Err = fprintf(pMatrixFile, "%d\t%d\t%-.15lg\t%-.15lg\n", - Row, Col, (double)pElement->Real, (double)pElement->Imag); - if (Err < 0) - return 0; - pElement = pElement->NextInCol; - } - } - /* Output terminator, a line of zeros. */ - if (Header) - if (fprintf(pMatrixFile, "0\t0\t0.0\t0.0\n") < 0) - return 0; - } -#endif /* spCOMPLEX */ #if REAL if (Data AND NOT Matrix->Complex) { @@ -511,11 +464,9 @@ int spFileMatrix(char* eMatrix, char* File, char* Label, int Reordered, int Data * File (char *) * Name of file into which matrix is to be written. * RHS (RealNumber []) - * Right-hand side vector. This is only the real portion if - * spSEPARATED_COMPLEX_VECTORS is true. + * Right-hand side vector. * iRHS (RealNumber []) - * Right-hand side vector, imaginary portion. Not necessary if matrix - * is real or if spSEPARATED_COMPLEX_VECTORS is set false. + * Right-hand side vector, imaginary portion. * * >>> Local variables: * pMatrixFile (FILE *) @@ -523,17 +474,12 @@ int spFileMatrix(char* eMatrix, char* File, char* Label, int Reordered, int Data * Size (int) * The size of the matrix. * - * >>> Obscure Macros - * IMAG_RHS - * Replaces itself with `, iRHS' if the options spCOMPLEX and - * spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears - * without a trace. */ int spFileVector(char* eMatrix, char* File, RealVector RHS, RealVector iRHS) { MatrixPtr Matrix = (MatrixPtr)eMatrix; - int I, Size, Err; + int I, Size; FILE* pMatrixFile; /* Begin `spFileVector'. */ @@ -545,44 +491,11 @@ int spFileVector(char* eMatrix, char* File, RealVector RHS, RealVector iRHS) /* Correct array pointers for ARRAY_OFFSET. */ #if NOT ARRAY_OFFSET -#if spCOMPLEX - if (Matrix->Complex) { -#if spSEPARATED_COMPLEX_VECTORS - ASSERT(iRHS != NULL) - --RHS; - --iRHS; -#else - RHS -= 2; -#endif - } else -#endif /* spCOMPLEX */ --RHS; #endif /* NOT ARRAY_OFFSET */ /* Output vector. */ Size = Matrix->Size; -#if spCOMPLEX - if (Matrix->Complex) { -#if spSEPARATED_COMPLEX_VECTORS - for (I = 1; I <= Size; I++) { - Err = fprintf(pMatrixFile, "%-.15lg\t%-.15lg\n", - (double)RHS[I], (double)iRHS[I]); - if (Err < 0) - return 0; - } -#else - for (I = 1; I <= Size; I++) { - Err = fprintf(pMatrixFile, "%-.15lg\t%-.15lg\n", - (double)RHS[2 * I], (double)RHS[2 * I + 1]); - if (Err < 0) - return 0; - } -#endif - } -#endif /* spCOMPLEX */ -#if REAL AND spCOMPLEX - else -#endif #if REAL { for (I = 1; I <= Size; I++) { @@ -656,10 +569,7 @@ int spFileStats(char* eMatrix, char* File, char* Label) fprintf(pStatsFile, "Matrix has not been factored.\n"); fprintf(pStatsFile, "||| Starting new matrix |||\n"); fprintf(pStatsFile, "%s\n", Label); - if (Matrix->Complex) - fprintf(pStatsFile, "Matrix is complex.\n"); - else - fprintf(pStatsFile, "Matrix is real.\n"); + fprintf(pStatsFile, "Matrix is real.\n"); fprintf(pStatsFile, " Size = %d\n", Size); /* Search matrix. */ diff --git a/src/sparse13/spsolve.cpp b/src/sparse13/spsolve.cpp index ee70d3b769..95bc4d436c 100644 --- a/src/sparse13/spsolve.cpp +++ b/src/sparse13/spsolve.cpp @@ -11,10 +11,6 @@ * >>> User accessible functions contained in this file: * spSolve * spSolveTransposed - * - * >>> Other functions contained in this file: - * SolveComplexMatrix - * SolveComplexTransposedMatrix */ /* @@ -33,11 +29,6 @@ * or implied warranty. */ -#ifndef lint -static char copyright[] = "Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert"; -static char RCSid[] = "@(#)$Header$"; -#endif - /* * IMPORTS * @@ -55,10 +46,6 @@ static char RCSid[] = "@(#)$Header$"; #include "spdefs.h" #include "spmatrix.h" -/* avoid "declared implicitly `extern' and later `static' " warnings. */ -static void SolveComplexMatrix(MatrixPtr Matrix, RealVector RHS, RealVector Solution, std::optional iRHS = std::nullopt, std::optional iSolution = std::nullopt); -static void SolveComplexTransposedMatrix(MatrixPtr Matrix, RealVector RHS, RealVector Solution, std::optional iRHS = std::nullopt, std::optional iSolution = std::nullopt); - /* * SOLVE MATRIX EQUATION * @@ -126,12 +113,6 @@ void spSolve(char* eMatrix, RealVector RHS, RealVector Solution, std::optionalComplex) { - SolveComplexMatrix(Matrix, RHS, Solution, iRHS, iSolution); - return; - } -#endif #if REAL Intermediate = Matrix->Intermediate; @@ -183,155 +164,6 @@ void spSolve(char* eMatrix, RealVector RHS, RealVector Solution, std::optional>> Arguments: - * Matrix (char *) - * Pointer to matrix. - * RHS (RealVector) - * RHS is the real portion of the input data array, the right hand - * side. This data is undisturbed and may be reused for other solves. - * Solution (RealVector) - * Solution is the real portion of the output data array. This routine - * is constructed such that RHS and Solution can be the same - * array. - * iRHS (RealVector) - * iRHS is the imaginary portion of the input data array, the right - * hand side. This data is undisturbed and may be reused for other solves. - * If spSEPARATED_COMPLEX_VECTOR is set false, there is no need to - * supply this array. - * iSolution (RealVector) - * iSolution is the imaginary portion of the output data array. This - * routine is constructed such that iRHS and iSolution can be - * the same array. If spSEPARATED_COMPLEX_VECTOR is set false, there is no - * need to supply this array. - * - * >>> Local variables: - * Intermediate (ComplexVector) - * Temporary storage for use in forward elimination and backward - * substitution. Commonly referred to as c, when the LU factorization - * equations are given as Ax = b, Lc = b, Ux = c. - * Local version of Matrix->Intermediate, which was created during - * the initial factorization in function CreateInternalVectors() in the - * matrix factorization module. - * pElement (ElementPtr) - * Pointer used to address elements in both the lower and upper triangle - * matrices. - * pExtOrder (int *) - * Pointer used to sequentially access each entry in IntToExtRowMap - * and IntToExtColMap arrays. Used to quickly scramble and unscramble - * RHS and Solution to account for row and column interchanges. - * pPivot (ElementPtr) - * Pointer that points to current pivot or diagonal element. - * Size (int) - * Size of matrix. Made local to reduce indirection. - * Temp (ComplexNumber) - * Temporary storage for entries in arrays. - */ - -static void SolveComplexMatrix(MatrixPtr Matrix, RealVector RHS, RealVector Solution, std::optional iRHS, std::optional iSolution) -{ - ElementPtr pElement; - ComplexVector Intermediate; - int I, *pExtOrder, Size; - ElementPtr pPivot; - ComplexVector ExtVector; - ComplexNumber Temp; - - /* Begin `SolveComplexMatrix'. */ - - Size = Matrix->Size; - Intermediate = (ComplexVector)Matrix->Intermediate; - -/* Correct array pointers for ARRAY_OFFSET. */ -#if NOT ARRAY_OFFSET -#if spSEPARATED_COMPLEX_VECTORS - --RHS; - --iRHS; - --Solution; - --iSolution; -#else - RHS -= 2; - Solution -= 2; -#endif -#endif - - /* Initialize Intermediate vector. */ - pExtOrder = &Matrix->IntToExtRowMap[Size]; - -#if spSEPARATED_COMPLEX_VECTORS - for (I = Size; I > 0; I--) { - Intermediate[I].Real = RHS[*(pExtOrder)]; - ASSERT(iRHS.has_value()); - Intermediate[I].Imag = iRHS.value()[*(pExtOrder--)]; - } -#else - ExtVector = (ComplexVector)RHS; - for (I = Size; I > 0; I--) - Intermediate[I] = ExtVector[*(pExtOrder--)]; -#endif - - /* Forward substitution. Solves Lc = b.*/ - for (I = 1; I <= Size; I++) { - Temp = Intermediate[I]; - - /* This step of the substitution is skipped if Temp equals zero. */ - if ((Temp.Real != 0.0) OR(Temp.Imag != 0.0)) { - pPivot = Matrix->Diag[I]; - /* Cmplx expr: Temp *= (1.0 / Pivot). */ - CMPLX_MULT_ASSIGN(Temp, *pPivot); - Intermediate[I] = Temp; - pElement = pPivot->NextInCol; - while (pElement != NULL) { - /* Cmplx expr: Intermediate[Element->Row] -= Temp * *Element. */ - CMPLX_MULT_SUBT_ASSIGN(Intermediate[pElement->Row], - Temp, *pElement); - pElement = pElement->NextInCol; - } - } - } - - /* Backward Substitution. Solves Ux = c.*/ - for (I = Size; I > 0; I--) { - Temp = Intermediate[I]; - pElement = Matrix->Diag[I]->NextInRow; - - while (pElement != NULL) { - /* Cmplx expr: Temp -= *Element * Intermediate[Element->Col]. */ - CMPLX_MULT_SUBT_ASSIGN(Temp, *pElement, Intermediate[pElement->Col]); - pElement = pElement->NextInRow; - } - Intermediate[I] = Temp; - } - - /* Unscramble Intermediate vector while placing data in to Solution vector. */ - pExtOrder = &Matrix->IntToExtColMap[Size]; - -#if spSEPARATED_COMPLEX_VECTORS - for (I = Size; I > 0; I--) { - Solution[*(pExtOrder)] = Intermediate[I].Real; - iSolution.value()[*(pExtOrder--)] = Intermediate[I].Imag; - } -#else - ExtVector = (ComplexVector)Solution; - for (I = Size; I > 0; I--) - ExtVector[*(pExtOrder--)] = Intermediate[I]; -#endif - - return; -} -#endif /* spCOMPLEX */ - #if TRANSPOSE /* * SOLVE TRANSPOSED MATRIX EQUATION @@ -402,13 +234,6 @@ void spSolveTransposed(char* eMatrix, RealVector RHS, RealVector Solution, std:: /* Begin `spSolveTransposed'. */ ASSERT(IS_VALID(Matrix) AND IS_FACTORED(Matrix)); -#if spCOMPLEX - if (Matrix->Complex) { - SolveComplexTransposedMatrix(Matrix, RHS, Solution, iRHS, iSolution); - return; - } -#endif - #if REAL Size = Matrix->Size; Intermediate = Matrix->Intermediate; @@ -457,156 +282,3 @@ void spSolveTransposed(char* eMatrix, RealVector RHS, RealVector Solution, std:: #endif /* REAL */ } #endif /* TRANSPOSE */ - -#if TRANSPOSE AND spCOMPLEX -/* - * SOLVE COMPLEX TRANSPOSED MATRIX EQUATION - * - * Performs forward elimination and back substitution to find the - * unknown vector from the RHS vector and transposed factored - * matrix. This routine is useful when performing sensitivity analysis - * on a circuit using the adjoint method. This routine assumes that - * the pivots are associated with the untransposed lower triangular - * (L) matrix and that the diagonal of the untransposed upper - * triangular (U) matrix consists of ones. - * - * >>> Arguments: - * Matrix (char *) - * Pointer to matrix. - * RHS (RealVector) - * RHS is the input data array, the right hand - * side. This data is undisturbed and may be reused for other solves. - * This vector is only the real portion if the matrix is complex and - * spSEPARATED_COMPLEX_VECTORS is set true. - * Solution (RealVector) - * Solution is the real portion of the output data array. This routine - * is constructed such that RHS and Solution can be the same array. - * This vector is only the real portion if the matrix is complex and - * spSEPARATED_COMPLEX_VECTORS is set true. - * iRHS (RealVector) - * iRHS is the imaginary portion of the input data array, the right - * hand side. This data is undisturbed and may be reused for other solves. - * If either spCOMPLEX or spSEPARATED_COMPLEX_VECTOR is set false, there - * is no need to supply this array. - * iSolution (RealVector) - * iSolution is the imaginary portion of the output data array. This - * routine is constructed such that iRHS and iSolution can be - * the same array. If spCOMPLEX or spSEPARATED_COMPLEX_VECTOR is set - * false, there is no need to supply this array. - * - * >>> Local variables: - * Intermediate (ComplexVector) - * Temporary storage for use in forward elimination and backward - * substitution. Commonly referred to as c, when the LU factorization - * equations are given as Ax = b, Lc = b, Ux = c. Local version of - * Matrix->Intermediate, which was created during - * the initial factorization in function CreateInternalVectors() in the - * matrix factorization module. - * pElement (ElementPtr) - * Pointer used to address elements in both the lower and upper triangle - * matrices. - * pExtOrder (int *) - * Pointer used to sequentially access each entry in IntToExtRowMap - * and IntToExtColMap arrays. Used to quickly scramble and unscramble - * RHS and Solution to account for row and column interchanges. - * pPivot (ElementPtr) - * Pointer that points to current pivot or diagonal element. - * Size (int) - * Size of matrix. Made local to reduce indirection. - * Temp (ComplexNumber) - * Temporary storage for entries in arrays. - * - */ - -static void -SolveComplexTransposedMatrix(MatrixPtr Matrix, RealVector RHS, RealVector Solution, std::optional iRHS, std::optional iSolution) -{ - ElementPtr pElement; - ComplexVector Intermediate; - int I, *pExtOrder, Size; - ComplexVector ExtVector; - ElementPtr pPivot; - ComplexNumber Temp; - - /* Begin `SolveComplexTransposedMatrix'. */ - - Size = Matrix->Size; - Intermediate = (ComplexVector)Matrix->Intermediate; - -/* Correct array pointers for ARRAY_OFFSET. */ -#if NOT ARRAY_OFFSET -#if spSEPARATED_COMPLEX_VECTORS - --RHS; - --iRHS; - --Solution; - --iSolution; -#else - RHS -= 2; - Solution -= 2; -#endif -#endif - - /* Initialize Intermediate vector. */ - pExtOrder = &Matrix->IntToExtColMap[Size]; - -#if spSEPARATED_COMPLEX_VECTORS - for (I = Size; I > 0; I--) { - Intermediate[I].Real = RHS[*(pExtOrder)]; - ASSERT(iRHS.has_value()); - Intermediate[I].Imag = iRHS.value()[*(pExtOrder--)]; - } -#else - ExtVector = (ComplexVector)RHS; - for (I = Size; I > 0; I--) - Intermediate[I] = ExtVector[*(pExtOrder--)]; -#endif - - /* Forward elimination. */ - for (I = 1; I <= Size; I++) { - Temp = Intermediate[I]; - - /* This step of the elimination is skipped if Temp equals zero. */ - if ((Temp.Real != 0.0) OR(Temp.Imag != 0.0)) { - pElement = Matrix->Diag[I]->NextInRow; - while (pElement != NULL) { - /* Cmplx expr: Intermediate[Element->Col] -= Temp * *Element. */ - CMPLX_MULT_SUBT_ASSIGN(Intermediate[pElement->Col], - Temp, *pElement); - pElement = pElement->NextInRow; - } - } - } - - /* Backward Substitution. */ - for (I = Size; I > 0; I--) { - pPivot = Matrix->Diag[I]; - Temp = Intermediate[I]; - pElement = pPivot->NextInCol; - - while (pElement != NULL) { - /* Cmplx expr: Temp -= Intermediate[Element->Row] * *Element. */ - CMPLX_MULT_SUBT_ASSIGN(Temp, Intermediate[pElement->Row], *pElement); - - pElement = pElement->NextInCol; - } - /* Cmplx expr: Intermediate = Temp * (1.0 / *pPivot). */ - CMPLX_MULT(Intermediate[I], Temp, *pPivot); - } - - /* Unscramble Intermediate vector while placing data in to Solution vector. */ - pExtOrder = &Matrix->IntToExtRowMap[Size]; - -#if spSEPARATED_COMPLEX_VECTORS - for (I = Size; I > 0; I--) { - Solution[*(pExtOrder)] = Intermediate[I].Real; - iSolution.value()[*(pExtOrder--)] = Intermediate[I].Imag; - } -#else - ExtVector = (ComplexVector)Solution; - for (I = Size; I > 0; I--) - ExtVector[*(pExtOrder--)] = Intermediate[I]; -#endif - - return; -} -#endif /* TRANSPOSE AND spCOMPLEX */ diff --git a/src/sparse13/sputils.cpp b/src/sparse13/sputils.cpp index 66984fef77..613d10216f 100644 --- a/src/sparse13/sputils.cpp +++ b/src/sparse13/sputils.cpp @@ -24,9 +24,6 @@ * >>> Other functions contained in this file: * CountTwins * SwapCols - * ScaleComplexMatrix - * ComplexMatrixMultiply - * ComplexCondition */ /* @@ -45,11 +42,6 @@ * or implied warranty. */ -#ifndef lint -static char copyright[] = "Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert"; -static char RCSid[] = "@(#)$Header$"; -#endif - /* * IMPORTS * @@ -76,10 +68,6 @@ extern ElementPtr spcFindElementInCol(MatrixPtr Matrix, ElementPtr* LastAddr, in /* avoid "declared implicitly `extern' and later `static' " warnings. */ static int CountTwins(MatrixPtr Matrix, int Col, ElementPtr* ppTwin1, ElementPtr* ppTwin2); static void SwapCols(MatrixPtr Matrix, ElementPtr pTwin1, ElementPtr pTwin2); -static void ScaleComplexMatrix(MatrixPtr Matrix, RealVector RHS_ScaleFactors, RealVector SolutionScaleFactors); -static void ComplexMatrixMultiply(MatrixPtr Matrix, RealVector RHS, RealVector Solution, std::optional iRHS = std::nullopt, std::optional iSolution = std::nullopt); -static void ComplexTransposedMatrixMultiply(MatrixPtr Matrix, RealVector RHS, RealVector Solution, std::optional iRHS = std::nullopt, std::optional iSolution = std::nullopt); -static RealNumber ComplexCondition(MatrixPtr Matrix, RealNumber NormOfMatrix, int* pError); #if MODIFIED_NODAL /* @@ -267,10 +255,6 @@ static void SwapCols(MatrixPtr Matrix, ElementPtr pTwin1, ElementPtr pTwin2) SWAP(ElementPtr, Matrix->FirstInCol[Col1], Matrix->FirstInCol[Col2]); SWAP(int, Matrix->IntToExtColMap[Col1], Matrix->IntToExtColMap[Col2]); -#if TRANSLATE - Matrix->ExtToIntColMap[Matrix->IntToExtColMap[Col2]] = Col2; - Matrix->ExtToIntColMap[Matrix->IntToExtColMap[Col1]] = Col1; -#endif Matrix->Diag[Col1] = pTwin2; Matrix->Diag[Col2] = pTwin1; @@ -348,13 +332,6 @@ void spScale(char* eMatrix, RealVector RHS_ScaleFactors, RealVector SolutionScal if (NOT Matrix->RowsLinked) spcLinkRows(Matrix); -#if spCOMPLEX - if (Matrix->Complex) { - ScaleComplexMatrix(Matrix, RHS_ScaleFactors, SolutionScaleFactors); - return; - } -#endif - #if REAL lSize = Matrix->Size; @@ -395,108 +372,6 @@ void spScale(char* eMatrix, RealVector RHS_ScaleFactors, RealVector SolutionScal } #endif /* SCALING */ -#if spCOMPLEX AND SCALING -/* - * SCALE COMPLEX MATRIX - * - * This function scales the matrix to enhance the possibility of - * finding a good pivoting order. Note that scaling enhances accuracy - * of the solution only if it affects the pivoting order, so it makes - * no sense to scale the matrix before spFactor(). If scaling is - * desired it should be done before spOrderAndFactor(). There - * are several things to take into account when choosing the scale - * factors. First, the scale factors are directly multiplied against - * the elements in the matrix. To prevent roundoff, each scale factor - * should be equal to an integer power of the number base of the - * machine. Since most machines operate in base two, scale factors - * should be a power of two. Second, the matrix should be scaled such - * that the matrix of element uncertainties is equilibrated. Third, - * this function multiplies the scale factors by the elements, so if - * one row tends to have uncertainties 1000 times smaller than the - * other rows, then its scale factor should be 1024, not 1/1024. - * Fourth, to save time, this function does not scale rows or columns - * if their scale factors are equal to one. Thus, the scale factors - * should be normalized to the most common scale factor. Rows and - * columns should be normalized separately. For example, if the size - * of the matrix is 100 and 10 rows tend to have uncertainties near - * 1e-6 and the remaining 90 have uncertainties near 1e-12, then the - * scale factor for the 10 should be 1/1,048,576 and the scale factors - * for the remaining 90 should be 1. Fifth, since this routine - * directly operates on the matrix, it is necessary to apply the scale - * factors to the RHS and Solution vectors. It may be easier to - * simply use spOrderAndFactor() on a scaled matrix to choose the - * pivoting order, and then throw away the matrix. Subsequent - * factorizations, performed with spFactor(), will not need to have - * the RHS and Solution vectors descaled. Lastly, this function - * should not be executed before the function spMNA_Preorder. - * - * >>> Arguments: - * Matrix (char *) - * Pointer to the matrix to be scaled. - * SolutionScaleFactors (RealVector) - * The array of Solution scale factors. These factors scale the columns. - * All scale factors are real valued. - * RHS_ScaleFactors (RealVector) - * The array of RHS scale factors. These factors scale the rows. - * All scale factors are real valued. - * - * >>> Local variables: - * lSize (int) - * Local version of the size of the matrix. - * pElement (ElementPtr) - * Pointer to an element in the matrix. - * pExtOrder (int *) - * Pointer into either IntToExtRowMap or IntToExtColMap vector. Used to - * compensate for any row or column swaps that have been performed. - * ScaleFactor (RealNumber) - * The scale factor being used on the current row or column. - */ - -static void ScaleComplexMatrix(MatrixPtr Matrix, RealVector RHS_ScaleFactors, RealVector SolutionScaleFactors) -{ - ElementPtr pElement; - int I, lSize, *pExtOrder; - RealNumber ScaleFactor; - - /* Begin `ScaleComplexMatrix'. */ - lSize = Matrix->Size; - -/* Correct pointers to arrays for ARRAY_OFFSET */ -#if NOT ARRAY_OFFSET - --RHS_ScaleFactors; - --SolutionScaleFactors; -#endif - - /* Scale Rows */ - pExtOrder = &Matrix->IntToExtRowMap[1]; - for (I = 1; I <= lSize; I++) { - if ((ScaleFactor = RHS_ScaleFactors[*(pExtOrder++)]) != 1.0) { - pElement = Matrix->FirstInRow[I]; - - while (pElement != NULL) { - pElement->Real *= ScaleFactor; - pElement->Imag *= ScaleFactor; - pElement = pElement->NextInRow; - } - } - } - - /* Scale Columns */ - pExtOrder = &Matrix->IntToExtColMap[1]; - for (I = 1; I <= lSize; I++) { - if ((ScaleFactor = SolutionScaleFactors[*(pExtOrder++)]) != 1.0) { - pElement = Matrix->FirstInCol[I]; - - while (pElement != NULL) { - pElement->Real *= ScaleFactor; - pElement->Imag *= ScaleFactor; - pElement = pElement->NextInCol; - } - } - } - return; -} -#endif /* SCALING AND spCOMPLEX */ #if MULTIPLICATION /* @@ -516,12 +391,10 @@ static void ScaleComplexMatrix(MatrixPtr Matrix, RealVector RHS_ScaleFactors, Re * Solution is the vector being multiplied by the matrix. * iRHS (RealVector) * iRHS is the imaginary portion of the right hand side. This is - * what is being solved for. This is only necessary if the matrix is - * complex and spSEPARATED_COMPLEX_VECTORS is true. + * what is being solved for. * iSolution (RealVector) * iSolution is the imaginary portion of the vector being multiplied - * by the matrix. This is only necessary if the matrix is - * complex and spSEPARATED_COMPLEX_VECTORS is true. + * by the matrix. * */ @@ -538,13 +411,6 @@ void spMultiply(char* eMatrix, RealVector RHS, RealVector Solution, std::optiona if (NOT Matrix->RowsLinked) spcLinkRows(Matrix); -#if spCOMPLEX - if (Matrix->Complex) { - ComplexMatrixMultiply(Matrix, RHS, Solution, iRHS, iSolution); - return; - } -#endif - #if REAL #if NOT ARRAY_OFFSET /* Correct array pointers for ARRAY_OFFSET. */ @@ -574,95 +440,6 @@ void spMultiply(char* eMatrix, RealVector RHS, RealVector Solution, std::optiona } #endif /* MULTIPLICATION */ -#if spCOMPLEX AND MULTIPLICATION -/* - * COMPLEX MATRIX MULTIPLICATION - * - * Multiplies matrix by solution vector to find source vector. - * Assumes matrix has not been factored. This routine can be used - * as a test to see if solutions are correct. - * - * >>> Arguments: - * Matrix (char *) - * Pointer to the matrix. - * RHS (RealVector) - * RHS is the right hand side. This is what is being solved for. - * This is only the real portion of the right-hand side if the matrix - * is complex and spSEPARATED_COMPLEX_VECTORS is set true. - * Solution (RealVector) - * Solution is the vector being multiplied by the matrix. This is only - * the real portion if the matrix is complex and - * spSEPARATED_COMPLEX_VECTORS is set true. - * iRHS (RealVector) - * iRHS is the imaginary portion of the right hand side. This is - * what is being solved for. This is only necessary if the matrix is - * complex and spSEPARATED_COMPLEX_VECTORS is true. - * iSolution (RealVector) - * iSolution is the imaginary portion of the vector being multiplied - * by the matrix. This is only necessary if the matrix is - * complex and spSEPARATED_COMPLEX_VECTORS is true. - * - */ - -static void ComplexMatrixMultiply(MatrixPtr Matrix, RealVector RHS, RealVector Solution, std::optional iRHS, std::optional iSolution) -{ - ElementPtr pElement; - ComplexVector Vector; - ComplexNumber Sum; - int I, *pExtOrder; - -/* Begin `ComplexMatrixMultiply'. */ - -/* Correct array pointers for ARRAY_OFFSET. */ -#if NOT ARRAY_OFFSET -#if spSEPARATED_COMPLEX_VECTORS - --RHS; - --iRHS; - --Solution; - --iSolution; -#else - RHS -= 2; - Solution -= 2; -#endif -#endif - - /* Initialize Intermediate vector with reordered Solution vector. */ - Vector = (ComplexVector)Matrix->Intermediate; - pExtOrder = &Matrix->IntToExtColMap[Matrix->Size]; - -#if spSEPARATED_COMPLEX_VECTORS - for (I = Matrix->Size; I > 0; I--) { - Vector[I].Real = Solution[*pExtOrder]; - ASSERT(iSolution.has_value()); - Vector[I].Imag = iSolution.value()[*(pExtOrder--)]; - } -#else - for (I = Matrix->Size; I > 0; I--) - Vector[I] = ((ComplexVector)Solution)[*(pExtOrder--)]; -#endif - - pExtOrder = &Matrix->IntToExtRowMap[Matrix->Size]; - for (I = Matrix->Size; I > 0; I--) { - pElement = Matrix->FirstInRow[I]; - Sum.Real = Sum.Imag = 0.0; - - while (pElement != NULL) { /* Cmplx expression : Sum += Element * Vector[Col] */ - CMPLX_MULT_ADD_ASSIGN(Sum, *pElement, Vector[pElement->Col]); - pElement = pElement->NextInRow; - } - -#if spSEPARATED_COMPLEX_VECTORS - RHS[*pExtOrder] = Sum.Real; - ASSERT(iRHS != std::nullopt); - iRHS.value()[*pExtOrder--] = Sum.Imag; -#else - ((ComplexVector)RHS)[*pExtOrder--] = Sum; -#endif - } - return; -} -#endif /* spCOMPLEX AND MULTIPLICATION */ - #if MULTIPLICATION AND TRANSPOSE /* * TRANSPOSED MATRIX MULTIPLICATION @@ -681,12 +458,10 @@ static void ComplexMatrixMultiply(MatrixPtr Matrix, RealVector RHS, RealVector S * Solution is the vector being multiplied by the matrix. * iRHS (RealVector) * iRHS is the imaginary portion of the right hand side. This is - * what is being solved for. This is only necessary if the matrix is - * complex and spSEPARATED_COMPLEX_VECTORS is true. + * what is being solved for. * iSolution (RealVector) * iSolution is the imaginary portion of the vector being multiplied - * by the matrix. This is only necessary if the matrix is - * complex and spSEPARATED_COMPLEX_VECTORS is true. + * by the matrix. * */ @@ -701,12 +476,6 @@ void spMultTransposed(char* eMatrix, RealVector RHS, RealVector Solution, std::o /* Begin `spMultTransposed'. */ ASSERT(IS_SPARSE(Matrix) AND NOT Matrix->Factored); -#if spCOMPLEX - if (Matrix->Complex) { - ComplexTransposedMatrixMultiply(Matrix, RHS, Solution, iRHS, iSolution); - return; - } -#endif #if REAL #if NOT ARRAY_OFFSET @@ -737,94 +506,6 @@ void spMultTransposed(char* eMatrix, RealVector RHS, RealVector Solution, std::o } #endif /* MULTIPLICATION AND TRANSPOSE */ -#if spCOMPLEX AND MULTIPLICATION AND TRANSPOSE -/* - * COMPLEX TRANSPOSED MATRIX MULTIPLICATION - * - * Multiplies transposed matrix by solution vector to find source vector. - * Assumes matrix has not been factored. This routine can be used - * as a test to see if solutions are correct. - * - * >>> Arguments: - * Matrix (char *) - * Pointer to the matrix. - * RHS (RealVector) - * RHS is the right hand side. This is what is being solved for. - * This is only the real portion of the right-hand side if the matrix - * is complex and spSEPARATED_COMPLEX_VECTORS is set true. - * Solution (RealVector) - * Solution is the vector being multiplied by the matrix. This is only - * the real portion if the matrix is complex and - * spSEPARATED_COMPLEX_VECTORS is set true. - * iRHS (RealVector) - * iRHS is the imaginary portion of the right hand side. This is - * what is being solved for. This is only necessary if the matrix is - * complex and spSEPARATED_COMPLEX_VECTORS is true. - * iSolution (RealVector) - * iSolution is the imaginary portion of the vector being multiplied - * by the matrix. This is only necessary if the matrix is - * complex and spSEPARATED_COMPLEX_VECTORS is true. - * - */ - -static void ComplexTransposedMatrixMultiply(MatrixPtr Matrix, RealVector RHS, RealVector Solution, std::optional iRHS, std::optional iSolution) -{ - ElementPtr pElement; - ComplexVector Vector; - ComplexNumber Sum; - int I, *pExtOrder; - -/* Begin `ComplexMatrixMultiply'. */ - -/* Correct array pointers for ARRAY_OFFSET. */ -#if NOT ARRAY_OFFSET -#if spSEPARATED_COMPLEX_VECTORS - --RHS; - --iRHS; - --Solution; - --iSolution; -#else - RHS -= 2; - Solution -= 2; -#endif -#endif - - /* Initialize Intermediate vector with reordered Solution vector. */ - Vector = (ComplexVector)Matrix->Intermediate; - pExtOrder = &Matrix->IntToExtRowMap[Matrix->Size]; - -#if spSEPARATED_COMPLEX_VECTORS - for (I = Matrix->Size; I > 0; I--) { - Vector[I].Real = Solution[*pExtOrder]; - ASSERT(iSolution.has_value()) - Vector[I].Imag = iSolution.value()[*(pExtOrder--)]; - } -#else - for (I = Matrix->Size; I > 0; I--) - Vector[I] = ((ComplexVector)Solution)[*(pExtOrder--)]; -#endif - - pExtOrder = &Matrix->IntToExtColMap[Matrix->Size]; - for (I = Matrix->Size; I > 0; I--) { - pElement = Matrix->FirstInCol[I]; - Sum.Real = Sum.Imag = 0.0; - - while (pElement != NULL) { /* Cmplx expression : Sum += Element * Vector[Row] */ - CMPLX_MULT_ADD_ASSIGN(Sum, *pElement, Vector[pElement->Row]); - pElement = pElement->NextInCol; - } - -#if spSEPARATED_COMPLEX_VECTORS - RHS[*pExtOrder] = Sum.Real; - ASSERT(iRHS != std::nullopt); - iRHS.value()[*pExtOrder--] = Sum.Imag; -#else - ((ComplexVector)RHS)[*pExtOrder--] = Sum; -#endif - } - return; -} -#endif /* spCOMPLEX AND MULTIPLICATION AND TRANSPOSE */ #if DETERMINANT /* @@ -859,8 +540,6 @@ static void ComplexTransposedMatrixMultiply(MatrixPtr Matrix, RealVector RHS, Re * number is scaled to be greater than or equal to 1.0 and less than 10.0. * * >>> Local variables: - * Norm (RealNumber) - * L-infinity norm of a complex number. * Size (int) * Local storage for Matrix->Size. * Temp (RealNumber) @@ -871,10 +550,6 @@ void spDeterminant(char* eMatrix, int* pExponent, RealNumber* pDeterminant, std: { MatrixPtr Matrix = (MatrixPtr)eMatrix; int I, Size; - RealNumber Norm, nr, ni; - ComplexNumber Pivot, cDeterminant; - -#define NORM(a) (nr = ABS((a).Real), ni = ABS((a).Imag), MAX(nr, ni)) /* Begin `spDeterminant'. */ ASSERT(IS_SPARSE(Matrix) AND IS_FACTORED(Matrix)); @@ -882,72 +557,12 @@ void spDeterminant(char* eMatrix, int* pExponent, RealNumber* pDeterminant, std: if (Matrix->Error == spSINGULAR) { *pDeterminant = 0.0; -#if spCOMPLEX - if (Matrix->Complex) { - ASSERT(piDeterminant != std::nullopt); - *piDeterminant.value() = 0.0; - } -#endif return; } Size = Matrix->Size; I = 0; -#if spCOMPLEX - if (Matrix->Complex) /* Complex Case. */ - { - cDeterminant.Real = 1.0; - cDeterminant.Imag = 0.0; - - while (++I <= Size) { - CMPLX_RECIPROCAL(Pivot, *Matrix->Diag[I]); - CMPLX_MULT_ASSIGN(cDeterminant, Pivot); - - /* Scale Determinant. */ - Norm = NORM(cDeterminant); - if (Norm != 0.0) { - while (Norm >= 1.0e12) { - cDeterminant.Real *= 1.0e-12; - cDeterminant.Imag *= 1.0e-12; - *pExponent += 12; - Norm = NORM(cDeterminant); - } - while (Norm < 1.0e-12) { - cDeterminant.Real *= 1.0e12; - cDeterminant.Imag *= 1.0e12; - *pExponent -= 12; - Norm = NORM(cDeterminant); - } - } - } - - /* Scale Determinant again, this time to be between 1.0 <= x < 10.0. */ - Norm = NORM(cDeterminant); - if (Norm != 0.0) { - while (Norm >= 10.0) { - cDeterminant.Real *= 0.1; - cDeterminant.Imag *= 0.1; - (*pExponent)++; - Norm = NORM(cDeterminant); - } - while (Norm < 1.0) { - cDeterminant.Real *= 10.0; - cDeterminant.Imag *= 10.0; - (*pExponent)--; - Norm = NORM(cDeterminant); - } - } - if (Matrix->NumberOfInterchangesIsOdd) - CMPLX_NEGATE(cDeterminant); - - *pDeterminant = cDeterminant.Real; - *piDeterminant.value() = cDeterminant.Imag; - } -#endif /* spCOMPLEX */ -#if REAL AND spCOMPLEX - else -#endif #if REAL { /* Real Case. */ *pDeterminant = 1.0; @@ -1074,123 +689,6 @@ void spStripFills(char* eMatrix) } #endif -#if TRANSLATE AND DELETE -/* - * DELETE A ROW AND COLUMN FROM THE MATRIX - * - * Deletes a row and a column from a matrix. - * - * Sparse will abort if an attempt is made to delete a row or column that - * doesn't exist. - * - * >>> Arguments: - * eMatrix (char *) - * Pointer to the matrix in which the row and column are to be deleted. - * Row (int) - * Row to be deleted. - * Col (int) - * Column to be deleted. - * - * >>> Local variables: - * ExtCol (int) - * The external column that is being deleted. - * ExtRow (int) - * The external row that is being deleted. - * pElement (ElementPtr) - * Pointer to an element in the matrix. Used when scanning rows and - * columns in order to eliminate elements from the last row or column. - * ppElement (ElementPtr *) - * Pointer to the location of an ElementPtr. This location will be - * filled with a NULL pointer if it is the new last element in its row - * or column. - * pElement (ElementPtr) - * Pointer to an element in the last row or column of the matrix. - * Size (int) - * The local version Matrix->Size, the size of the matrix. - */ - -void spDeleteRowAndCol(char* eMatrix, int Row, int Col) -{ - MatrixPtr Matrix = (MatrixPtr)eMatrix; - ElementPtr pElement, *ppElement, pLastElement; - int Size, ExtRow, ExtCol; - - /* Begin `spDeleteRowAndCol'. */ - - ASSERT(IS_SPARSE(Matrix) AND Row > 0 AND Col > 0); - ASSERT(Row <= Matrix->ExtSize AND Col <= Matrix->ExtSize); - - Size = Matrix->Size; - ExtRow = Row; - ExtCol = Col; - if (NOT Matrix->RowsLinked) - spcLinkRows(Matrix); - - Row = Matrix->ExtToIntRowMap[Row]; - Col = Matrix->ExtToIntColMap[Col]; - ASSERT(Row > 0 AND Col > 0); - - /* Move Row so that it is the last row in the matrix. */ - if (Row != Size) - spcRowExchange(Matrix, Row, Size); - - /* Move Col so that it is the last column in the matrix. */ - if (Col != Size) - spcColExchange(Matrix, Col, Size); - - /* Correct Diag pointers. */ - if (Row == Col) - SWAP(ElementPtr, Matrix->Diag[Row], Matrix->Diag[Size]) - else { - Matrix->Diag[Row] = spcFindElementInCol(Matrix, Matrix->FirstInCol + Row, - Row, Row, NO); - Matrix->Diag[Col] = spcFindElementInCol(Matrix, Matrix->FirstInCol + Col, - Col, Col, NO); - } - - /* - * Delete last row and column of the matrix. - */ - /* Break the column links to every element in the last row. */ - pLastElement = Matrix->FirstInRow[Size]; - while (pLastElement != NULL) { - ppElement = &(Matrix->FirstInCol[pLastElement->Col]); - while ((pElement = *ppElement) != NULL) { - if (pElement == pLastElement) - *ppElement = NULL; /* Unlink last element in column. */ - else - ppElement = &pElement->NextInCol; /* Skip element. */ - } - pLastElement = pLastElement->NextInRow; - } - - /* Break the row links to every element in the last column. */ - pLastElement = Matrix->FirstInCol[Size]; - while (pLastElement != NULL) { - ppElement = &(Matrix->FirstInRow[pLastElement->Row]); - while ((pElement = *ppElement) != NULL) { - if (pElement == pLastElement) - *ppElement = NULL; /* Unlink last element in row. */ - else - ppElement = &pElement->NextInRow; /* Skip element. */ - } - pLastElement = pLastElement->NextInCol; - } - - /* Clean up some details. */ - Matrix->Size = Size - 1; - Matrix->Diag[Size] = NULL; - Matrix->FirstInRow[Size] = NULL; - Matrix->FirstInCol[Size] = NULL; - Matrix->CurrentSize--; - Matrix->ExtToIntRowMap[ExtRow] = -1; - Matrix->ExtToIntColMap[ExtCol] = -1; - Matrix->NeedsOrdering = YES; - - return; -} -#endif - #if PSEUDOCONDITION /* * CALCULATE PSEUDOCONDITION @@ -1321,23 +819,14 @@ RealNumber spCondition(char* eMatrix, RealNumber NormOfMatrix, int* pError) return 0.0; } -#if spCOMPLEX - if (Matrix->Complex) - return ComplexCondition(Matrix, NormOfMatrix, pError); -#endif - #if REAL Size = Matrix->Size; T = Matrix->Intermediate; -#if spCOMPLEX - Tm = Matrix->Intermediate + Size; -#else Tm = ALLOC(RealNumber, Size + 1); if (Tm == NULL) { *pError = spNO_MEMORY; return 0.0; } -#endif for (I = Size; I > 0; I--) T[I] = 0.0; @@ -1484,9 +973,7 @@ RealNumber spCondition(char* eMatrix, RealNumber NormOfMatrix, int* pError) for (ASz = 0.0, I = Size; I > 0; I--) ASz += ABS(T[I]); -#if NOT spCOMPLEX FREE(Tm); -#endif Linpack = ASy / ASz; OLeary = E / MaxY; @@ -1495,210 +982,6 @@ RealNumber spCondition(char* eMatrix, RealNumber NormOfMatrix, int* pError) #endif /* REAL */ } -#if spCOMPLEX -/* - * ESTIMATE CONDITION NUMBER - * - * Complex version of spCondition(). - * - * >>> Returns: - * The reciprocal of the condition number. - * - * >>> Arguments: - * Matrix (MatrixPtr) - * Pointer to the matrix. - * NormOfMatrix (RealNumber) - * The L-infinity norm of the unfactored matrix as computed by - * spNorm(). - * pError (int *) - * Used to return error code. - * - * >>> Possible errors: - * spNO_MEMORY - */ - -static RealNumber ComplexCondition(MatrixPtr Matrix, RealNumber NormOfMatrix, int* pError) -{ - ElementPtr pElement; - ComplexVector T, Tm; - int I, K, Row; - ElementPtr pPivot; - int Size; - RealNumber E, Em, ASp, ASm, ASw, ASy, ASv, ASz, MaxY, ScaleFactor; - RealNumber Linpack, OLeary, InvNormOfInverse; - ComplexNumber Wp, Wm; - - /* Begin `ComplexCondition'. */ - - Size = Matrix->Size; - T = (ComplexVector)Matrix->Intermediate; - Tm = ALLOC(ComplexNumber, Size + 1); - if (Tm == NULL) { - *pError = spNO_MEMORY; - return 0.0; - } - for (I = Size; I > 0; I--) - T[I].Real = T[I].Imag = 0.0; - - /* - * Part 1. Ay = e. - * Solve Ay = LUy = e where e consists of +1 and -1 terms with the sign - * chosen to maximize the size of w in Lw = e. Since the terms in w can - * get very large, scaling is used to avoid overflow. - */ - - /* Forward elimination. Solves Lw = e while choosing e. */ - E = 1.0; - for (I = 1; I <= Size; I++) { - pPivot = Matrix->Diag[I]; - if (T[I].Real < 0.0) - Em = -E; - else - Em = E; - Wm = T[I]; - Wm.Real += Em; - ASm = CMPLX_1_NORM(Wm); - CMPLX_MULT_ASSIGN(Wm, *pPivot); - if (CMPLX_1_NORM(Wm) > SLACK) { - ScaleFactor = 1.0 / MAX(SQR(SLACK), CMPLX_1_NORM(Wm)); - for (K = Size; K > 0; K--) - SCLR_MULT_ASSIGN(T[K], ScaleFactor); - E *= ScaleFactor; - Em *= ScaleFactor; - ASm *= ScaleFactor; - SCLR_MULT_ASSIGN(Wm, ScaleFactor); - } - Wp = T[I]; - Wp.Real -= Em; - ASp = CMPLX_1_NORM(Wp); - CMPLX_MULT_ASSIGN(Wp, *pPivot); - - /* Update T for both values of W, minus value is placed in Tm. */ - pElement = pPivot->NextInCol; - while (pElement != NULL) { - Row = pElement->Row; - /* Cmplx expr: Tm[Row] = T[Row] - (Wp * *pElement). */ - CMPLX_MULT_SUBT(Tm[Row], Wm, *pElement, T[Row]); - /* Cmplx expr: T[Row] -= Wp * *pElement. */ - CMPLX_MULT_SUBT_ASSIGN(T[Row], Wm, *pElement); - ASp += CMPLX_1_NORM(T[Row]); - ASm += CMPLX_1_NORM(Tm[Row]); - pElement = pElement->NextInCol; - } - - /* If minus value causes more growth, overwrite T with its values. */ - if (ASm > ASp) { - T[I] = Wm; - pElement = pPivot->NextInCol; - while (pElement != NULL) { - T[pElement->Row] = Tm[pElement->Row]; - pElement = pElement->NextInCol; - } - } else - T[I] = Wp; - } - - /* Compute 1-norm of T, which now contains w, and scale ||T|| to 1/SLACK. */ - for (ASw = 0.0, I = Size; I > 0; I--) - ASw += CMPLX_1_NORM(T[I]); - ScaleFactor = 1.0 / (SLACK * ASw); - if (ScaleFactor < 0.5) { - for (I = Size; I > 0; I--) - SCLR_MULT_ASSIGN(T[I], ScaleFactor); - E *= ScaleFactor; - } - - /* Backward Substitution. Solves Uy = w.*/ - for (I = Size; I >= 1; I--) { - pElement = Matrix->Diag[I]->NextInRow; - while (pElement != NULL) { /* Cmplx expr: T[I] -= T[pElement->Col] * *pElement. */ - CMPLX_MULT_SUBT_ASSIGN(T[I], T[pElement->Col], *pElement); - pElement = pElement->NextInRow; - } - if (CMPLX_1_NORM(T[I]) > SLACK) { - ScaleFactor = 1.0 / MAX(SQR(SLACK), CMPLX_1_NORM(T[I])); - for (K = Size; K > 0; K--) - SCLR_MULT_ASSIGN(T[K], ScaleFactor); - E *= ScaleFactor; - } - } - - /* Compute 1-norm of T, which now contains y, and scale ||T|| to 1/SLACK. */ - for (ASy = 0.0, I = Size; I > 0; I--) - ASy += CMPLX_1_NORM(T[I]); - ScaleFactor = 1.0 / (SLACK * ASy); - if (ScaleFactor < 0.5) { - for (I = Size; I > 0; I--) - SCLR_MULT_ASSIGN(T[I], ScaleFactor); - ASy = 1.0 / SLACK; - E *= ScaleFactor; - } - - /* Compute infinity-norm of T for O'Leary's estimate. */ - for (MaxY = 0.0, I = Size; I > 0; I--) - if (MaxY < CMPLX_1_NORM(T[I])) - MaxY = CMPLX_1_NORM(T[I]); - - /* - * Part 2. A* z = y where the * represents the transpose. - * Recall that A = LU implies A* = U* L*. - */ - - /* Forward elimination, U* v = y. */ - for (I = 1; I <= Size; I++) { - pElement = Matrix->Diag[I]->NextInRow; - while (pElement != NULL) { /* Cmplx expr: T[pElement->Col] -= T[I] * *pElement. */ - CMPLX_MULT_SUBT_ASSIGN(T[pElement->Col], T[I], *pElement); - pElement = pElement->NextInRow; - } - if (CMPLX_1_NORM(T[I]) > SLACK) { - ScaleFactor = 1.0 / MAX(SQR(SLACK), CMPLX_1_NORM(T[I])); - for (K = Size; K > 0; K--) - SCLR_MULT_ASSIGN(T[K], ScaleFactor); - ASy *= ScaleFactor; - } - } - - /* Compute 1-norm of T, which now contains v, and scale ||T|| to 1/SLACK. */ - for (ASv = 0.0, I = Size; I > 0; I--) - ASv += CMPLX_1_NORM(T[I]); - ScaleFactor = 1.0 / (SLACK * ASv); - if (ScaleFactor < 0.5) { - for (I = Size; I > 0; I--) - SCLR_MULT_ASSIGN(T[I], ScaleFactor); - ASy *= ScaleFactor; - } - - /* Backward Substitution, L* z = v. */ - for (I = Size; I >= 1; I--) { - pPivot = Matrix->Diag[I]; - pElement = pPivot->NextInCol; - while (pElement != NULL) { /* Cmplx expr: T[I] -= T[pElement->Row] * *pElement. */ - CMPLX_MULT_SUBT_ASSIGN(T[I], T[pElement->Row], *pElement); - pElement = pElement->NextInCol; - } - CMPLX_MULT_ASSIGN(T[I], *pPivot); - if (CMPLX_1_NORM(T[I]) > SLACK) { - ScaleFactor = 1.0 / MAX(SQR(SLACK), CMPLX_1_NORM(T[I])); - for (K = Size; K > 0; K--) - SCLR_MULT_ASSIGN(T[K], ScaleFactor); - ASy *= ScaleFactor; - } - } - - /* Compute 1-norm of T, which now contains z. */ - for (ASz = 0.0, I = Size; I > 0; I--) - ASz += CMPLX_1_NORM(T[I]); - - FREE(Tm); - - Linpack = ASy / ASz; - OLeary = E / MaxY; - InvNormOfInverse = MIN(Linpack, OLeary); - return InvNormOfInverse / NormOfMatrix; -} -#endif /* spCOMPLEX */ - /* * L-INFINITY MATRIX NORM * @@ -1727,33 +1010,17 @@ RealNumber spNorm(char* eMatrix) if (NOT Matrix->RowsLinked) spcLinkRows(Matrix); -/* Compute row sums. */ + /* Compute row sums. */ #if REAL - if (NOT Matrix->Complex) { - for (I = Matrix->Size; I > 0; I--) { - pElement = Matrix->FirstInRow[I]; - AbsRowSum = 0.0; - while (pElement != NULL) { - AbsRowSum += ABS(pElement->Real); - pElement = pElement->NextInRow; - } - if (Max < AbsRowSum) - Max = AbsRowSum; - } - } -#endif -#if spCOMPLEX - if (Matrix->Complex) { - for (I = Matrix->Size; I > 0; I--) { - pElement = Matrix->FirstInRow[I]; - AbsRowSum = 0.0; - while (pElement != NULL) { - AbsRowSum += CMPLX_1_NORM(*pElement); - pElement = pElement->NextInRow; - } - if (Max < AbsRowSum) - Max = AbsRowSum; + for (I = Matrix->Size; I > 0; I--) { + pElement = Matrix->FirstInRow[I]; + AbsRowSum = 0.0; + while (pElement != NULL) { + AbsRowSum += ABS(pElement->Real); + pElement = pElement->NextInRow; } + if (Max < AbsRowSum) + Max = AbsRowSum; } #endif return Max; @@ -1833,14 +1100,13 @@ RealNumber spLargestElement(char* eMatrix) int I; RealNumber Mag, AbsColSum, Max = 0.0, MaxRow = 0.0, MaxCol = 0.0; RealNumber Pivot; - ComplexNumber cPivot; ElementPtr pElement, pDiag; /* Begin `spLargestElement'. */ ASSERT(IS_SPARSE(Matrix)); #if REAL - if (Matrix->Factored AND NOT Matrix->Complex) { + if (Matrix->Factored) { if (Matrix->Error == spSINGULAR) return 0.0; @@ -1871,7 +1137,7 @@ RealNumber spLargestElement(char* eMatrix) if (AbsColSum > MaxCol) MaxCol = AbsColSum; } - } else if (NOT Matrix->Complex) { + } else { for (I = 1; I <= Matrix->Size; I++) { pElement = Matrix->FirstInCol[I]; while (pElement != NULL) { @@ -1883,51 +1149,6 @@ RealNumber spLargestElement(char* eMatrix) } return Max; } -#endif -#if spCOMPLEX - if (Matrix->Factored AND Matrix->Complex) { - if (Matrix->Error == spSINGULAR) - return 0.0; - - /* Find the bound on the size of the largest element over all factorization. */ - for (I = 1; I <= Matrix->Size; I++) { - pDiag = Matrix->Diag[I]; - - /* Lower triangular matrix. */ - CMPLX_RECIPROCAL(cPivot, *pDiag); - Mag = CMPLX_1_NORM(cPivot); - if (Mag > MaxRow) - MaxRow = Mag; - pElement = Matrix->FirstInRow[I]; - while (pElement != pDiag) { - Mag = CMPLX_1_NORM(*pElement); - if (Mag > MaxRow) - MaxRow = Mag; - pElement = pElement->NextInRow; - } - - /* Upper triangular matrix. */ - pElement = Matrix->FirstInCol[I]; - AbsColSum = 1.0; /* Diagonal of U is unity. */ - while (pElement != pDiag) { - AbsColSum += CMPLX_1_NORM(*pElement); - pElement = pElement->NextInCol; - } - if (AbsColSum > MaxCol) - MaxCol = AbsColSum; - } - } else if (Matrix->Complex) { - for (I = 1; I <= Matrix->Size; I++) { - pElement = Matrix->FirstInCol[I]; - while (pElement != NULL) { - Mag = CMPLX_1_NORM(*pElement); - if (Mag > Max) - Max = Mag; - pElement = pElement->NextInCol; - } - } - return Max; - } #endif return MaxRow * MaxCol; } From 7f8ba2d4a3389ac69a54f308bd6b090ddea24915 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Mon, 21 Oct 2024 10:54:26 +0200 Subject: [PATCH 08/17] We should not used save-always anymore (#3136) Here is the depreciation: https://github.com/actions/cache/pull/1452 --- .github/workflows/neuron-ci.yml | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/.github/workflows/neuron-ci.yml b/.github/workflows/neuron-ci.yml index dfe8a9327c..c77d5a6b45 100644 --- a/.github/workflows/neuron-ci.yml +++ b/.github/workflows/neuron-ci.yml @@ -237,10 +237,10 @@ jobs: echo ----- - name: Restore compiler cache - uses: actions/cache@v4 + uses: actions/cache/restore@v4 + id: restore-compiler-cache with: path: ${{runner.workspace}}/ccache - save-always: true key: ${{matrix.os}}-${{hashfiles('matrix.json')}}-${{github.ref}}-${{github.sha}} restore-keys: | ${{matrix.os}}-${{hashfiles('matrix.json')}}-${{github.ref}}- @@ -441,6 +441,15 @@ jobs: INSTALL_DIR : ${{ runner.workspace }}/install MATRIX_EVAL: ${{ matrix.config.matrix_eval }} + - name: Save compiler cache + uses: actions/cache/save@v4 + if: always() && steps.restore-compiler-cache.outputs.cache-hit != 'true' + with: + path: ${{runner.workspace}}/ccache + key: | + ${{matrix.os}}-${{hashfiles('matrix.json')}}-${{github.ref}}- + ${{matrix.os}}-${{hashfiles('matrix.json')}}- + # This step will set up an SSH connection on tmate.io for live debugging. # To enable it, you have to: # * add 'live-debug-ci' to your PR title From ceedb15fc50214aed566af5ded5382fe835b6d9a Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Mon, 21 Oct 2024 13:37:52 +0200 Subject: [PATCH 09/17] nrn_py_CallObject take a nb::object for args (#3127) This object can be any iterable like list or tuple. --- src/nrnpython/nrnpy_hoc.cpp | 1 - src/nrnpython/nrnpy_p2h.cpp | 21 +++++++++++---------- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/nrnpython/nrnpy_hoc.cpp b/src/nrnpython/nrnpy_hoc.cpp index 311139567c..96ae36c38b 100644 --- a/src/nrnpython/nrnpy_hoc.cpp +++ b/src/nrnpython/nrnpy_hoc.cpp @@ -21,7 +21,6 @@ #include #include - namespace nb = nanobind; extern PyTypeObject* psection_type; diff --git a/src/nrnpython/nrnpy_p2h.cpp b/src/nrnpython/nrnpy_p2h.cpp index fc0a4d085b..6bfae26924 100644 --- a/src/nrnpython/nrnpy_p2h.cpp +++ b/src/nrnpython/nrnpy_p2h.cpp @@ -18,7 +18,7 @@ namespace nb = nanobind; static char* nrnpyerr_str(); -static nb::object nrnpy_pyCallObject(nb::callable, PyObject*); +static nb::object nrnpy_pyCallObject(nb::callable, nb::object); static PyObject* main_module; static PyObject* main_namespace; @@ -46,7 +46,7 @@ static void call_python_with_section(Object* pyact, Section* sec) { nanobind::gil_scoped_acquire lock{}; nb::tuple args = nb::make_tuple(reinterpret_cast(newpysechelp(sec))); - nb::object r = nrnpy_pyCallObject(po, args.ptr()); + nb::object r = nrnpy_pyCallObject(po, args); if (!r.is_valid()) { char* mes = nrnpyerr_str(); if (mes) { @@ -106,12 +106,13 @@ Object* nrnpy_pyobject_in_obj(PyObject* po) { return on; } -static nb::object nrnpy_pyCallObject(nb::callable callable, PyObject* args) { +static nb::object nrnpy_pyCallObject(nb::callable callable, nb::object args) { // When hoc calls a PythonObject method, then in case python // calls something back in hoc, the hoc interpreter must be // at the top level HocTopContextSet - PyObject* p = PyObject_CallObject(callable.ptr(), args); + nb::tuple tup(args); + nb::object p = nb::steal(PyObject_CallObject(callable.ptr(), tup.ptr())); #if 0 printf("PyObject_CallObject callable\n"); PyObject_Print(callable, stdout, 0); @@ -139,7 +140,7 @@ printf("\nreturn %p\n", p); } } **/ - return nb::steal(p); + return p; } static void py2n_component(Object* ob, Symbol* sym, int nindex, int isfunc) { @@ -196,7 +197,7 @@ static void py2n_component(Object* ob, Symbol* sym, int nindex, int isfunc) { } } // printf("PyObject_CallObject %s %p\n", sym->name, tail); - result = nrnpy_pyCallObject(nb::borrow(tail), args).release().ptr(); + result = nrnpy_pyCallObject(nb::borrow(tail), nb::borrow(args)).release().ptr(); Py_DECREF(args); // PyObject_Print(result, stdout, 0); // printf(" result of call\n"); @@ -333,9 +334,9 @@ static nb::object hoccommand_exec_help1(nb::object po) { if (!nb::tuple::check_(args)) { args = nb::make_tuple(args); } - return nrnpy_pyCallObject(po[0], args.ptr()); + return nrnpy_pyCallObject(po[0], args); } else { - return nrnpy_pyCallObject(nb::borrow(po), nb::tuple().ptr()); + return nrnpy_pyCallObject(nb::borrow(po), nb::tuple()); } } @@ -417,7 +418,7 @@ static void grphcmdtool(Object* ho, int type, double x, double y, int key) { nanobind::gil_scoped_acquire lock{}; nb::tuple args = nb::make_tuple(type, x, y, key); - nb::object r = nrnpy_pyCallObject(po, args.ptr()); + nb::object r = nrnpy_pyCallObject(po, args); if (!r.is_valid()) { char* mes = nrnpyerr_str(); if (mes) { @@ -483,7 +484,7 @@ static double func_call(Object* ho, int narg, int* err) { } } - nb::object r = nrnpy_pyCallObject(po, args); + nb::object r = nrnpy_pyCallObject(po, nb::borrow(args)); Py_XDECREF(args); double rval = 0.0; if (!r.is_valid()) { From 510b9bc94592e77d282d28921db483e25a1400f6 Mon Sep 17 00:00:00 2001 From: nrnhines Date: Mon, 21 Oct 2024 20:50:56 -0400 Subject: [PATCH 10/17] NRN_ENABLE_DIGEST and NRN_ENABLE_ARCH_INDEP_EXP_POW (#3135) * Selected messages from old digest-debug branch Print digest of cvode f(y,t) and solvex(b,...) cmake -DNRN_DIGEST=ON nrn_digest() starts the accumulation of cvode digest info. nrn_digest("filename") prints accumulated cvode digest info. filename format is: message threadid index t digest where digest is the first 16 hex characters of the SHA1 hash of the double* array indicated by the message. -DNRN_ENABLE_ARCH_INDEP_EXP_POW=ON (default OFF) Provides h.use_exp_pow_precision(style) style = 0 means use normal machine IEEE precision for exp(x) and pow(x,y) style = 1 means use 53 bit mpfr style = 2 means use IEEE but truncate to 32 bit precision. sundials uses hoc_pow. digest format has a f(y,t) call count. * nrn_digest(tid, i) will print the details of the i'th digest item of thread tid. * sundial RPowerR uses hoc_pow. nrn_digest(tid, i, "abort") calls abort() on reaching index i of thread tid * Documentation for nrn_digest and use_exp_pow_precision --- CMakeLists.txt | 11 +++ cmake/BuildOptionDefaults.cmake | 2 + cmake_nrnconf.h.in | 7 ++ docs/cmake_doc/options.rst | 20 +++++ docs/python/programming/internals.rst | 112 +++++++++++++++++++++++++ src/nocmodl/nocpout.cpp | 2 + src/nrncvode/occvode.cpp | 31 +++++-- src/nrniv/CMakeLists.txt | 20 +++++ src/nrniv/nmodlrandom.cpp | 8 +- src/nrniv/nrnmenu.cpp | 3 +- src/nrnpython/CMakeLists.txt | 3 +- src/oc/debug.cpp | 92 +++++++++++++++++++++ src/oc/hoc_init.cpp | 4 + src/oc/math.cpp | 114 +++++++++++++++++++++++++- src/oc/nrndigest.h | 19 +++++ src/oc/oc_ansi.h | 1 + src/oc/ocfunc.h | 6 ++ src/sundials/shared/sundialsmath.c | 6 +- 18 files changed, 444 insertions(+), 17 deletions(-) create mode 100644 src/oc/nrndigest.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 87d423a7b4..57bd453036 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -85,6 +85,13 @@ option( NRN_ENABLE_PERFORMANCE_TESTS "Enable tests that measure performance. These are known to be unreliable when run on busy/oversubscribed machines such as CI runners." ${NRN_ENABLE_PERFORMANCE_TESTS_DEFAULT}) +option(NRN_ENABLE_DIGEST + "Provides nrn_digest function for debugging cross platform floating result differences." + ${NRN_ENABLE_DIGEST_DEFAULT}) +option( + NRN_ENABLE_ARCH_INDEP_EXP_POW + "Provides use_exp_pow_precision(style) function so that exp and pow produce same results on all platforms" + ${NRN_ENABLE_ARCH_INDEP_EXP_POW_DEFAULT}) # This can be helpful in very specific CI build configurations, where ccache is used *and* different # CI builds are built under different directories. option(NRN_AVOID_ABSOLUTE_PATHS @@ -1056,6 +1063,10 @@ if(NRN_ENABLE_PROFILING) message(STATUS " Caliper | ${caliper_DIR}") endif() endif() +if(NRN_ENABLE_DIGEST OR NRN_ARCH_INDEP_EXP_POW) + message(STATUS "NRN_ENABLE_DIGEST | ${NRN_ENABLE_DIGEST}") + message(STATUS "NRN_ENABLE_ARCH_INDEP_EXP_POW | ${NRN_ENABLE_ARCH_INDEP_EXP_POW}") +endif() message(STATUS "--------------+--------------------------------------------------------------") message(STATUS " See documentation : https://www.neuron.yale.edu/neuron/") message(STATUS "--------------+--------------------------------------------------------------") diff --git a/cmake/BuildOptionDefaults.cmake b/cmake/BuildOptionDefaults.cmake index e3e1f3dd63..b43338e5ac 100644 --- a/cmake/BuildOptionDefaults.cmake +++ b/cmake/BuildOptionDefaults.cmake @@ -29,6 +29,8 @@ set(NRN_AVOID_ABSOLUTE_PATHS_DEFAULT OFF) set(NRN_NMODL_CXX_FLAGS_DEFAULT "-O0") set(NRN_SANITIZERS_DEFAULT "") set(NRN_ENABLE_MATH_OPT_DEFAULT OFF) +set(NRN_ENABLE_DIGEST_DEFAULT OFF) +set(NRN_ENABLE_ARCH_INDEP_EXP_POW_DEFAULT OFF) # Some distributions may set the prefix. To avoid errors, unset it set(NRN_PYTHON_DYNAMIC_DEFAULT "") diff --git a/cmake_nrnconf.h.in b/cmake_nrnconf.h.in index 5bbd435916..3c7851c6d5 100644 --- a/cmake_nrnconf.h.in +++ b/cmake_nrnconf.h.in @@ -1,5 +1,12 @@ #pragma once +/* Define to one if want to debug using sha1 hashes of data */ +#cmakedefine01 NRN_ENABLE_DIGEST + +/* Define to one if want to allow selection of architecture independent */ +/* 53 bit double precision of exp and pow from mpfr */ +#cmakedefine01 NRN_ENABLE_ARCH_INDEP_EXP_POW + /* Define if building universal (internal helper macro) */ #cmakedefine AC_APPLE_UNIVERSAL_BUILD @AC_APPLE_UNIVERSAL_BUILD@ diff --git a/docs/cmake_doc/options.rst b/docs/cmake_doc/options.rst index 087b477ab1..fcb296682a 100644 --- a/docs/cmake_doc/options.rst +++ b/docs/cmake_doc/options.rst @@ -674,3 +674,23 @@ NRN_ENABLE_MATH_OPT:BOOL=OFF Note: Compilers like Intel, NVHPC, Cray etc enable such optimisations by default. + +NRN_ENABLE_DIGEST:BOOL=OFF +------------------------------ + Provides \ :func:`nrn_digest` function for debugging cross platform floating + result differences. + + Requires libcrypto + +NRN_ENABLE_ARCH_INDEP_EXP_POW:BOOL=OFF +--------------------------------- + Provides \ :func:`use_exp_pow_precision` function so that exp and pow produce + same results on all platforms. + + Requires mpfr (multiple precision floating-point computation). eg. + ``sudo apt install libmpfr-dev`` + + To get platform independent floating point results with clang, + also consider using + ``-DCMAKE_C_FLAGS="-ffp-contract=off" -DCMAKE_CXX_FLAGS="-ffp-contract=off"`` + or, alternatively, ``"-fp-model=strict`` diff --git a/docs/python/programming/internals.rst b/docs/python/programming/internals.rst index fab28b96e0..5f6096bae6 100755 --- a/docs/python/programming/internals.rst +++ b/docs/python/programming/internals.rst @@ -261,3 +261,115 @@ Miscellaneous the variable from its interpreter name. Not needed by or useful for the user; returns 1.0 on success. +---- + +Debugging +~~~~~~~~~~~ + +.. function:: nrn_digest + + Syntax: + ``h.nrn_digest()`` + + ``h.nrn_digest(tid, i)`` + + ``h.nrn_digest(tid, i, "abort")`` + + ``h.nrn_digest(filename)`` + + Description: + Available when configured with the cmake option ``-DNRN_ENABLE_DIGEST=ON`` + + If the same simulation gives different results on different machines, + this function can help isolate the statement that generates the + first difference during the simulation. + I think :meth:`ParallelContext.prcellstate` is generally better, but in rare + situations, nrn_digest can be very helpful. + + The first three forms begin digest gathering. The last form + prints the gathered digest information to the filename. + With just the two ``tid, i`` arguments, the i gathered item of the + tid thread is printed (for single thread simulations, use ``tid = 0``), + to the terminal as well as the individual values of the array + for that digest item. With the third ``"abort"`` argument, the + ith gathered item is printed and ``abort()`` is called (dropping + into gdb if that is being used so that one can observe the backtrace). + + Lines are inserted into the digest by calling the C function declared + in ``src/oc/nrndigest.h``. + ``void nrn_digest_dbl_array(const char* msg, int tid, double t, double* array, size_t sz);`` + at the moment, such lines are present in ``src/nrncvode/occvode.cpp`` + to instrument the cvode callbacks that compute ``y' = f(y, t)`` and the + approximate jacobian matrix solver ``M*x = b``. I.e in part + + .. code-block:: + + #include "nrndigest.h" + ... + void Cvode::fun_thread(neuron::model_sorted_token const& sorted_token, + double tt, + double* y, + double* ydot, + NrnThread* nt) { + CvodeThreadData& z = CTD(nt->id); + #if NRN_DIGEST + if (nrn_digest_) { + nrn_digest_dbl_array("y", nt->id, tt, y, z.nvsize_); + } + #endif + ... + #if NRN_DIGEST + if (nrn_digest_ && ydot) { + nrn_digest_dbl_array("ydot", nt->id, tt, ydot, z.nvsize_); + } + #endif + + Note: when manually adding such lines, the conditional compilation and + nrn\_digest\_ test are not needed. The arguments to + ``nrn_digest_dbl_array`` determine the line added to the digest. + The 5th arg is the size of the 4th arg double array. The double array + is processed by SHA1 and the first 16 hex digits are appended to the line. + An example of the first few lines of output in a digest file is + .. code-block:: + + tid=0 size=1344 + y 0 0 0 e1f6a372856b45e6 + y 0 1 0 e1f6a372856b45e6 + ydot 0 2 0 523c9694c335e458 + y 0 3 4.7121609153871379e-09 fabb4bc469447404 + ydot 0 4 4.7121609153871379e-09 60bcff174645fc29 + + The first line is thread id and number of lines for that thread. + Other thread groups, if any, follow the end of each thread group. + The digest lines consist of thread id, line identifier (start from 0 + for each group), double value of the 3rd arg, hash of the array. + +---- + +.. function:: use_exp_pow_precision + + Syntax: + ``h.use_exp_pow_precision(istyle)`` + + Description: + Works when configured with the cmake option + ``-DNRN_ENABLE_ARCH_INDEP_EXP_POW=ON`` and otherwise does nothing. + + * istyle = 1 + All calls to :func:`exp` and :func:`pow` as well as their use + internally, in mod files, and by cvode, are computed on mac, linux, + windows so that double precision floating point results are + cross platform consistent. (Makes use of a + multiple precision floating-point computation library.) + + * istyle = 2 + exp and pow are rounded to 32 bits of mantissa + + * istyle = 0 + Default. + exp and pow calcualted natively (cross platform values can have + round off error differences.) + + When using clang (eg. on a mac) cross platform floating point + identity is often attainable with C and C++ flag option + ``"-ffp-contract=off"``. diff --git a/src/nocmodl/nocpout.cpp b/src/nocmodl/nocpout.cpp index 92c18d6b2d..c11d63e24d 100644 --- a/src/nocmodl/nocpout.cpp +++ b/src/nocmodl/nocpout.cpp @@ -256,6 +256,8 @@ void parout() { \n#if !NRNGPU\ \n#undef exp\ \n#define exp hoc_Exp\ +\n#undef pow\ +\n#define pow hoc_pow\ \n#endif\n\ "); if (protect_include_) { diff --git a/src/nrncvode/occvode.cpp b/src/nrncvode/occvode.cpp index 67e24d1c87..a2f6bd24eb 100644 --- a/src/nrncvode/occvode.cpp +++ b/src/nrncvode/occvode.cpp @@ -12,6 +12,7 @@ #include "vrecitem.h" #include "membfunc.h" #include "nonvintblock.h" +#include "nrndigest.h" #include #include @@ -323,7 +324,7 @@ void Cvode::new_no_cap_memb(CvodeThreadData& z, NrnThread* _nt) { } } } - assert(ncm->ml.size() == n); + assert(ncm->ml.size() == std::size_t(n)); } } @@ -456,7 +457,7 @@ extern void nrn_extra_scatter_gather(int, int); void Cvode::scatter_y(neuron::model_sorted_token const& sorted_token, double* y, int tid) { CvodeThreadData& z = CTD(tid); - assert(z.nonvint_extra_offset_ == z.pv_.size()); + assert(std::size_t(z.nonvint_extra_offset_) == z.pv_.size()); for (int i = 0; i < z.nonvint_extra_offset_; ++i) { // TODO: understand why this wasn't needed before if (z.pv_[i]) { @@ -494,7 +495,7 @@ void Cvode::gather_y(N_Vector y) { void Cvode::gather_y(double* y, int tid) { CvodeThreadData& z = CTD(tid); nrn_extra_scatter_gather(1, tid); - assert(z.nonvint_extra_offset_ == z.pv_.size()); + assert(std::size_t(z.nonvint_extra_offset_) == z.pv_.size()); for (int i = 0; i < z.nonvint_extra_offset_; ++i) { // TODO: understand why this wasn't needed before if (z.pv_[i]) { @@ -565,6 +566,12 @@ int Cvode::solvex_thread(neuron::model_sorted_token const& sorted_token, if (z.nvsize_ == 0) { return 0; } +#if NRN_DIGEST + if (nrn_digest_) { + nrn_digest_dbl_array("solvex enter b", nt->id, t_, b, z.nvsize_); + nrn_digest_dbl_array("solvex enter y", nt->id, t_, y, z.nvsize_); + } +#endif lhs(sorted_token, nt); // special version for cvode. scatter_ydot(b, nt->id); if (z.cmlcap_) { @@ -597,6 +604,11 @@ int Cvode::solvex_thread(neuron::model_sorted_token const& sorted_token, // printf("\texit b\n"); // for (i=0; i < neq_; ++i) { printf("\t\t%d %g\n", i, b[i]);} nrn_nonvint_block_ode_solve(z.nvsize_, b, y, nt->id); +#if NRN_DIGEST + if (nrn_digest_) { + nrn_digest_dbl_array("solvex leave b", nt->id, t_, b, z.nvsize_); + } +#endif return 0; } @@ -670,9 +682,20 @@ void Cvode::fun_thread(neuron::model_sorted_token const& sorted_token, double* ydot, NrnThread* nt) { CvodeThreadData& z = CTD(nt->id); +#if NRN_DIGEST + if (nrn_digest_) { + nrn_digest_dbl_array("y", nt->id, tt, y, z.nvsize_); + } +#endif fun_thread_transfer_part1(sorted_token, tt, y, nt); nrn_nonvint_block_ode_fun(z.nvsize_, y, ydot, nt->id); fun_thread_transfer_part2(sorted_token, ydot, nt); + +#if NRN_DIGEST + if (nrn_digest_ && ydot) { + nrn_digest_dbl_array("ydot", nt->id, tt, ydot, z.nvsize_); + } +#endif } void Cvode::fun_thread_transfer_part1(neuron::model_sorted_token const& sorted_token, @@ -740,7 +763,6 @@ void Cvode::fun_thread_transfer_part2(neuron::model_sorted_token const& sorted_t } void Cvode::fun_thread_ms_part1(double tt, double* y, NrnThread* nt) { - CvodeThreadData& z = ctd_[nt->id]; nt->_t = tt; // fix this!!! @@ -1001,7 +1023,6 @@ void Cvode::error_weights(double* pd) { void Cvode::acor(double* pd) { int i, id; - NrnThread* nt; for (id = 0; id < nctd_; ++id) { CvodeThreadData& z = ctd_[id]; double* s = n_vector_data(acorvec(), id); diff --git a/src/nrniv/CMakeLists.txt b/src/nrniv/CMakeLists.txt index f1f5830dde..0f9cd2564a 100644 --- a/src/nrniv/CMakeLists.txt +++ b/src/nrniv/CMakeLists.txt @@ -414,6 +414,26 @@ endif() if(NRN_ENABLE_THREADS) target_link_libraries(nrniv_lib Threads::Threads) endif() + +if(${NRN_ENABLE_DIGEST}) + if(NRN_MACOS_BUILD) + # where to get openssl/sha.h after brew install openssl + set_property( + SOURCE ${NRN_OC_SRC_DIR}/debug.cpp + APPEND + PROPERTY INCLUDE_DIRECTORIES /opt/homebrew/Cellar/openssl@3/3.1.0/include) + find_library(LIB_CRYPTO crypto PATHS /opt/homebrew/Cellar/openssl@3/3.1.0/lib REQUIRED) + target_link_libraries(nrniv_lib ${LIB_CRYPTO}) + else() + target_link_libraries(nrniv_lib crypto) + endif() +endif() + +if(${NRN_ENABLE_ARCH_INDEP_EXP_POW}) + find_library(LIB_MPFR mpfr REQUIRED) + target_link_libraries(nrniv_lib ${LIB_MPFR}) +endif() + if(NRN_WINDOWS_BUILD) target_link_libraries(nrniv_lib ${TERMCAP_LIBRARIES} ${Readline_LIBRARY}) else() diff --git a/src/nrniv/nmodlrandom.cpp b/src/nrniv/nmodlrandom.cpp index 66be7a9dc1..88806ddde9 100644 --- a/src/nrniv/nmodlrandom.cpp +++ b/src/nrniv/nmodlrandom.cpp @@ -100,12 +100,8 @@ static void nmodlrandom_destruct(void* v) { } void NMODLRandom_reg() { - class2oc("NMODLRandom", - nmodlrandom_cons, - nmodlrandom_destruct, - members, - retobj_members, - nullptr); + class2oc( + "NMODLRandom", nmodlrandom_cons, nmodlrandom_destruct, members, retobj_members, nullptr); if (!nmodlrandom_sym) { nmodlrandom_sym = hoc_lookup("NMODLRandom"); assert(nmodlrandom_sym); diff --git a/src/nrniv/nrnmenu.cpp b/src/nrniv/nrnmenu.cpp index f4103033cb..0f8e87f750 100644 --- a/src/nrniv/nrnmenu.cpp +++ b/src/nrniv/nrnmenu.cpp @@ -1129,8 +1129,7 @@ static Member_ret_obj_func mt_retobj_members[] = {{"pp_begin", mt_pp_begin}, {0, 0}}; static Member_ret_str_func mt_retstr_func[] = {{"code", mt_code}, {"file", mt_file}, {0, 0}}; void MechanismType_reg() { - class2oc( - "MechanismType", mt_cons, mt_destruct, mt_members, mt_retobj_members, mt_retstr_func); + class2oc("MechanismType", mt_cons, mt_destruct, mt_members, mt_retobj_members, mt_retstr_func); mt_class_sym_ = hoc_lookup("MechanismType"); } diff --git a/src/nrnpython/CMakeLists.txt b/src/nrnpython/CMakeLists.txt index c857e3c37b..0996b99c04 100644 --- a/src/nrnpython/CMakeLists.txt +++ b/src/nrnpython/CMakeLists.txt @@ -74,7 +74,8 @@ else() target_link_libraries(nrnpython ${NRN_DEFAULT_PYTHON_LIBRARIES}) endif() target_link_libraries(nrnpython fmt::fmt) - target_include_directories(nrnpython SYSTEM PUBLIC ${PROJECT_SOURCE_DIR}/${NRN_3RDPARTY_DIR}/eigen) + target_include_directories(nrnpython SYSTEM + PUBLIC ${PROJECT_SOURCE_DIR}/${NRN_3RDPARTY_DIR}/eigen) target_include_directories(nrnpython PUBLIC ${PROJECT_BINARY_DIR}/src/nrniv/oc_generated) make_nanobind_target(nanobind ${NRN_DEFAULT_PYTHON_INCLUDES}) target_link_libraries(nrnpython nanobind) diff --git a/src/oc/debug.cpp b/src/oc/debug.cpp index 2f8aeb5626..c9e0bb9988 100644 --- a/src/oc/debug.cpp +++ b/src/oc/debug.cpp @@ -1,12 +1,21 @@ #include <../../nrnconf.h> /* /local/src/master/nrn/src/oc/debug.cpp,v 1.7 1996/04/09 16:39:14 hines Exp */ + #include "hocdec.h" #include "code.h" #include "equation.h" +#include "multicore.h" #include #include "utils/logger.hpp" +#include "nrndigest.h" +#if NRN_DIGEST +#include +#include +#include +#endif + int hoc_zzdebug; #define prcod(c1, c2) else if (p->pf == c1) Printf("%p %p %s", fmt::ptr(p), fmt::ptr(p->pf), c2) @@ -144,3 +153,86 @@ void debugzz(Inst* p) { } #endif /*OCSMALL*/ } + +#if NRN_DIGEST + +int nrn_digest_; +static std::vector> digest; // nthread string vectors +static std::vector digest_cnt; // nthread counts. +static int nrn_digest_print_item_ = -1; +static int nrn_digest_print_tid_ = 0; +static bool nrn_digest_abort_ = false; + +void nrn_digest() { + if (ifarg(1) && hoc_is_str_arg(1)) { + // print the digest to the file and turn off accumulation + const char* fname = gargstr(1); + FILE* f = fopen(fname, "w"); + if (!f) { + hoc_execerr_ext("Could not open %s for writing", fname); + } + + int tid = 0; + for (auto& d: digest) { + fprintf(f, "tid=%d size=%zd\n", tid, digest[tid].size()); + for (auto& s: d) { + fprintf(f, "%s\n", s.c_str()); + } + tid++; + } + fclose(f); + nrn_digest_ = 0; + } else { // start accumulating digest info + nrn_digest_ = 1; + nrn_digest_print_item_ = -1; + nrn_digest_print_tid_ = 0; + if (ifarg(2)) { + nrn_digest_print_tid_ = int(chkarg(1, 0., nrn_nthread - 1)); + nrn_digest_print_item_ = int(chkarg(2, 0., 1e9)); + } + nrn_digest_abort_ = (ifarg(3) && strcmp(gargstr(3), "abort") == 0); + } + size_t size = digest.size() ? digest[0].size() : 0; + digest.clear(); // in any case, start over. + digest.resize(nrn_nthread); + digest_cnt.clear(); + digest_cnt.resize(nrn_nthread); + hoc_ret(); + hoc_pushx(double(size)); +} + +void nrn_digest_dbl_array(const char* msg, int tid, double t, double* array, size_t sz) { + if (!nrn_digest_) { + return; + } + unsigned char md[SHA_DIGEST_LENGTH]; + size_t n = sz * sizeof(double); + unsigned char* d = (unsigned char*) array; + SHA1(d, n, md); + + std::string s(msg); + char buf[100]; + int ix = int(digest_cnt[tid]); + digest_cnt[tid]++; + sprintf(buf, " %d %d %.17g ", tid, ix, t); + s += buf; + + for (int i = 0; i < 8; ++i) { + sprintf(buf, "%02x", (int) md[i]); + s += buf; + } + + digest[tid].push_back(s); + + if (nrn_digest_print_item_ == ix && nrn_digest_print_tid_ == tid) { + printf("ZZ %s\n", s.c_str()); + if (nrn_digest_abort_) { + abort(); + } + for (size_t i = 0; i < sz; ++i) { + printf("Z %zd %.20g\n", i, array[i]); + } + } +} + +#endif // NRN_DIGEST diff --git a/src/oc/hoc_init.cpp b/src/oc/hoc_init.cpp index bfcd5a9342..e8bd07f46c 100644 --- a/src/oc/hoc_init.cpp +++ b/src/oc/hoc_init.cpp @@ -195,6 +195,10 @@ static struct { /* Builtin functions with multiple or variable args */ #if defined(WIN32) {"WinExec", hoc_win_exec}, #endif +#if NRN_DIGEST + {"nrn_digest", nrn_digest}, +#endif + {"use_exp_pow_precision", hoc_use_exp_pow_precision}, {0, 0}}; static struct { /* functions that return a string */ diff --git a/src/oc/math.cpp b/src/oc/math.cpp index a568087679..2f2321fcf7 100644 --- a/src/oc/math.cpp +++ b/src/oc/math.cpp @@ -9,12 +9,16 @@ #include "nrnmpiuse.h" #include "ocfunc.h" +#include "nrnassrt.h" #include #include #include #include +#if NRN_ARCH_INDEP_EXP_POW +#include +#endif #define EPS hoc_epsilon #define MAXERRCOUNT 5 @@ -61,6 +65,109 @@ double hoc_Log10(double x) { return errcheck(log10(x), "log10"); } +#if NRN_ARCH_INDEP_EXP_POW + +static double accuracy32(double val) { + int ex; + double mant = frexp(val, &ex); + // round to about 32 bits after . + double prec = 4294967296.0; + double result = mant * prec; + result = round(result); + result /= prec; + result = ldexp(result, ex); + return result; +} + +static double pow_arch_indep(double x, double y) { + mpfr_prec_t prec = 53; + mpfr_rnd_t rnd = MPFR_RNDN; + + mpfr_t x_, y_; + mpfr_init2(x_, prec); + mpfr_init2(y_, prec); + + mpfr_set_d(x_, x, rnd); + mpfr_set_d(y_, y, rnd); + + mpfr_pow(x_, x_, y_, rnd); + + double r = mpfr_get_d(x_, rnd); + + mpfr_clear(x_); + mpfr_clear(y_); + + return r; +} + +static double exp_arch_indep(double x) { + mpfr_prec_t prec = 53; + mpfr_rnd_t rnd = MPFR_RNDN; + + mpfr_t x_; + mpfr_init2(x_, prec); + mpfr_set_d(x_, x, rnd); + mpfr_exp(x_, x_, rnd); + double r = mpfr_get_d(x_, rnd); + mpfr_clear(x_); + return r; +} + +static double pow_precision32(double x, double y) { + return accuracy32(pow(x, y)); +} + +static double exp_precision32(double x) { + return accuracy32(exp(x)); +} + +static double (*pow_ptr)(double x, double y) = pow; +static double (*pow_ieee_ptr)(double x, + double y) = pow; // avoid error! Maybe because of c++ overloading. +static double (*exp_ptr)(double x) = exp; + +int nrn_use_exp_pow_precision(int style) { + if (style == 0) { // default IEEE + pow_ptr = pow; + exp_ptr = exp; + } else if (style == 1) { // 53 bit mpfr + pow_ptr = pow_arch_indep; + exp_ptr = exp_arch_indep; + } else if (style == 2) { // 32 bit truncation + pow_ptr = pow_precision32; + exp_ptr = exp_precision32; + } + return style; +} + +#endif // NRN_ARCH_INDEP_EXP_POW + +void hoc_use_exp_pow_precision() { + int style = 0; + if (ifarg(1)) { + style = chkarg(1, 0.0, 2.0); + } +#if NRN_ARCH_INDEP_EXP_POW + style = nrn_use_exp_pow_precision(style); +#else + style = 0; +#endif // NRN_ARCH_INDEP_EXP_POW + hoc_ret(); + hoc_pushx(double(style)); +} + +// Try to overcome difference between linux and windows. +// by rounding mantissa to 32 bit accuracy or using +// hopefully arch independent mpfr for exp and pow. + +double hoc_pow(double x, double y) { +#if NRN_ARCH_INDEP_EXP_POW + return (*pow_ptr)(x, y); +#else + return pow(x, y); +#endif +} + /* used by nmodl and other c, c++ code */ double hoc_Exp(double x) { if (x < -700.) { @@ -78,7 +185,12 @@ double hoc_Exp(double x) { } return exp(700.); } + +#if NRN_ARCH_INDEP_EXP_POW + return (*exp_ptr)(x); +#else return exp(x); +#endif } /* used by interpreter */ @@ -103,7 +215,7 @@ double hoc_Sqrt(double x) { double hoc_Pow(double x, double y) { clear_fe_except(); - return errcheck(pow(x, y), "exponentiation"); + return errcheck(hoc_pow(x, y), "exponentiation"); } double hoc_integer(double x) { diff --git a/src/oc/nrndigest.h b/src/oc/nrndigest.h new file mode 100644 index 0000000000..0709460d0c --- /dev/null +++ b/src/oc/nrndigest.h @@ -0,0 +1,19 @@ + +#pragma once + +#include + +#ifdef __cplusplus +extern "C" { +#endif + +#define NRN_DIGEST NRN_ENABLE_DIGEST + +#if NRN_DIGEST +extern int nrn_digest_; // debugging differences on different machines. +extern void nrn_digest_dbl_array(const char* msg, int tid, double t, double* array, size_t sz); +#endif + +#ifdef __cplusplus +} +#endif diff --git a/src/oc/oc_ansi.h b/src/oc/oc_ansi.h index 0b95fc8732..70466b93bb 100644 --- a/src/oc/oc_ansi.h +++ b/src/oc/oc_ansi.h @@ -100,6 +100,7 @@ void install_vector_method(const char*, double (*)(void*)); int vector_arg_px(int i, double** p); double hoc_Exp(double); +extern "C" double hoc_pow(double, double); int hoc_is_tempobj_arg(int narg); std::FILE* hoc_obj_file_arg(int i); void hoc_reg_nmodl_text(int type, const char* txt); diff --git a/src/oc/ocfunc.h b/src/oc/ocfunc.h index b30006011a..67bf910d86 100644 --- a/src/oc/ocfunc.h +++ b/src/oc/ocfunc.h @@ -38,6 +38,12 @@ extern void hoc_Setcolor(void); extern void hoc_init_space(void); extern void hoc_install_hoc_obj(void); extern void nrn_feenableexcept(void); + +#if NRN_DIGEST +extern void nrn_digest(); +#endif +extern void hoc_use_exp_pow_precision(); + void hoc_coreneuron_handle(); void hoc_get_config_key(); void hoc_get_config_val(); diff --git a/src/sundials/shared/sundialsmath.c b/src/sundials/shared/sundialsmath.c index 972d3d82f1..84fb7e555e 100755 --- a/src/sundials/shared/sundialsmath.c +++ b/src/sundials/shared/sundialsmath.c @@ -39,12 +39,14 @@ realtype RPowerI(realtype base, int exponent) return(prod); } +extern double hoc_pow(double, double); + realtype RPowerR(realtype base, realtype exponent) { if (base <= ZERO) return(ZERO); #if defined(SUNDIALS_USE_GENERIC_MATH) - return((realtype) pow((double) base, (double) exponent)); + return((realtype) hoc_pow((double) base, (double) exponent)); #elif defined(SUNDIALS_SINGLE_PRECISION) return(powf(base, exponent)); #elif defined(SUNDIALS_EXTENDED_PRECISION) @@ -79,7 +81,7 @@ realtype RAbs(realtype x) realtype RPower2(realtype x) { #if defined(SUNDIALS_USE_GENERIC_MATH) - return((realtype) pow((double) x, 2.0)); + return((realtype) hoc_pow((double) x, 2.0)); #elif defined(SUNDIALS_SINGLE_PRECISION) return(powf(x, TWO)); #elif defined(SUNDIALS_EXTENDED_PRECISION) From 13e183b09419c1e064ed72dab7bf005d8c134b39 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 24 Oct 2024 10:28:38 +0200 Subject: [PATCH 11/17] Arguments for nrnpy_pyCallObject are created with nanobind (#3140) --- src/nrnpython/nrnpy_p2h.cpp | 43 +++++++++++++------------------------ 1 file changed, 15 insertions(+), 28 deletions(-) diff --git a/src/nrnpython/nrnpy_p2h.cpp b/src/nrnpython/nrnpy_p2h.cpp index 6bfae26924..c9dc98db92 100644 --- a/src/nrnpython/nrnpy_p2h.cpp +++ b/src/nrnpython/nrnpy_p2h.cpp @@ -177,28 +177,22 @@ static void py2n_component(Object* ob, Symbol* sym, int nindex, int isfunc) { PyErr_Print(); hoc_execerror("No attribute:", sym->name); } - PyObject* args = 0; Object* on; PyObject* result = 0; if (isfunc) { - args = PyTuple_New(nindex); + nb::list args{}; for (i = 0; i < nindex; ++i) { - PyObject* arg = nrnpy_hoc_pop("isfunc py2n_component"); + nb::object arg = nb::steal(nrnpy_hoc_pop("isfunc py2n_component")); if (!arg) { PyErr2NRNString e; e.get_pyerr(); - Py_DECREF(args); hoc_execerr_ext("arg %d error: %s", i, e.c_str()); } - // PyObject_Print(arg, stdout, 0); - // printf(" %d arg %d\n", arg->ob_refcnt, i); - if (PyTuple_SetItem(args, nindex - 1 - i, arg)) { - assert(0); - } + args.append(arg); } + args.reverse(); // printf("PyObject_CallObject %s %p\n", sym->name, tail); - result = nrnpy_pyCallObject(nb::borrow(tail), nb::borrow(args)).release().ptr(); - Py_DECREF(args); + result = nrnpy_pyCallObject(nb::borrow(tail), args).release().ptr(); // PyObject_Print(result, stdout, 0); // printf(" result of call\n"); if (!result) { @@ -216,7 +210,7 @@ static void py2n_component(Object* ob, Symbol* sym, int nindex, int isfunc) { return; } } else if (nindex) { - PyObject* arg; + nb::object arg; int n = hoc_pop_ndim(); if (n > 1) { hoc_execerr_ext( @@ -226,15 +220,15 @@ static void py2n_component(Object* ob, Symbol* sym, int nindex, int isfunc) { n); } if (hoc_stack_type() == NUMBER) { - arg = Py_BuildValue("l", (long) hoc_xpop()); + arg = nb::int_((long) hoc_xpop()); } else { // I don't think it is syntactically possible // for this to be a VAR. It is possible for it to // be an Object but the GetItem below will raise // TypeError: list indices must be integers or slices, not hoc.HocObject - arg = nrnpy_hoc_pop("nindex py2n_component"); + arg = nb::steal(nrnpy_hoc_pop("nindex py2n_component")); } - result = PyObject_GetItem(tail, arg); + result = PyObject_GetItem(tail, arg.ptr()); if (!result) { PyErr_Print(); hoc_execerror("Python get item failed:", hoc_object_name(ob)); @@ -468,24 +462,17 @@ static double func_call(Object* ho, int narg, int* err) { nb::callable po = nb::borrow(((Py2Nrn*) ho->u.this_pointer)->po_); nanobind::gil_scoped_acquire lock{}; - PyObject* args = PyTuple_New((Py_ssize_t) narg); - if (args == NULL) { - hoc_execerror("PyTuple_New failed", 0); - } + nb::list args{}; for (int i = 0; i < narg; ++i) { - PyObject* item = nrnpy_hoc_pop("func_call"); - if (item == NULL) { - Py_XDECREF(args); + nb::object item = nb::steal(nrnpy_hoc_pop("func_call")); + if (!item) { hoc_execerror("nrnpy_hoc_pop failed", 0); } - if (PyTuple_SetItem(args, (Py_ssize_t) (narg - i - 1), item) != 0) { - Py_XDECREF(args); - hoc_execerror("PyTuple_SetItem failed", 0); - } + args.append(item); } + args.reverse(); - nb::object r = nrnpy_pyCallObject(po, nb::borrow(args)); - Py_XDECREF(args); + nb::object r = nrnpy_pyCallObject(po, args); double rval = 0.0; if (!r.is_valid()) { if (!err || *err) { From f8e94c53da9d59570e6922acef6715f85776ff2e Mon Sep 17 00:00:00 2001 From: adamjhn Date: Thu, 24 Oct 2024 14:59:12 -0400 Subject: [PATCH 12/17] Move reference data from test_rxd.py to test_rxd.json. (#3143) * Added a --save flag to regenerate the data. --------- Co-authored-by: nrnhines --- setup.py | 2 +- share/lib/python/neuron/tests/test_rxd.json | 666 +++++++++++++++++ share/lib/python/neuron/tests/test_rxd.py | 750 ++------------------ 3 files changed, 730 insertions(+), 688 deletions(-) create mode 100644 share/lib/python/neuron/tests/test_rxd.json diff --git a/setup.py b/setup.py index ba0922b273..05e8ca67e1 100644 --- a/setup.py +++ b/setup.py @@ -500,7 +500,7 @@ def setup_package(): name=package_name, package_dir={"": NRN_PY_ROOT}, packages=py_packages, - package_data={"neuron": ["*.dat"]}, + package_data={"neuron": ["*.dat", "tests/*.json"]}, ext_modules=extensions, scripts=[ os.path.join(NRN_PY_SCRIPTS, f) diff --git a/share/lib/python/neuron/tests/test_rxd.json b/share/lib/python/neuron/tests/test_rxd.json new file mode 100644 index 0000000000..10c681c18f --- /dev/null +++ b/share/lib/python/neuron/tests/test_rxd.json @@ -0,0 +1,666 @@ +{ + "trivial_ecs_data": { + "False": [ + 1.0, + 0.9999975013886804, + 0.9999774378669442, + 0.9998977298459816, + 0.999683249239208, + 0.999233095122319, + 0.9984342775161097, + 0.9971750000657644, + 0.9953548976762606, + 0.9928916564339986, + 0.9897243754423741, + 0.985814368310196, + 0.9811441475925942, + 0.9757152507508027, + 0.9695454356132381, + 0.9626656387440524, + 0.9551169704970437, + 0.9469479227048566, + 0.9382118896237476, + 0.928965047713157, + 0.9192646017261445, + 0.9091673798703476, + 0.8987287461912615, + 0.8880017910312052, + 0.8770367581645542, + 0.8658806682496021, + 0.8545771012638642, + 0.8431661046665001, + 0.8316841985149556, + 0.8201644532317675, + 0.8086366199145236, + 0.7971272968665187, + 0.7856601193366658, + 0.7742559622858297, + 0.7629331483660402, + 0.7517076552500634, + 0.7405933180308322, + 0.729602023674572, + 0.7187438955075627, + 0.7080274664897305, + 0.6974598406190147, + 0.6870468422535746, + 0.6767931534640339, + 0.666702439759881, + 0.6567774646932262, + 0.6470201939460375, + 0.6374318895671022, + 0.6280131950530249, + 0.618764211972033, + 0.6096845688168512, + 0.6007734827485834, + 0.592029814861357, + 0.5834521195604474, + 0.5750386886069674, + 0.566787590341624, + 0.5586967045597302, + 0.5507637534704564, + 0.542986329135765, + 0.5353619177489981, + 0.5278879210798492, + 0.520561675381579, + 0.5133804680278274, + 0.5063415521201762, + 0.4994421592836728, + 0.4926795108456907, + 0.4860508275736653, + 0.47955333812926787, + 0.47318428638032217, + 0.4669409376970896, + 0.46082058434632084, + 0.4548205500845658, + 0.44893819404152296, + 0.4431709139745962, + 0.437516148967192, + 0.4319713816355469, + 0.4265341399019367, + 0.42120199838589845, + 0.4159725794595267, + 0.41084355400792044, + 0.40581264193139277, + 0.40087761242206477, + 0.3960362840438892, + 0.3912865246419576, + 0.38662625110408755, + 0.3820534289951289, + 0.3775660720821487, + 0.37316224176661034, + 0.3688400464378394, + 0.3645976407604392, + 0.3604332249068622, + 0.356345043745043, + 0.3523313859898391, + 0.34839058332598377, + 0.3445210095093366, + 0.3407210794523858, + 0.33698924829922156, + 0.33332401049454546, + 0.3297238988506902, + 0.32618748361610794, + 0.3227133715483166, + 0.319300204993885 + ], + "0.01": [ + 1.0, + 1.0, + 1.0, + 0.9999999999993757, + 0.999999999994894, + 0.9999999999684935, + 0.9999999997476527, + 0.9999999928933891, + 0.9999999611564773, + 0.9999998797767268, + 0.9999996998881439, + 0.9999993374567407, + 0.9999984853833065, + 0.9999969740580187, + 0.9999944882384337, + 0.999990631594901, + 0.9999849509581288, + 0.9999769206381163, + 0.9999582122018839, + 0.99992971528816, + 0.9998885520830331, + 0.9998315279906315, + 0.9996818635176349, + 0.9994522087548297, + 0.9991215537513378, + 0.9986688188931263, + 0.9980737640971887, + 0.9973169807956457, + 0.9963810719470815, + 0.9952501731851596, + 0.9939108695612834, + 0.9923517494757513, + 0.990563762843154, + 0.9885402079277492, + 0.9862765182500737, + 0.9837701825347182, + 0.9810207339245086, + 0.9780291325294735, + 0.9747989764631944, + 0.9713343413786614, + 0.9676399370120405, + 0.961066451151739, + 0.9539287981955907, + 0.9462557577841545, + 0.9380908351574285, + 0.9295048630248105, + 0.9272938750656939, + 0.9250579717836777, + 0.9227970935562565, + 0.916466538640017, + 0.9099911643714976, + 0.903362145740033, + 0.896607665570261, + 0.8897506341029571, + 0.882783639668767, + 0.875743025331908, + 0.8686235564310014, + 0.8614266253277882, + 0.8541863694353786, + 0.8468993821613305, + 0.8395705411269264, + 0.8322240727536439, + 0.8248512178576504, + 0.8174717412783322, + 0.8101012803420092, + 0.8027295128162794, + 0.795382570912665, + 0.7880520562862553, + 0.7807347758418224, + 0.7734449586662135, + 0.7661975818043348, + 0.7589799271016704, + 0.7518171027828684, + 0.7446977692713531, + 0.7376160661849058, + 0.7305954685552188, + 0.7236266087584206, + 0.7167068363521265, + 0.7140993659873981, + 0.7115003569088247, + 0.7053989759513706, + 0.6937157033260776, + 0.6892286429757224, + 0.6847738597616794, + 0.680346043973446, + 0.6736943892225341, + 0.6671180984407125, + 0.662236424428377, + 0.6573969190546712, + 0.6525933206930011, + 0.6478347503609198, + 0.6431166994074956, + 0.6384369406609263, + 0.629051615757503, + 0.6198222188768284, + 0.6175428046075583, + 0.6152735586037428, + 0.6111170811525735, + 0.6031212235932468, + 0.5952593537962105, + 0.5921305073271925, + 0.589021543706074, + 0.5843049292559478, + 0.5796335964525046, + 0.57501603797955, + 0.5704420014938749, + 0.5659178448897603, + 0.5614455517225058, + 0.5591788147005415, + 0.5569248277948933, + 0.5506074431348664, + 0.5478382336692825, + 0.5450878675432425, + 0.5387888979259772, + 0.5362404522913171, + 0.5337094641249713, + 0.5282374991063938, + 0.5259248453123049, + 0.5236262813729924, + 0.5182614252948405, + 0.5158555311460749, + 0.5134660721175922, + 0.5077709278071862, + 0.5053179819941784, + 0.5028817741460606, + 0.49744076421889366, + 0.49520127183112733, + 0.4929770110029178, + 0.4880637158772831, + 0.4859684631915766, + 0.4838861487508862, + 0.4790294757925702, + 0.47687065162757547, + 0.47472669028331466, + 0.4696849124713687, + 0.4675266379954105, + 0.4653827515472591, + 0.46059622403775785, + 0.45860850681188087, + 0.4566341643290534, + 0.45222826180913167, + 0.4503369914749967, + 0.4484573582920104, + 0.4440766399353272, + 0.44214365209440704, + 0.4402239033393764, + 0.4357529446107943, + 0.4338466468246234, + 0.4319526748954548, + 0.42771940093115685, + 0.4259480291636831, + 0.42418839478752657, + 0.4202300998323772, + 0.4185239107729718, + 0.4168280791675472, + 0.4128821813954716, + 0.411151841187283, + 0.40943320175716824, + 0.4054591322631254, + 0.4037685048216922, + 0.40208843869666, + 0.39832558823662206, + 0.39674079574945603, + 0.3951663234845347, + 0.3916025661247801, + 0.3900627879465062, + 0.38853217065630374, + 0.3849784483231059, + 0.38342810283841006, + 0.38188806435279576, + 0.3783453864065237, + 0.37683971858666915, + 0.3753431215943966, + 0.37198272430919915, + 0.3705596473205069, + 0.369145673011817, + 0.36593017118156906, + 0.36453932880580625, + 0.363156559348617, + 0.35995389116692217, + 0.358562523467025, + 0.3571802314864285, + 0.35401195406372227, + 0.3526654576925955, + 0.35132677928092004, + 0.34831289231094187, + 0.34703074630323383, + 0.3457566587383628, + 0.3428492728160761, + 0.34159142512656304, + 0.34034067750180974, + 0.3374508889913977, + 0.33619956600663337, + 0.3349562466940476, + 0.33211336362253635, + 0.3309044582410871, + 0.3297023177069934, + 0.32698870975793815, + 0.32583004833643353, + 0.3246785486149972, + 0.322044467574018, + 0.3209052709390766, + 0.31977231337146067 + ], + "1e-05": [ + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 0.9999999999999994, + 0.999999999999973, + 0.9999999999998286, + 0.9999999999994353, + 0.9999999999985202, + 0.9999999999965848, + 0.999999999991421, + 0.9999999999815654, + 0.999999999964278, + 0.999999999935795, + 0.9999999998244792, + 0.9999999995977931, + 0.9999999987558622, + 0.9999999969054565, + 0.9999999933717028, + 0.9999999871894636, + 0.9999999770945371, + 0.9999999615055236, + 0.9999999384936957, + 0.9999999057511308, + 0.9999998222204166, + 0.9999996889877318, + 0.9999994874837334, + 0.999999195494715, + 0.9999987869759616, + 0.9999978778927708, + 0.9999965094234947, + 0.9999945377974364, + 0.9999917976154723, + 0.9999881018395208, + 0.9999832421044165, + 0.99997698930736, + 0.9999690944252602, + 0.9999592894154012, + 0.9999401848217758, + 0.9999149868024333, + 0.9998825523031787, + 0.9998416677700266, + 0.999791059592578, + 0.9997294043566004, + 0.9996553389790277, + 0.9995674703998882, + 0.9994643850266823, + 0.9993446578768971, + 0.9992068613209312, + 0.9990495732666511, + 0.9988713847468418, + 0.9986709069017091, + 0.9984467774205655, + 0.9981976664057219, + 0.9979222816860815, + 0.9976193735297934, + 0.9972877388223238, + 0.9969262246908328, + 0.996533731653732, + 0.9961092162614316, + 0.9956516933112639, + 0.995160237588464, + 0.9946339852522628, + 0.9940721347820276, + 0.9934739476576494, + 0.9928387486094244, + 0.9921659257088754, + 0.9914549300102379, + 0.9907052751690002, + 0.9899165365563947, + 0.9890883505202034, + 0.9882204130405648, + 0.9873124787416938, + 0.9863643591450846, + 0.9853759215363472, + 0.9843470868598637, + 0.9832778285393545, + 0.9821681700454686, + 0.9810181837632249, + 0.9798279882523508, + 0.9785977472498147, + 0.9773276666315711, + 0.9760179936337368, + 0.9746690135138534, + 0.973281049070232, + 0.9718544574226932, + 0.970389628347756, + 0.9688869823638351, + 0.9662560899384262, + 0.9635203094788997, + 0.9606821993780256, + 0.9577444920812167, + 0.9547100690683178, + 0.9515819411366752, + 0.9483632227757351, + 0.9450571205495613, + 0.9416669020994433, + 0.9381959000999215, + 0.9346474739514014, + 0.9310250019996876, + 0.9273318651890746, + 0.9235714575523157, + 0.9197471440736972, + 0.9158622774695346, + 0.9119201741117782, + 0.9079241163948638, + 0.9038773101093882, + 0.8997829152661703, + 0.8956440506833048, + 0.8914637634703556, + 0.8872450198182428, + 0.8829907105064231, + 0.8787036576202755, + 0.8707492329414495, + 0.8687467372991504, + 0.8667391782827082, + 0.8647268018964069, + 0.862709850134354, + 0.8585978518713738, + 0.8544698540373321, + 0.8503277143424198, + 0.8461732532162092, + 0.8420082064910313, + 0.8378342644221053, + 0.8336530592192254, + 0.8294661713938566, + 0.8252750957253966, + 0.8210812732955917, + 0.8168861005800381, + 0.8126909232910099, + 0.808497028390498, + 0.801870527941527, + 0.7952548736200997, + 0.7886543612287468, + 0.7820730076362225, + 0.775514525626251, + 0.7689823476263041, + 0.76247965836786, + 0.756009378236416, + 0.749574210390875, + 0.743176642418555, + 0.7368189633320615, + 0.7305032399754438, + 0.724231346521673, + 0.7180049651772312, + 0.7118256478753264, + 0.7056948259233441, + 0.7032583360819944, + 0.700829860418756, + 0.6984094703765288, + 0.6959972339923495, + 0.6935932151265922, + 0.6889317429731345, + 0.6817668349582598, + 0.6746798126979188, + 0.6676718056770632, + 0.6607437321602112, + 0.6538963342916061, + 0.6471301405555971, + 0.6404455882961573, + 0.638787241606903, + 0.6371340113764097, + 0.6354859019485614, + 0.632317967018681, + 0.6291691359121044, + 0.624349241554268, + 0.619574789517913, + 0.614845708115443, + 0.6101619207501211, + 0.6055232995513015, + 0.6009296921089062, + 0.5963809255909179, + 0.5884265027887375, + 0.5806112287452613, + 0.5729337321697969, + 0.5653924576034643, + 0.5579857024960585, + 0.5561548262393686, + 0.5543322220611878, + 0.5525178614766851, + 0.5497619522652349, + 0.5470251014903538, + 0.5421958421956508, + 0.5374261214848611, + 0.5327153071436265, + 0.5280627574430141, + 0.5234678232476324, + 0.518929844460602, + 0.5144481596200524, + 0.5100220936639072, + 0.5056509822652775, + 0.501334134582734, + 0.497070898714232, + 0.492860583575257, + 0.488702496136303, + 0.4845959780644786, + 0.4805403496661883, + 0.4765349345681554, + 0.47257907275051975, + 0.46867209134329857, + 0.46481334795798307, + 0.4610021777856322, + 0.45723791746718145, + 0.4535199242983486, + 0.44984757092630373, + 0.4462202100082472, + 0.44263723745342903, + 0.4390980178177817, + 0.4356019204361891, + 0.43214835473501945, + 0.4287367135374875, + 0.4253663990336738, + 0.4220368326096732, + 0.41874742793163006, + 0.4154976298562691, + 0.41228686587129565, + 0.4091145679414439, + 0.40436315170089987, + 0.3996975358470319, + 0.3951158478239174, + 0.39061627014508404, + 0.3861970390853348, + 0.3841086119821148, + 0.3820380229802458, + 0.37998508572555983, + 0.3779496187417315, + 0.3745426919167079, + 0.3711844376968717, + 0.3678740151003091, + 0.3646106064268745, + 0.36139339169632445, + 0.3582215708412834, + 0.35509431859041374, + 0.35201085923465375, + 0.34897048346845855, + 0.345972428841367, + 0.3430159436132327, + 0.3401003270324417, + 0.3372248854371575, + 0.33438895086942516, + 0.33159183891956634, + 0.32883286088925856, + 0.32611135696096716, + 0.32342673756054513, + 0.320778362376132 + ] + }, + "scalar_bistable_data": [ + 4.666144368739565e-24, + 2.888704007378301e-23, + 1.9865049531958455e-22, + 1.3417088797559409e-21, + 8.872570814175612e-21, + 5.740880124662936e-20, + 3.632196361482048e-19, + 2.2456041210043948e-18, + 1.3557115052023306e-17, + 7.986451339137776e-17, + 4.587323676899687e-16, + 2.567045965818934e-15, + 1.398326895049454e-14, + 7.407949505967694e-14, + 3.8132509319176133e-13, + 1.905347304599464e-12, + 9.231798364410495e-12, + 4.33273265966826e-11, + 1.9674719569026998e-10, + 8.63399030466642e-10, + 3.657079221471436e-09, + 1.4932070567139546e-08, + 5.869401066025267e-08, + 2.218029569259483e-07, + 8.04721279962928e-07, + 2.79921382956858e-06, + 9.323183477925758e-06, + 2.969562407156746e-05, + 9.035408030381584e-05, + 0.0002623954841897342, + 0.0007269141545185264, + 0.0019207269099111785, + 0.004841879243431059, + 0.01164965343224167, + 0.02674863273052536, + 0.05846777500048251, + 0.12067994530088386, + 0.23084596756508963, + 0.3962758789592411, + 0.5900229199038892, + 0.7586218889021915, + 0.8722981880509894, + 0.9370823930113995, + 0.9705058492437157, + 0.9867204567968437, + 0.9942426940783657, + 0.9975977799681778, + 0.9990359504416326, + 0.9996249801006252, + 0.9998477834074035, + 0.9999041758657206, + 0.9998477834074035, + 0.9996249801006252, + 0.9990359504416326, + 0.9975977799681778, + 0.9942426940783657, + 0.9867204567968437, + 0.9705058492437147, + 0.9370823930113918, + 0.8722981880509718, + 0.75862188890219, + 0.5900229199038813, + 0.39627587895923605, + 0.2308459675650862, + 0.12067994530088234, + 0.058467775000481934, + 0.026748632730525093, + 0.011649653432241593, + 0.0048418792434310265, + 0.0019207269099111698, + 0.0007269141545185239, + 0.00026239548418973194, + 9.035408030381517e-05, + 2.9695624071567297e-05, + 9.323183477925707e-06, + 2.799213829568571e-06, + 8.04721279962927e-07, + 2.2180295692594684e-07, + 5.869401066025244e-08, + 1.4932070567139493e-08, + 3.6570792214714333e-09, + 8.633990304666437e-10, + 1.967471956902688e-10, + 4.332732659668245e-11, + 9.231798364410489e-12, + 1.9053473045994673e-12, + 3.813250931917625e-13, + 7.407949505967726e-14, + 1.3983268950494506e-14, + 2.5670459658189344e-15, + 4.587323676899693e-16, + 7.986451339137798e-17, + 1.3557115052023379e-17, + 2.2456041210043886e-18, + 3.632196361482047e-19, + 5.740880124662945e-20, + 8.872570814175644e-21, + 1.341708879755948e-21, + 1.9865049531958394e-22, + 2.8887040073782976e-23, + 4.66614436873957e-24 + ] +} \ No newline at end of file diff --git a/share/lib/python/neuron/tests/test_rxd.py b/share/lib/python/neuron/tests/test_rxd.py index db85696c90..dc6cd2aef8 100644 --- a/share/lib/python/neuron/tests/test_rxd.py +++ b/share/lib/python/neuron/tests/test_rxd.py @@ -2,7 +2,9 @@ import neuron import unittest import sys -from multiprocessing import Process +import os +import json +from multiprocessing import Process, Lock try: import multiprocessing as mp @@ -11,674 +13,12 @@ except: pass -scalar_bistable_data = [ - 4.666144368739553e-24, - 2.888704007378294e-23, - 1.986504953195841e-22, - 1.341708879755938e-21, - 8.872570814175589e-21, - 5.740880124662921e-20, - 3.632196361482038e-19, - 2.245604121004388e-18, - 1.355711505202327e-17, - 7.986451339137754e-17, - 4.587323676899676e-16, - 2.567045965818926e-15, - 1.398326895049450e-14, - 7.407949505967670e-14, - 3.813250931917600e-13, - 1.905347304599457e-12, - 9.231798364410461e-12, - 4.332732659668245e-11, - 1.967471956902693e-10, - 8.633990304666386e-10, - 3.657079221471421e-09, - 1.493207056713948e-08, - 5.869401066025243e-08, - 2.218029569259474e-07, - 8.047212799629250e-07, - 2.799213829568570e-06, - 9.323183477925731e-06, - 2.969562407156739e-05, - 9.035408030381566e-05, - 2.623954841897339e-04, - 7.269141545185255e-04, - 1.920726909911178e-03, - 4.841879243431064e-03, - 1.164965343224173e-02, - 2.674863273052559e-02, - 5.846777500048252e-02, - 1.206799453008834e-01, - 2.308459675650935e-01, - 3.962758789592548e-01, - 5.900229199039158e-01, - 7.586218889022415e-01, - 8.722981880510015e-01, - 9.370823930114011e-01, - 9.705058492437171e-01, - 9.867204567968444e-01, - 9.942426940783661e-01, - 9.975977799681778e-01, - 9.990359504416327e-01, - 9.996249801006252e-01, - 9.998477834074035e-01, - 9.999041758657206e-01, - 9.998477834074035e-01, - 9.996249801006252e-01, - 9.990359504416326e-01, - 9.975977799681777e-01, - 9.942426940783655e-01, - 9.867204567968437e-01, - 9.705058492437160e-01, - 9.370823930113995e-01, - 8.722981880509845e-01, - 7.586218889021992e-01, - 5.900229199038706e-01, - 3.962758789592359e-01, - 2.308459675650852e-01, - 1.206799453008814e-01, - 5.846777500048142e-02, - 2.674863273052497e-02, - 1.164965343224158e-02, - 4.841879243431020e-03, - 1.920726909911166e-03, - 7.269141545185224e-04, - 2.623954841897313e-04, - 9.035408030381501e-05, - 2.969562407156726e-05, - 9.323183477925702e-06, - 2.799213829568569e-06, - 8.047212799629269e-07, - 2.218029569259469e-07, - 5.869401066025247e-08, - 1.493207056713951e-08, - 3.657079221471437e-09, - 8.633990304666446e-10, - 1.967471956902691e-10, - 4.332732659668252e-11, - 9.231798364410503e-12, - 1.905347304599471e-12, - 3.813250931917631e-13, - 7.407949505967741e-14, - 1.398326895049453e-14, - 2.567045965818940e-15, - 4.587323676899705e-16, - 7.986451339137816e-17, - 1.355711505202341e-17, - 2.245604121004394e-18, - 3.632196361482056e-19, - 5.740880124662959e-20, - 8.872570814175665e-21, - 1.341708879755951e-21, - 1.986504953195844e-22, - 2.888704007378305e-23, - 4.666144368739581e-24, -] +# load the reference data from rxd_data.json +fdir = os.path.dirname(os.path.abspath(__file__)) +test_data = json.load(open(os.path.join(fdir, "test_rxd.json"), "r")) -trivial_ecs_data = { - False: [ - 1.000000000000000e00, - 9.999975013886804e-01, - 9.999774378669442e-01, - 9.998977298459814e-01, - 9.996832492392076e-01, - 9.992330951223182e-01, - 9.984342775161091e-01, - 9.971750000657639e-01, - 9.953548976762590e-01, - 9.928916564339932e-01, - 9.897243754423555e-01, - 9.858143683101370e-01, - 9.811441475924241e-01, - 9.757152507503439e-01, - 9.695454356120868e-01, - 9.626656387413414e-01, - 9.551169704910168e-01, - 9.469479226921377e-01, - 9.382118895981384e-01, - 9.289650476637475e-01, - 9.192646016344460e-01, - 9.091673797060952e-01, - 8.987287459064514e-01, - 8.880017905518656e-01, - 8.770367573796769e-01, - 8.658806669966089e-01, - 8.545770993099358e-01, - 8.431661016850746e-01, - 8.316841940567001e-01, - 8.201644466893430e-01, - 8.086366104805101e-01, - 7.971272834839320e-01, - 7.856601006415908e-01, - 7.742559365417961e-01, - 7.629331133899011e-01, - 7.517076083292696e-01, - 7.405932558321451e-01, - 7.296019421444131e-01, - 7.187437897643434e-01, - 7.080273307086645e-01, - 6.974596679100502e-01, - 6.870466245332152e-01, - 6.767928813219394e-01, - 6.667021023212317e-01, - 6.567770494779278e-01, - 6.470196867258985e-01, - 6.374312742221647e-01, - 6.280124534282580e-01, - 6.187633237355973e-01, - 6.096835113211366e-01, - 6.007722308951912e-01, - 5.920283409711493e-01, - 5.834503932497362e-01, - 5.750366766708355e-01, - 5.667852566452943e-01, - 5.586940099388136e-01, - 5.507606556408047e-01, - 5.429827826135633e-01, - 5.353578737816238e-01, - 5.278833275879240e-01, - 5.205564769125375e-01, - 5.133746057212193e-01, - 5.063349636848303e-01, - 4.994347789867492e-01, - 4.926712695135610e-01, - 4.860416526044784e-01, - 4.795431535169798e-01, - 4.731730127498988e-01, - 4.669284923505265e-01, - 4.608068813190689e-01, - 4.548055002118984e-01, - 4.489217050343421e-01, - 4.431528905041363e-01, - 4.374964927580476e-01, - 4.319499915664357e-01, - 4.265109121135835e-01, - 4.211768263954196e-01, - 4.159453542806875e-01, - 4.108141642766345e-01, - 4.057809740358395e-01, - 4.008435506368045e-01, - 3.959997106673694e-01, - 3.912473201368166e-01, - 3.865842942396779e-01, - 3.820085969917059e-01, - 3.775182407561858e-01, - 3.731112856767305e-01, - 3.687858390308746e-01, - 3.645400545171554e-01, - 3.603721314869148e-01, - 3.562803141307546e-01, - 3.522628906284160e-01, - 3.483181922698216e-01, - 3.444445925540838e-01, - 3.406405062724689e-01, - 3.369043885805584e-01, - 3.332347340641985e-01, - 3.296300758032397e-01, - 3.260889844365475e-01, - 3.226100672312980e-01, - 3.191919671591613e-01, - ], - 1e-2: [ - 1.000000000000000e00, - 1.000000000000000e00, - 1.000000000000000e00, - 9.999999999993757e-01, - 9.999999999948940e-01, - 9.999999999684935e-01, - 9.999999997476527e-01, - 9.999999928933891e-01, - 9.999999611564773e-01, - 9.999998797767268e-01, - 9.999996998881439e-01, - 9.999993374567406e-01, - 9.999984853833063e-01, - 9.999969740580184e-01, - 9.999944882384333e-01, - 9.999906315949002e-01, - 9.999849509581277e-01, - 9.999769206381147e-01, - 9.999582122018814e-01, - 9.999297152881566e-01, - 9.998885520830283e-01, - 9.998315279906250e-01, - 9.996818635176246e-01, - 9.994522087548142e-01, - 9.991215537513162e-01, - 9.986688188930973e-01, - 9.980737640971513e-01, - 9.973169807955987e-01, - 9.963810719470240e-01, - 9.952501731850908e-01, - 9.939108695612021e-01, - 9.923517494756579e-01, - 9.905637628430480e-01, - 9.885402079276301e-01, - 9.862765182499404e-01, - 9.837701825345700e-01, - 9.810207339243457e-01, - 9.780291325292961e-01, - 9.747989764630028e-01, - 9.713343413784561e-01, - 9.676399370118218e-01, - 9.610664511514549e-01, - 9.539287981952346e-01, - 9.462557577837218e-01, - 9.380908351569157e-01, - 9.295048630242134e-01, - 9.272938750650753e-01, - 9.250579717830375e-01, - 9.227970935555948e-01, - 9.164665386395120e-01, - 9.099911643711576e-01, - 9.033621457398684e-01, - 8.966076655702790e-01, - 8.897506341031615e-01, - 8.827836396691656e-01, - 8.757430253325005e-01, - 8.686235564317917e-01, - 8.614266253287830e-01, - 8.541863694365750e-01, - 8.468993821627302e-01, - 8.395705411285332e-01, - 8.322240727554534e-01, - 8.248512178596673e-01, - 8.174717412805539e-01, - 8.101012803444277e-01, - 8.027295128188973e-01, - 7.953825709154709e-01, - 7.880520562892473e-01, - 7.807347758450052e-01, - 7.734449586695845e-01, - 7.661975818078829e-01, - 7.589799271053987e-01, - 7.518171027867613e-01, - 7.446977692754075e-01, - 7.376160661891280e-01, - 7.305954685595947e-01, - 7.236266087629475e-01, - 7.167068363568081e-01, - 7.140993659912696e-01, - 7.115003569118910e-01, - 7.053989759515575e-01, - 6.937157033233743e-01, - 6.892286429728892e-01, - 6.847738597587193e-01, - 6.803460439703601e-01, - 6.736943892201644e-01, - 6.671180984390447e-01, - 6.622364244386257e-01, - 6.573969190766189e-01, - 6.525933207264888e-01, - 6.478347504056896e-01, - 6.431166994633450e-01, - 6.384369407276960e-01, - 6.290516157893640e-01, - 6.198222188748919e-01, - 6.175428045973979e-01, - 6.152735585854211e-01, - 6.111170811451881e-01, - 6.031212236479893e-01, - 5.952593539110380e-01, - 5.921305074751056e-01, - 5.890215438866859e-01, - 5.843049294344698e-01, - 5.796335966290896e-01, - 5.750160381540556e-01, - 5.704420016665088e-01, - 5.659178450605010e-01, - 5.614455518912613e-01, - 5.591788148460601e-01, - 5.569248279174298e-01, - 5.506074431743505e-01, - 5.478382336868175e-01, - 5.450878675391179e-01, - 5.387888979310328e-01, - 5.362404523147913e-01, - 5.337094641666161e-01, - 5.282374991972423e-01, - 5.259248454152894e-01, - 5.236262814879684e-01, - 5.182614254027890e-01, - 5.158555312376837e-01, - 5.134660721930772e-01, - 5.077709278344216e-01, - 5.053179820104490e-01, - 5.028817741515068e-01, - 4.974407642335685e-01, - 4.952012718592664e-01, - 4.929770110443390e-01, - 4.880637159524582e-01, - 4.859684632733423e-01, - 4.838861488391641e-01, - 4.790294758699501e-01, - 4.768706516921771e-01, - 4.747266903353097e-01, - 4.696849124915995e-01, - 4.675266380099869e-01, - 4.653827515562493e-01, - 4.605962240576634e-01, - 4.586085068420688e-01, - 4.566341643693850e-01, - 4.522282618718538e-01, - 4.503369915407668e-01, - 4.484573583607938e-01, - 4.440766399922996e-01, - 4.421436521416536e-01, - 4.402239033770282e-01, - 4.357529446279363e-01, - 4.338466468392896e-01, - 4.319526749076715e-01, - 4.277194009540304e-01, - 4.259480291942094e-01, - 4.241883948256027e-01, - 4.202300998848154e-01, - 4.185239108263153e-01, - 4.168280792217872e-01, - 4.128821814385824e-01, - 4.111518412231812e-01, - 4.094332017859506e-01, - 4.054591322790473e-01, - 4.037685048369686e-01, - 4.020884387112946e-01, - 3.983255882607288e-01, - 3.967407957791139e-01, - 3.951663235196693e-01, - 3.916025661686941e-01, - 3.900627879901100e-01, - 3.885321706996029e-01, - 3.849784483567140e-01, - 3.834281028668022e-01, - 3.818880643760414e-01, - 3.783453864220350e-01, - 3.768397186025170e-01, - 3.753431216105738e-01, - 3.719827243333056e-01, - 3.705596473485314e-01, - 3.691456730437073e-01, - 3.659301712184497e-01, - 3.645393288417587e-01, - 3.631565593836533e-01, - 3.599538911939479e-01, - 3.585625234903733e-01, - 3.571802315061479e-01, - 3.540119540790766e-01, - 3.526654577087503e-01, - 3.513267792978637e-01, - 3.483128923342382e-01, - 3.470307463292109e-01, - 3.457566587669847e-01, - 3.428492728471916e-01, - 3.415914251565042e-01, - 3.403406775305906e-01, - 3.374508890138111e-01, - 3.361995660265230e-01, - 3.349562467114465e-01, - 3.321133636377239e-01, - 3.309044582572311e-01, - 3.297023177240813e-01, - 3.269887097799309e-01, - 3.258300483601927e-01, - 3.246785486404988e-01, - 3.220444676004350e-01, - 3.209052709642934e-01, - 3.197723133954909e-01, - ], - 1e-5: [ - 1.000000000000000e00, - 1.000000000000000e00, - 1.000000000000000e00, - 1.000000000000000e00, - 1.000000000000000e00, - 1.000000000000000e00, - 1.000000000000000e00, - 1.000000000000000e00, - 9.999999999999994e-01, - 9.999999999999730e-01, - 9.999999999998286e-01, - 9.999999999994353e-01, - 9.999999999985202e-01, - 9.999999999965848e-01, - 9.999999999914210e-01, - 9.999999999815654e-01, - 9.999999999642780e-01, - 9.999999999357950e-01, - 9.999999998244792e-01, - 9.999999995977931e-01, - 9.999999987558622e-01, - 9.999999969054565e-01, - 9.999999933717028e-01, - 9.999999871894636e-01, - 9.999999770945371e-01, - 9.999999615055236e-01, - 9.999999384936957e-01, - 9.999999057511308e-01, - 9.999998222204165e-01, - 9.999996889877317e-01, - 9.999994874837332e-01, - 9.999991954947148e-01, - 9.999987869759613e-01, - 9.999978778927705e-01, - 9.999965094234943e-01, - 9.999945377974357e-01, - 9.999917976154713e-01, - 9.999881018395195e-01, - 9.999832421044147e-01, - 9.999769893073577e-01, - 9.999690944252573e-01, - 9.999592894153977e-01, - 9.999401848217710e-01, - 9.999149868024271e-01, - 9.998825523031706e-01, - 9.998416677700163e-01, - 9.997910595925650e-01, - 9.997294043565842e-01, - 9.996553389790078e-01, - 9.995674703998642e-01, - 9.994643850266534e-01, - 9.993446578768627e-01, - 9.992068613208905e-01, - 9.990495732666037e-01, - 9.988713847467867e-01, - 9.986709069016458e-01, - 9.984467774204933e-01, - 9.981976664056399e-01, - 9.979222816859891e-01, - 9.976193735296895e-01, - 9.972877388222079e-01, - 9.969262246907040e-01, - 9.965337316535895e-01, - 9.961092162612745e-01, - 9.956516933110915e-01, - 9.951602375882754e-01, - 9.946339852520572e-01, - 9.940721347818040e-01, - 9.934739476574072e-01, - 9.928387486091628e-01, - 9.921659257085935e-01, - 9.914549300099351e-01, - 9.907052751686757e-01, - 9.899165365560477e-01, - 9.890883505198331e-01, - 9.882204130401707e-01, - 9.873124787412751e-01, - 9.863643591446407e-01, - 9.853759215358775e-01, - 9.843470868593674e-01, - 9.832778285388312e-01, - 9.821681700449177e-01, - 9.810181837626457e-01, - 9.798279882517430e-01, - 9.785977472491776e-01, - 9.773276666309043e-01, - 9.760179936330398e-01, - 9.746690135131258e-01, - 9.732810490694735e-01, - 9.718544574219034e-01, - 9.703896283469345e-01, - 9.688869823629818e-01, - 9.662560899379753e-01, - 9.635203094788860e-01, - 9.606821993784824e-01, - 9.577444920821762e-01, - 9.547100690698106e-01, - 9.515819411387303e-01, - 9.483632227783803e-01, - 9.450571205528219e-01, - 9.416669021033438e-01, - 9.381959001044843e-01, - 9.346474739566473e-01, - 9.310250020056358e-01, - 9.273318651957428e-01, - 9.235714575597201e-01, - 9.197471440818522e-01, - 9.158622774784529e-01, - 9.119201741214706e-01, - 9.079241164053400e-01, - 9.038773101206564e-01, - 8.997829152782374e-01, - 8.956440506961761e-01, - 8.914637634840353e-01, - 8.872450198327338e-01, - 8.829907105217275e-01, - 8.787036576363941e-01, - 8.707492329602121e-01, - 8.687467373185769e-01, - 8.667391783027991e-01, - 8.647268019171629e-01, - 8.627098501557753e-01, - 8.585978518951599e-01, - 8.544698540634895e-01, - 8.503277143709528e-01, - 8.461732532471191e-01, - 8.420082065243172e-01, - 8.378342644577640e-01, - 8.336530592572513e-01, - 8.294661714342418e-01, - 8.252750957681314e-01, - 8.210812733406646e-01, - 8.168861006274361e-01, - 8.126909233407180e-01, - 8.084970284424999e-01, - 8.018705279973948e-01, - 7.952548736797830e-01, - 7.886543612921902e-01, - 7.820730077033664e-01, - 7.755145256970315e-01, - 7.689823477006539e-01, - 7.624796584457091e-01, - 7.560093783176929e-01, - 7.495742104755058e-01, - 7.431766425064636e-01, - 7.368189634231701e-01, - 7.305032400696735e-01, - 7.242313466189446e-01, - 7.180049652774664e-01, - 7.118256479784464e-01, - 7.056948260292686e-01, - 7.032583361895295e-01, - 7.008298605278855e-01, - 6.984094704872370e-01, - 6.959972341046202e-01, - 6.935932152404094e-01, - 6.889317430774708e-01, - 6.817668350489530e-01, - 6.746798127752269e-01, - 6.676718057412527e-01, - 6.607437322115557e-01, - 6.538963343303841e-01, - 6.471301405820917e-01, - 6.404455883106537e-01, - 6.387872416184449e-01, - 6.371340113850150e-01, - 6.354859019542481e-01, - 6.323179670207283e-01, - 6.291691359105558e-01, - 6.243492415503886e-01, - 6.195747895117473e-01, - 6.148457081070352e-01, - 6.101619207395151e-01, - 6.055232995385412e-01, - 6.009296920940354e-01, - 5.963809255739797e-01, - 5.884265027623434e-01, - 5.806112287097513e-01, - 5.729337321255051e-01, - 5.653924575507177e-01, - 5.579857024351769e-01, - 5.561548261765027e-01, - 5.543322219963571e-01, - 5.525178614099090e-01, - 5.497619521958399e-01, - 5.470251014183799e-01, - 5.421958421339816e-01, - 5.374261214332550e-01, - 5.327153071018452e-01, - 5.280627574108234e-01, - 5.234678232248013e-01, - 5.189298444469039e-01, - 5.144481596152644e-01, - 5.100220936678099e-01, - 5.056509822776560e-01, - 5.013341346033774e-01, - 4.970708987429328e-01, - 4.928605836118120e-01, - 4.887024961805134e-01, - 4.845959781161491e-01, - 4.805403497251276e-01, - 4.765349346341766e-01, - 4.725790728234394e-01, - 4.686720914229373e-01, - 4.648133480441644e-01, - 4.610021778781840e-01, - 4.572379175659361e-01, - 4.535199244031417e-01, - 4.498475710369739e-01, - 4.462202101246375e-01, - 4.426372375753840e-01, - 4.390980179451507e-01, - 4.356019205688260e-01, - 4.321483548727800e-01, - 4.287367136802313e-01, - 4.253663991812640e-01, - 4.220368327619753e-01, - 4.187474280885136e-01, - 4.154976300176047e-01, - 4.122868660369579e-01, - 4.091145681113114e-01, - 4.043631518312791e-01, - 3.996975359393022e-01, - 3.951158478794149e-01, - 3.906162701651016e-01, - 3.861970390711226e-01, - 3.841086119708872e-01, - 3.820380229719526e-01, - 3.799850857201517e-01, - 3.779496187391600e-01, - 3.745426919402392e-01, - 3.711844377457629e-01, - 3.678740151738357e-01, - 3.646106065243293e-01, - 3.613933918170190e-01, - 3.582215709845467e-01, - 3.550943187555968e-01, - 3.520108594211252e-01, - 3.489704836755959e-01, - 3.459724290685680e-01, - 3.430159438599149e-01, - 3.401003272980346e-01, - 3.372248857211048e-01, - 3.343889511711801e-01, - 3.315918392385992e-01, - 3.288328612250609e-01, - 3.261113573130471e-01, - 3.234267379284118e-01, - 3.207783627593125e-01, - ], -} - -def scalar_bistable(): +def scalar_bistable(lock, path=None): from neuron import rxd h.load_file("stdrun.hoc") @@ -694,15 +34,25 @@ def scalar_bistable(): # check the results result = h.Vector(c.nodes.concentration) - cmpV = h.Vector(scalar_bistable_data) - cmpV.sub(result) - cmpV.abs() - if cmpV.sum() < 1e-6: - sys.exit(0) - sys.exit(-1) - - -def trivial_ecs(scale): + if path is not None: + lock.acquire() + if os.path.exists(path): + data = json.load(open(path, "r")) + else: + data = {} + data["scalar_bistable_data"] = list(result) + json.dump(data, open(path, "w"), indent=4) + lock.release() + else: + cmpV = h.Vector(test_data["scalar_bistable_data"]) + cmpV.sub(result) + cmpV.abs() + if cmpV.sum() >= 1e-6: + sys.exit(-1) + sys.exit(0) + + +def trivial_ecs(scale, lock, path=None): from neuron import h, crxd as rxd import numpy import warnings @@ -714,6 +64,7 @@ def trivial_ecs(scale): h.CVode().active(True) h.CVode().event(tstop) else: # fixed step case + h.CVode().active(False) h.dt = 0.1 sec = h.Section() # NEURON requires at least 1 section @@ -757,40 +108,63 @@ def trivial_ecs(scale): h.finitialize() h.continuerun(tstop) # run the simulation - # compare with previous solution - ecs_vec.sub(h.Vector(trivial_ecs_data[scale])) - ecs_vec.abs() - if ecs_vec.sum() > 1e-9: - return -1 - return 0 + if path is not None: + lock.acquire() + if os.path.exists(path): + data = json.load(open(path, "r")) + else: + data = {} + if "trivial_ecs_data" not in data: + data["trivial_ecs_data"] = {} + data["trivial_ecs_data"][str(scale)] = list(ecs_vec) + json.dump(data, open(path, "w"), indent=4) + lock.release() + else: + # compare with previous solution + ecs_vec.sub(h.Vector(test_data["trivial_ecs_data"][str(scale)])) + ecs_vec.abs() + if ecs_vec.sum() > 1e-9: + sys.exit(-1) + sys.exit(0) class RxDTestCase(unittest.TestCase): """Tests of rxd""" + @classmethod + def setUpClass(cls): + # Check for --save in command-line arguments + cls.path = None + cls.lock = Lock() + if len(sys.argv) > 1: + for arg1, arg2 in zip(sys.argv, sys.argv[1:]): + if arg1 == "--save": + cls.path = arg2 + break + def test_rxd(self): - p = Process(target=scalar_bistable) + p = Process(target=scalar_bistable, args=(self.lock, self.path)) p.start() p.join() assert p.exitcode == 0 return 0 def test_ecs_diffusion_fixed_step(self): - p = Process(target=trivial_ecs, args=(False,)) + p = Process(target=trivial_ecs, args=(False, self.lock, self.path)) p.start() p.join() assert p.exitcode == 0 return 0 def test_ecs_diffusion_variable_step_coarse(self): - p = Process(target=trivial_ecs, args=(1e-2,)) + p = Process(target=trivial_ecs, args=(1e-2, self.lock, self.path)) p.start() p.join() assert p.exitcode == 0 return 0 def test_ecs_diffusion_variable_step_fine(self): - p = Process(target=trivial_ecs, args=(1e-5,)) + p = Process(target=trivial_ecs, args=(1e-5, self.lock, self.path)) p.start() p.join() assert p.exitcode == 0 @@ -807,5 +181,7 @@ def test(): if __name__ == "__main__": - # unittest.main() + """Run the rxd tests unless --save is given, in which case save the + test data to that path without comparing the results to the reference data. + """ test() From cd7ac7dd9635384e3cce28cc7dcc482afa98ea54 Mon Sep 17 00:00:00 2001 From: nrnhines Date: Thu, 24 Oct 2024 20:41:03 -0400 Subject: [PATCH 13/17] Code change for ion and ion style semantics. Unlimited numbers. (#3097) --------- Co-authored-by: Luc Grosheintz --- cmake/NeuronFileLists.cmake | 1 + src/coreneuron/io/nrn_checkpoint.cpp | 6 ++++-- src/coreneuron/io/phase2.cpp | 5 +++-- src/coreneuron/mechanism/register_mech.cpp | 9 +++++---- src/coreneuron/permute/node_permute.cpp | 5 +++-- src/neuron/cache/mechanism_range.hpp | 2 +- .../nrncore_write/callbacks/nrncore_callbacks.cpp | 6 +++--- src/nrniv/nrncore_write/data/cell_group.cpp | 6 +++--- src/nrniv/nrncore_write/data/cell_group.h | 4 ++-- src/nrnoc/init.cpp | 4 ++-- src/nrnoc/ion_semantics.h | 14 ++++++++++++++ src/nrnoc/membfunc.h | 1 + 12 files changed, 42 insertions(+), 21 deletions(-) create mode 100644 src/nrnoc/ion_semantics.h diff --git a/cmake/NeuronFileLists.cmake b/cmake/NeuronFileLists.cmake index d8b03dae3a..47dd1dcc93 100644 --- a/cmake/NeuronFileLists.cmake +++ b/cmake/NeuronFileLists.cmake @@ -21,6 +21,7 @@ set(HEADER_FILES_TO_INSTALL nrnoc/md1redef.h nrnoc/md2redef.h nrnoc/membdef.h + nrnoc/ion_semantics.h nrnoc/membfunc.h nrnoc/multicore.h nrnoc/multisplit.h diff --git a/src/coreneuron/io/nrn_checkpoint.cpp b/src/coreneuron/io/nrn_checkpoint.cpp index ed91d7de81..c7c7743ed9 100644 --- a/src/coreneuron/io/nrn_checkpoint.cpp +++ b/src/coreneuron/io/nrn_checkpoint.cpp @@ -11,6 +11,7 @@ #include #include +#include "nrnoc/ion_semantics.h" #include "coreneuron/sim/multicore.hpp" #include "coreneuron/nrniv/nrniv_decl.h" #include "coreneuron/io/nrn_filehandler.hpp" @@ -296,8 +297,9 @@ void CheckPoints::write_phase2(NrnThread& nt) const { // out into the following function. d[ix] = nrn_original_aos_index(ptype, d[ix], nt, ml_pinv); } - } else if (s >= 0 && s < 1000) { // ion - d[ix] = nrn_original_aos_index(s, d[ix], nt, ml_pinv); + } else if (nrn_semantics_is_ion(s)) { // ion + auto type = nrn_semantics_ion_type(s); + d[ix] = nrn_original_aos_index(type, d[ix], nt, ml_pinv); } #if CHKPNTDEBUG if (s != -8) { // WATCH values change diff --git a/src/coreneuron/io/phase2.cpp b/src/coreneuron/io/phase2.cpp index 19fb12e8a3..3c712f1f5b 100644 --- a/src/coreneuron/io/phase2.cpp +++ b/src/coreneuron/io/phase2.cpp @@ -6,6 +6,7 @@ # ============================================================================= */ +#include "nrnoc/ion_semantics.h" #include "coreneuron/io/phase2.hpp" #include "coreneuron/coreneuron.hpp" #include "coreneuron/sim/multicore.hpp" @@ -699,8 +700,8 @@ void Phase2::pdata_relocation(const NrnThread& nt, const std::vector& **/ break; default: - if (s >= 0 && s < 1000) { // ion - int etype = s; + if (nrn_semantics_is_ion(s)) { // ion + int etype = nrn_semantics_ion_type(s); /* if ion is SoA, must recalculate pdata values */ /* if ion is AoS, have to deal with offset */ Memb_list* eml = nt._ml_list[etype]; diff --git a/src/coreneuron/mechanism/register_mech.cpp b/src/coreneuron/mechanism/register_mech.cpp index 260d420411..2ef95e976f 100644 --- a/src/coreneuron/mechanism/register_mech.cpp +++ b/src/coreneuron/mechanism/register_mech.cpp @@ -9,6 +9,7 @@ #include #include "coreneuron/nrnconf.h" +#include "nrnoc/ion_semantics.h" #include "coreneuron/sim/multicore.hpp" #include "coreneuron/membrane_definitions.h" #include "coreneuron/mechanism/eion.hpp" @@ -212,7 +213,7 @@ void hoc_register_dparam_semantics(int type, int ix, const char* name) { xx_ion and #xx_ion which will get a semantics value of -1, -2, -3, -4, -5, -6, -7, -8, -9, -10, - type, and type+1000 respectively + 2*type, and 2*type+1 respectively */ auto& memb_func = corenrn.get_memb_funcs(); if (strcmp(name, "area") == 0) { @@ -240,7 +241,7 @@ void hoc_register_dparam_semantics(int type, int ix, const char* name) { } else { int i = name[0] == '#' ? 1 : 0; int etype = nrn_get_mechtype(name + i); - memb_func[type].dparam_semantics[ix] = etype + i * 1000; + memb_func[type].dparam_semantics[ix] = nrn_semantics_from_ion(etype, i); /* note that if style is needed (i==1), then we are writing a concentration */ if (i) { ion_write_depend(type, etype); @@ -300,8 +301,8 @@ int nrn_mech_depend(int type, int* dependencies) { int idep = 0; if (ds) for (int i = 0; i < dpsize; ++i) { - if (ds[i] > 0 && ds[i] < 1000) { - int deptype = ds[i]; + if (nrn_semantics_is_ion(ds[i])) { + int deptype = nrn_semantics_ion_type(ds[i]); int idepnew = depend_append(idep, dependencies, deptype, type); if ((idepnew > idep) && !corenrn.get_ion_write_dependency().empty() && !corenrn.get_ion_write_dependency()[deptype].empty()) { diff --git a/src/coreneuron/permute/node_permute.cpp b/src/coreneuron/permute/node_permute.cpp index a19230b4fc..eb33da0991 100644 --- a/src/coreneuron/permute/node_permute.cpp +++ b/src/coreneuron/permute/node_permute.cpp @@ -94,6 +94,7 @@ so pdata_m(k, isz) = inew + data_t #include "coreneuron/nrniv/nrniv_decl.h" #include "coreneuron/utils/nrn_assert.h" #include "coreneuron/coreneuron.hpp" +#include "nrnoc/ion_semantics.h" #else #include "nrnoc/multicore.h" #include "oc/nrnassrt.h" @@ -333,8 +334,8 @@ static void update_pdata_values(Memb_list* ml, int type, NrnThread& nt) { nrn_assert(0); } } - } else if (s >= 0 && s < 1000) { // ion - int etype = s; + } else if (nrn_semantics_is_ion(s)) { // ion + int etype = nrn_semantics_ion_type(s); int elayout = corenrn.get_mech_data_layout()[etype]; Memb_list* eml = nt._ml_list[etype]; int edata0 = eml->data - nt._data; diff --git a/src/neuron/cache/mechanism_range.hpp b/src/neuron/cache/mechanism_range.hpp index 4a332ebb06..8cca0ce44d 100644 --- a/src/neuron/cache/mechanism_range.hpp +++ b/src/neuron/cache/mechanism_range.hpp @@ -22,7 +22,7 @@ void indices_to_cache(short type, Callable callable) { auto const sem = dparam_semantics[field]; // See https://github.com/neuronsimulator/nrn/issues/2312 for discussion of possible // extensions to caching. - if ((sem > 0 && sem < 1000) || sem == -1 /* area */ || sem == -9 /* diam */) { + if (nrn_semantics_is_ion(sem) || sem == -1 /* area */ || sem == -9 /* diam */) { std::invoke(callable, field); } } diff --git a/src/nrniv/nrncore_write/callbacks/nrncore_callbacks.cpp b/src/nrniv/nrncore_write/callbacks/nrncore_callbacks.cpp index 158ae0c028..a91b69ec14 100644 --- a/src/nrniv/nrncore_write/callbacks/nrncore_callbacks.cpp +++ b/src/nrniv/nrncore_write/callbacks/nrncore_callbacks.cpp @@ -670,10 +670,10 @@ int* datum2int(int type, } } else if (etype == -9) { pdata[jj] = eindex; - } else if (etype > 0 && etype < 1000) { // ion pointer + } else if (nrn_semantics_is_ion(etype)) { // ion pointer pdata[jj] = eindex; - } else if (etype > 1000 && etype < 2000) { // ionstyle can be explicit instead of - // pointer to int* + } else if (nrn_semantics_is_ionstyle(etype)) { + // ionstyle can be explicit instead of pointer to int* pdata[jj] = eindex; } else if (etype == -2) { // an ion and this is the iontype pdata[jj] = eindex; diff --git a/src/nrniv/nrncore_write/data/cell_group.cpp b/src/nrniv/nrncore_write/data/cell_group.cpp index e034359387..8ee75d4191 100644 --- a/src/nrniv/nrncore_write/data/cell_group.cpp +++ b/src/nrniv/nrncore_write/data/cell_group.cpp @@ -343,15 +343,15 @@ void CellGroup::datumindex_fill(int ith, CellGroup& cg, DatumIndices& di, Memb_l } assert(etype != 0); // pointer into one of the tml types? - } else if (dmap[j] > 0 && dmap[j] < 1000) { // double* into eion type data - etype = dmap[j]; + } else if (nrn_semantics_is_ion(dmap[j])) { // double* into eion type data + etype = nrn_semantics_ion_type(dmap[j]); Memb_list* eml = cg.type2ml[etype]; assert(eml); auto* const pval = dparam[j].get(); auto const legacy_index = eml->legacy_index(pval); assert(legacy_index >= 0); eindex = legacy_index; - } else if (dmap[j] > 1000) { // int* into ion dparam[xxx][0] + } else if (nrn_semantics_is_ionstyle(dmap[j])) { // int* into ion dparam[xxx][0] // store the actual ionstyle etype = dmap[j]; eindex = *dparam[j].get(); diff --git a/src/nrniv/nrncore_write/data/cell_group.h b/src/nrniv/nrncore_write/data/cell_group.h index 374c2a2510..e86b1500cb 100644 --- a/src/nrniv/nrncore_write/data/cell_group.h +++ b/src/nrniv/nrncore_write/data/cell_group.h @@ -32,8 +32,8 @@ class CellGroup { // following three are parallel arrays std::vector output_ps; // n_presyn of these, real are first, tml order for acell. std::vector output_gid; // n_presyn of these, -(type + 1000*index) if no gid - std::vector output_vindex; // n_presyn of these. >=0 if associated with voltage, -(type + - // 1000*index) for acell. + std::vector output_vindex; // n_presyn of these. >=0 if associated with voltage, + // -(type + 1000*index) for acell. int n_netcon; // all that have targets associated with this threads Point_process. std::vector netcons; int* netcon_srcgid; // -(type + 1000*index) refers to acell with no gid diff --git a/src/nrnoc/init.cpp b/src/nrnoc/init.cpp index 267eac5c55..f77c107137 100644 --- a/src/nrnoc/init.cpp +++ b/src/nrnoc/init.cpp @@ -765,7 +765,7 @@ namespace { */ // name to int map for the negative types -// xx_ion and #xx_ion will get values of type and type+1000 respectively +// xx_ion and #xx_ion will get values of type*2 and type*2+1 respectively static std::unordered_map name_to_negint = {{"area", -1}, {"iontype", -2}, {"cvodeieq", -3}, @@ -785,7 +785,7 @@ int dparam_semantics_to_int(std::string_view name) { bool const i{name[0] == '#'}; Symbol* s = hoc_lookup(std::string{name.substr(i)}.c_str()); if (s && s->type == MECHANISM) { - return s->subtype + i * 1000; + return nrn_semantics_from_ion(s->subtype, i); } throw std::runtime_error("unknown dparam semantics: " + std::string{name}); } diff --git a/src/nrnoc/ion_semantics.h b/src/nrnoc/ion_semantics.h new file mode 100644 index 0000000000..2a3dfc8253 --- /dev/null +++ b/src/nrnoc/ion_semantics.h @@ -0,0 +1,14 @@ +#pragma once + +inline int nrn_semantics_from_ion(int type, int i) { + return 2 * type + i; +} +inline bool nrn_semantics_is_ion(int i) { + return i >= 0 && (i & 1) == 0; +} +inline bool nrn_semantics_is_ionstyle(int i) { + return i >= 0 && (i & 1) == 1; +} +inline int nrn_semantics_ion_type(int i) { + return i / 2; +} diff --git a/src/nrnoc/membfunc.h b/src/nrnoc/membfunc.h index ad80ed9d7e..18647c19b7 100644 --- a/src/nrnoc/membfunc.h +++ b/src/nrnoc/membfunc.h @@ -5,6 +5,7 @@ extern void hoc_register_prop_size(int type, int psize, int dpsize); #include "nrnoc_ml.h" #include "oc_ansi.h" // neuron::model_sorted_token #include "options.h" // EXTRACELLULAR +#include "ion_semantics.h" #include #include From d52f6ae252c966e564fc9631df315158c827a214 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 30 Oct 2024 15:12:07 +0100 Subject: [PATCH 14/17] Constify OcMatrix (#3148) Minify usage of mep as it get a non stable pointer on the matrix itself. --- src/ivoc/matrix.cpp | 17 ++-- src/ivoc/ocmatrix.cpp | 83 ++++++++++--------- src/ivoc/ocmatrix.h | 123 +++++++++++++++------------- src/nrniv/matrixmap.cpp | 16 ++-- src/nrniv/matrixmap.h | 2 +- src/nrnpython/rxd.cpp | 14 ++-- src/nrnpython/rxd_extracellular.cpp | 55 ++++++------- test/unit_tests/matrix.cpp | 5 ++ 8 files changed, 162 insertions(+), 153 deletions(-) diff --git a/src/ivoc/matrix.cpp b/src/ivoc/matrix.cpp index 8fa098d61d..ddc8d41cfa 100644 --- a/src/ivoc/matrix.cpp +++ b/src/ivoc/matrix.cpp @@ -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; } @@ -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.; } @@ -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); @@ -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); } diff --git a/src/ivoc/ocmatrix.cpp b/src/ivoc/ocmatrix.cpp index 7afa87498a..593bdef793 100644 --- a/src/ivoc/ocmatrix.cpp +++ b/src/ivoc/ocmatrix.cpp @@ -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& m, std::vector& n) { - m.clear(); - n.clear(); +std::vector> OcMatrix::nonzeros() const { + std::vector> 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() { @@ -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(); } @@ -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); } @@ -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 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 svd(m_, Eigen::ComputeFullU | Eigen::ComputeFullV); v1 = svd.singularValues(); @@ -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) { @@ -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(); } @@ -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) { @@ -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() { @@ -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; @@ -330,7 +329,7 @@ void OcSparseMatrix::setcol(int k, double in) { } } -void OcSparseMatrix::ident(void) { +void OcSparseMatrix::ident() { m_.setIdentity(); } @@ -348,7 +347,7 @@ 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; @@ -356,7 +355,7 @@ int OcSparseMatrix::sprowlen(int i) { 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) { @@ -368,13 +367,13 @@ double OcSparseMatrix::spgetrowval(int i, int jindx, int* j) { return 0; } -void OcSparseMatrix::nonzeros(std::vector& m, std::vector& n) { - m.clear(); - n.clear(); +std::vector> OcSparseMatrix::nonzeros() const { + std::vector> 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; } diff --git a/src/ivoc/ocmatrix.h b/src/ivoc/ocmatrix.h index 15960fe51f..9e874c86a3 100644 --- a/src/ivoc/ocmatrix.h +++ b/src/ivoc/ocmatrix.h @@ -1,6 +1,7 @@ #pragma once #include +#include #include #include @@ -20,23 +21,35 @@ class OcMatrix { static OcMatrix* instance(int nrow, int ncol, int type = MFULL); virtual ~OcMatrix() = default; - virtual double* mep(int i, int j) { + // This function is deprecated and should not be used! + // mep stands for 'matrix element pointer' + inline double* mep(int i, int j) { + return &coeff(i, j); + } + + inline double operator()(int i, int j) const { + return getval(i, j); + }; + + virtual double& coeff(int i, int j) { + static double zero = 0.0; unimp(); - return nullptr; - } // matrix element pointer + return zero; + } + inline double& operator()(int i, int j) { - return *mep(i, j); + return coeff(i, j); }; - virtual double getval(int i, int j) { + virtual double getval(int i, int j) const { unimp(); return 0.; } - virtual int nrow() { + virtual int nrow() const { unimp(); return 0; } - virtual int ncol() { + virtual int ncol() const { unimp(); return 0; } @@ -44,32 +57,32 @@ class OcMatrix { unimp(); } - virtual void nonzeros(std::vector& m, std::vector& n); + virtual std::vector> nonzeros() const; OcFullMatrix* full(); - inline void mulv(Vect& in, Vect& out) { + inline void mulv(Vect& in, Vect& out) const { mulv(&in, &out); }; - virtual void mulv(Vect* in, Vect* out) { + virtual void mulv(Vect* in, Vect* out) const { unimp(); } - virtual void mulm(Matrix* in, Matrix* out) { + virtual void mulm(Matrix* in, Matrix* out) const { unimp(); } - virtual void muls(double, Matrix* out) { + virtual void muls(double, Matrix* out) const { unimp(); } - virtual void add(Matrix*, Matrix* out) { + virtual void add(Matrix*, Matrix* out) const { unimp(); } - virtual void getrow(int, Vect* out) { + virtual void getrow(int, Vect* out) const { unimp(); } - virtual void getcol(int, Vect* out) { + virtual void getcol(int, Vect* out) const { unimp(); } - virtual void getdiag(int, Vect* out) { + virtual void getdiag(int, Vect* out) const { unimp(); } virtual void setrow(int, Vect* in) { @@ -96,47 +109,47 @@ class OcMatrix { virtual void ident() { unimp(); } - virtual void exp(Matrix* out) { + virtual void exp(Matrix* out) const { unimp(); } - virtual void pow(int, Matrix* out) { + virtual void pow(int, Matrix* out) const { unimp(); } - virtual void inverse(Matrix* out) { + virtual void inverse(Matrix* out) const { unimp(); } virtual void solv(Vect* vin, Vect* vout, bool use_lu) { unimp(); } - virtual void copy(Matrix* out) { + virtual void copy(Matrix* out) const { unimp(); } - virtual void bcopy(Matrix* mout, int i0, int j0, int n0, int m0, int i1, int j1) { + virtual void bcopy(Matrix* mout, int i0, int j0, int n0, int m0, int i1, int j1) const { unimp(); } virtual void transpose(Matrix* out) { unimp(); } - virtual void symmeigen(Matrix* mout, Vect* vout) { + virtual void symmeigen(Matrix* mout, Vect* vout) const { unimp(); } - virtual void svd1(Matrix* u, Matrix* v, Vect* d) { + virtual void svd1(Matrix* u, Matrix* v, Vect* d) const { unimp(); } - virtual double det(int* e) { + virtual double det(int* e) const { unimp(); return 0.0; } - virtual int sprowlen(int) { + virtual int sprowlen(int) const { unimp(); return 0; } - virtual double spgetrowval(int i, int jindx, int* j) { + virtual double spgetrowval(int i, int jindx, int* j) const { unimp(); return 0.; } - void unimp(); + void unimp() const; protected: OcMatrix(int type); @@ -155,19 +168,19 @@ class OcFullMatrix final: public OcMatrix { // type 1 OcFullMatrix(int, int); ~OcFullMatrix() override = default; - double* mep(int, int) override; - double getval(int i, int j) override; - int nrow() override; - int ncol() override; + double& coeff(int, int) override; + double getval(int i, int j) const override; + int nrow() const override; + int ncol() const 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 mulv(Vect* in, Vect* out) const override; + void mulm(Matrix* in, Matrix* out) const override; + void muls(double, Matrix* out) const override; + void add(Matrix*, Matrix* out) const override; + void getrow(int, Vect* out) const override; + void getcol(int, Vect* out) const override; + void getdiag(int, Vect* out) const override; void setrow(int, Vect* in) override; void setcol(int, Vect* in) override; void setdiag(int, Vect* in) override; @@ -176,16 +189,16 @@ class OcFullMatrix final: public OcMatrix { // type 1 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 exp(Matrix* out) const override; + void pow(int, Matrix* out) const override; + void inverse(Matrix* out) const 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 copy(Matrix* out) const override; + void bcopy(Matrix* mout, int i0, int j0, int n0, int m0, int i1, int j1) const 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; + void symmeigen(Matrix* mout, Vect* vout) const override; + void svd1(Matrix* u, Matrix* v, Vect* d) const override; + double det(int* exponent) const override; private: Eigen::Matrix m_{}; @@ -197,12 +210,12 @@ class OcSparseMatrix final: public OcMatrix { // type 2 OcSparseMatrix(int, int); ~OcSparseMatrix() override = default; - 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; + double& coeff(int, int) override; + int nrow() const override; + int ncol() const override; + double getval(int, int) const override; + void ident() override; + void mulv(Vect* in, Vect* out) const override; void solv(Vect* vin, Vect* vout, bool use_lu) override; void setrow(int, Vect* in) override; @@ -212,10 +225,10 @@ class OcSparseMatrix final: public OcMatrix { // type 2 void setcol(int, double in) override; void setdiag(int, double in) override; - void nonzeros(std::vector& m, std::vector& n) override; + std::vector> nonzeros() const override; - int sprowlen(int) override; // how many elements in row - double spgetrowval(int i, int jindx, int* j) override; + int sprowlen(int) const override; // how many elements in row + double spgetrowval(int i, int jindx, int* j) const override; void zero() override; diff --git a/src/nrniv/matrixmap.cpp b/src/nrniv/matrixmap.cpp index f07cfbed2b..658e66c5c0 100644 --- a/src/nrniv/matrixmap.cpp +++ b/src/nrniv/matrixmap.cpp @@ -19,7 +19,8 @@ void MatrixMap::mmfree() { void MatrixMap::add(double fac) { for (int i = 0; i < plen_; ++i) { - *ptree_[i] += fac * (*pm_[i]); + auto [it, jt] = pm_[i]; + *ptree_[i] += fac * m_(it, jt); } } @@ -28,13 +29,10 @@ void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) { mmfree(); plen_ = 0; - std::vector nonzero_i, nonzero_j; - m_.nonzeros(nonzero_i, nonzero_j); - pm_.resize(nonzero_i.size()); - ptree_.resize(nonzero_i.size()); - for (int k = 0; k < nonzero_i.size(); k++) { - const int i = nonzero_i[k]; - const int j = nonzero_j[k]; + std::vector> nzs = m_.nonzeros(); + pm_.resize(nzs.size()); + ptree_.resize(nzs.size()); + for (const auto [i, j]: nzs) { int it; if (i < nnode) { it = nodes[i]->eqn_index_ + layer[i]; @@ -45,7 +43,7 @@ void MatrixMap::alloc(int start, int nnode, Node** nodes, int* layer) { it = start + i - nnode; } int jt; - pm_[plen_] = m_.mep(i, j); + pm_[plen_] = std::make_pair(i, j); if (j < nnode) { jt = nodes[j]->eqn_index_ + layer[j]; if (layer[j] > 0 && !nodes[j]->extnode) { diff --git a/src/nrniv/matrixmap.h b/src/nrniv/matrixmap.h index f216dcc9c3..c20b774327 100644 --- a/src/nrniv/matrixmap.h +++ b/src/nrniv/matrixmap.h @@ -116,6 +116,6 @@ class MatrixMap { // the map int plen_ = 0; - std::vector pm_{}; + std::vector> pm_{}; std::vector ptree_{}; }; diff --git a/src/nrnpython/rxd.cpp b/src/nrnpython/rxd.cpp index 7cf61f7513..ff548a7f4a 100644 --- a/src/nrnpython/rxd.cpp +++ b/src/nrnpython/rxd.cpp @@ -1516,7 +1516,7 @@ void solve_reaction(ICSReactions* react, double pd; double dt = *dt_ptr; double dx = FLT_EPSILON; - auto jacobian = std::make_unique(N, N); + OcFullMatrix jacobian(N, N); auto b = std::make_unique(N); auto x = std::make_unique(N); @@ -1655,7 +1655,7 @@ void solve_reaction(ICSReactions* react, if (react->state_idx[segment][jac_i][jac_j] != SPECIES_ABSENT) { pd = (result_array_dx[jac_i][jac_j] - result_array[jac_i][jac_j]) / dx; - *jacobian->mep(jac_idx, idx) = (idx == jac_idx) - dt * pd; + jacobian(jac_idx, idx) = (idx == jac_idx) - dt * pd; jac_idx += 1; } result_array_dx[jac_i][jac_j] = 0; @@ -1665,7 +1665,7 @@ void solve_reaction(ICSReactions* react, // pd is our Jacobian approximated if (react->ecs_state[segment][jac_i] != NULL) { pd = (ecs_result_dx[jac_i] - ecs_result[jac_i]) / dx; - *jacobian->mep(jac_idx, idx) = -dt * pd; + jacobian(jac_idx, idx) = -dt * pd; jac_idx += 1; } ecs_result_dx[jac_i] = 0; @@ -1707,7 +1707,7 @@ void solve_reaction(ICSReactions* react, // pd is our Jacobian approximated if (react->state_idx[segment][jac_i][jac_j] != SPECIES_ABSENT) { pd = (result_array_dx[jac_i][jac_j] - result_array[jac_i][jac_j]) / dx; - *jacobian->mep(jac_idx, idx) = -dt * pd; + jacobian(jac_idx, idx) = -dt * pd; jac_idx += 1; } } @@ -1716,10 +1716,10 @@ void solve_reaction(ICSReactions* react, // pd is our Jacobian approximated if (react->ecs_state[segment][jac_i] != NULL) { pd = (ecs_result_dx[jac_i] - ecs_result[jac_i]) / dx; - *jacobian->mep(jac_idx, idx) = (idx == jac_idx) - dt * pd; + jacobian(jac_idx, idx) = (idx == jac_idx) - dt * pd; jac_idx += 1; } else { - *jacobian->mep(idx, idx) = 1.0; + jacobian(idx, idx) = 1.0; } // reset dx array ecs_states_for_reaction_dx[i] -= dx; @@ -1728,7 +1728,7 @@ void solve_reaction(ICSReactions* react, } } // solve for x, destructively - jacobian->solv(b.get(), x.get(), false); + jacobian.solv(b.get(), x.get(), false); if (bval != NULL) // variable-step { diff --git a/src/nrnpython/rxd_extracellular.cpp b/src/nrnpython/rxd_extracellular.cpp index aeb377eed4..a26d61eeb2 100644 --- a/src/nrnpython/rxd_extracellular.cpp +++ b/src/nrnpython/rxd_extracellular.cpp @@ -290,7 +290,6 @@ void* ecs_do_reactions(void* dataptr) { double* mc_mults_array = NULL; double dx = FLT_EPSILON; double pd; - std::unique_ptr jacobian; std::vector x{}; std::vector b{}; @@ -316,8 +315,7 @@ void* ecs_do_reactions(void* dataptr) { if (react->num_species_involved == 0) continue; /*allocate data structures*/ - jacobian = std::make_unique(react->num_species_involved, - react->num_species_involved); + OcFullMatrix jacobian(react->num_species_involved, react->num_species_involved); b.resize(react->num_species_involved); x.resize(react->num_species_involved); states_cache = (double*) malloc(sizeof(double) * react->num_species_involved); @@ -358,23 +356,23 @@ void* ecs_do_reactions(void* dataptr) { for (k = 0; k < react->num_species_involved; k++) { pd = (results_array_dx[k] - results_array[k]) / dx; - *jacobian->mep(k, j) = (j == k) - dt * pd; + jacobian(k, j) = (j == k) - dt * pd; } states_cache_dx[j] -= dx; } // solve for x if (react->num_species_involved == 1) { - react->species_states[0][i] += b[0] / jacobian->getval(0, 0); + react->species_states[0][i] += b[0] / jacobian(0, 0); } else { // find entry in leftmost column with largest absolute value // Pivot for (j = 0; j < react->num_species_involved; j++) { for (k = j + 1; k < react->num_species_involved; k++) { - if (abs(jacobian->getval(j, j)) < abs(jacobian->getval(k, j))) { + if (abs(jacobian(j, j)) < abs(jacobian(k, j))) { for (n = 0; n < react->num_species_involved; n++) { - temp = jacobian->getval(j, n); - *jacobian->mep(j, n) = jacobian->getval(k, n); - *jacobian->mep(k, n) = temp; + temp = jacobian(j, n); + jacobian(j, n) = jacobian(k, n); + jacobian(k, n) = temp; } } } @@ -382,11 +380,11 @@ void* ecs_do_reactions(void* dataptr) { for (j = 0; j < react->num_species_involved - 1; j++) { for (k = j + 1; k < react->num_species_involved; k++) { - ge_value = jacobian->getval(k, j) / jacobian->getval(j, j); + ge_value = jacobian(k, j) / jacobian(j, j); for (n = 0; n < react->num_species_involved; n++) { - val_to_set = jacobian->getval(k, n) - - ge_value * jacobian->getval(j, n); - *jacobian->mep(k, n) = val_to_set; + val_to_set = jacobian(k, n) - + ge_value * jacobian(j, n); + jacobian(k, n) = val_to_set; } b[k] = b[k] - ge_value * b[j]; } @@ -396,10 +394,10 @@ void* ecs_do_reactions(void* dataptr) { x[j] = b[j]; for (k = j + 1; k < react->num_species_involved; k++) { if (k != j) { - x[j] = x[j] - jacobian->getval(j, k) * x[k]; + x[j] = x[j] - jacobian(j, k) * x[k]; } } - x[j] = x[j] / jacobian->getval(j, j); + x[j] = x[j] / jacobian(j, j); } for (j = 0; j < react->num_species_involved; j++) { // I think this should be something like @@ -446,8 +444,7 @@ void* ecs_do_reactions(void* dataptr) { if (react->num_species_involved == 0) continue; /*allocate data structures*/ - jacobian = std::make_unique(react->num_species_involved, - react->num_species_involved); + OcFullMatrix jacobian(react->num_species_involved, react->num_species_involved); b.resize(react->num_species_involved); x.resize(react->num_species_involved); states_cache = (double*) malloc(sizeof(double) * react->num_species_involved); @@ -480,23 +477,23 @@ void* ecs_do_reactions(void* dataptr) { for (k = 0; k < react->num_species_involved; k++) { pd = (results_array_dx[k] - results_array[k]) / dx; - *jacobian->mep(k, j) = (j == k) - dt * pd; + jacobian(k, j) = (j == k) - dt * pd; } states_cache_dx[j] -= dx; } // solve for x if (react->num_species_involved == 1) { - react->species_states[0][i] += b[0] / jacobian->getval(0, 0); + react->species_states[0][i] += b[0] / jacobian(0, 0); } else { // find entry in leftmost column with largest absolute value // Pivot for (j = 0; j < react->num_species_involved; j++) { for (k = j + 1; k < react->num_species_involved; k++) { - if (abs(jacobian->getval(j, j)) < abs(jacobian->getval(k, j))) { + if (abs(jacobian(j, j)) < abs(jacobian(k, j))) { for (n = 0; n < react->num_species_involved; n++) { - temp = jacobian->getval(j, n); - *jacobian->mep(j, n) = jacobian->getval(k, n); - *jacobian->mep(k, n) = temp; + temp = jacobian(j, n); + jacobian(j, n) = jacobian(k, n); + jacobian(k, n) = temp; } } } @@ -504,11 +501,11 @@ void* ecs_do_reactions(void* dataptr) { for (j = 0; j < react->num_species_involved - 1; j++) { for (k = j + 1; k < react->num_species_involved; k++) { - ge_value = jacobian->getval(k, j) / jacobian->getval(j, j); + ge_value = jacobian(k, j) / jacobian(j, j); for (n = 0; n < react->num_species_involved; n++) { - val_to_set = jacobian->getval(k, n) - - ge_value * jacobian->getval(j, n); - *jacobian->mep(k, n) = val_to_set; + val_to_set = jacobian(k, n) - + ge_value * jacobian(j, n); + jacobian(k, n) = val_to_set; } b[k] = b[k] - ge_value * b[j]; } @@ -518,10 +515,10 @@ void* ecs_do_reactions(void* dataptr) { x[j] = b[j]; for (k = j + 1; k < react->num_species_involved; k++) { if (k != j) { - x[j] = x[j] - jacobian->getval(j, k) * x[k]; + x[j] = x[j] - jacobian(j, k) * x[k]; } } - x[j] = x[j] / jacobian->getval(j, j); + x[j] = x[j] / jacobian(j, j); } for (j = 0; j < react->num_species_involved; j++) { // I think this should be something like diff --git a/test/unit_tests/matrix.cpp b/test/unit_tests/matrix.cpp index 27386501ca..115dbfe120 100644 --- a/test/unit_tests/matrix.cpp +++ b/test/unit_tests/matrix.cpp @@ -94,6 +94,11 @@ SCENARIO("A Matrix", "[neuron_ivoc][OcMatrix]") { m.setrow(3, 2.0); REQUIRE(compareMatrix(m, {{3., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}, {2., 2., 2.}})); } + { + std::vector> nzs = m.nonzeros(); + std::vector> res = {{0, 0}, {1, 1} , {2, 2}, {3, 0}, {3, 1}, {3, 2}}; + REQUIRE(nzs == res); + } { std::vector x, y; m.nonzeros(x, y); From a6510feaff8c66dd611fce6b783adda6cb15ff29 Mon Sep 17 00:00:00 2001 From: nrnhines Date: Wed, 30 Oct 2024 19:21:58 -0400 Subject: [PATCH 15/17] Allow parallel test coverage (e.g. ctest -j 8 ) (#3162) --- cmake/Coverage.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/Coverage.cmake b/cmake/Coverage.cmake index aa2bae3bf5..9178c17636 100644 --- a/cmake/Coverage.cmake +++ b/cmake/Coverage.cmake @@ -34,7 +34,7 @@ if(NRN_ENABLE_COVERAGE) if(NOT BUILD_TYPE_UPPER STREQUAL "DEBUG") message(WARNING "Using CMAKE_BUILD_TYPE=Debug is recommended with NRN_ENABLE_COVERAGE") endif() - set(NRN_COVERAGE_FLAGS_UNQUOTED --coverage -fno-inline) + set(NRN_COVERAGE_FLAGS_UNQUOTED --coverage -fno-inline -fprofile-update=atomic) string(JOIN " " NRN_COVERAGE_FLAGS ${NRN_COVERAGE_FLAGS_UNQUOTED}) set(NRN_COVERAGE_LINK_FLAGS --coverage) From ebe3d3ba9a22beafe6f9b6d1099977f67d24f19d Mon Sep 17 00:00:00 2001 From: Luc Grosheintz Date: Thu, 31 Oct 2024 11:16:58 +0100 Subject: [PATCH 16/17] Update `external/nmodl`. (#3151) --- external/nmodl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/nmodl b/external/nmodl index e6250014dd..800d098cc3 160000 --- a/external/nmodl +++ b/external/nmodl @@ -1 +1 @@ -Subproject commit e6250014dd2c5ae5d9cb30af6920668835b31137 +Subproject commit 800d098cc3bf658cf7eed314df5870aaf96d9dd8 From 7653d2154fc33c90a74e048bcef42b25c6d2155b Mon Sep 17 00:00:00 2001 From: Emmanuel Ferdman Date: Fri, 1 Nov 2024 20:06:21 +0200 Subject: [PATCH 17/17] Update Python MacOS download page (#3169) Signed-off-by: Emmanuel Ferdman --- docs/install/mac_pkg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/install/mac_pkg.md b/docs/install/mac_pkg.md index 86cee22d88..cc93bd80c5 100644 --- a/docs/install/mac_pkg.md +++ b/docs/install/mac_pkg.md @@ -162,7 +162,7 @@ generally only for the architecture indicated by ```uname -m```. That is ok for openmpi but since the various python libraries are linked against during build to create the version specific neuron modules, those python installers also have to be universal. Fortunately, universal python versions -can be found at [python.org](http://python.org/Downloads/macOS) at least for +can be found at [python.org](https://www.python.org/downloads/macos) at least for (as of 2022-01-01) python3.8, python3.9, and python3.10. - ```xcode-select --install```: