diff --git a/.circleci/config.yml b/.circleci/config.yml index c2d9b65aa7..6ec6f56ba6 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -1,7 +1,7 @@ version: 2.1 orbs: - python: circleci/python@0.3.2 + python: circleci/python@2.1.1 jobs: manylinux2014-aarch64: @@ -13,7 +13,7 @@ jobs: type: string machine: - image: ubuntu-2004:202101-01 + image: default resource_class: arm.medium @@ -49,25 +49,18 @@ jobs: # choose available python versions from pyenv pyenv_py_ver="" case << parameters.NRN_PYTHON_VERSION >> in - 38) pyenv_py_ver="3.8.7" ;; - 39) pyenv_py_ver="3.9.1" ;; - 310) pyenv_py_ver="3.10.1" ;; - 311) pyenv_py_ver="3.11.0" ;; + 38) pyenv_py_ver="3.8" ;; + 39) pyenv_py_ver="3.9" ;; + 310) pyenv_py_ver="3.10" ;; + 311) pyenv_py_ver="3.11" ;; + 312) pyenv_py_ver="3.12" ;; *) echo "Error: pyenv python version not specified!" && exit 1;; esac - # install python dependencies: .10 is not available pyenv - if [ "<< parameters.NRN_PYTHON_VERSION >>" == "310" ]; then - sudo apt install software-properties-common -y - sudo add-apt-repository ppa:deadsnakes/ppa -y - sudo apt install python3.10 libpython3.10 python3.10-venv - export PYTHON_EXE=$(which python3.10) - else - cd /opt/circleci/.pyenv/plugins/python-build/../.. && git pull && cd - - env PYTHON_CONFIGURE_OPTS="--enable-shared" pyenv install $pyenv_py_ver --force - pyenv global $pyenv_py_ver - export PYTHON_EXE=$(which python) - fi + cd /opt/circleci/.pyenv/plugins/python-build/../.. && git fetch --all && git checkout -B master origin/master && cd - + env PYTHON_CONFIGURE_OPTS="--enable-shared" pyenv install $pyenv_py_ver --force + pyenv global $pyenv_py_ver + export PYTHON_EXE=$(which python) # test wheel packaging/python/test_wheels.sh $PYTHON_EXE $(ls -t wheelhouse/*.whl) @@ -95,7 +88,7 @@ workflows: - /circleci\/.*/ matrix: parameters: - NRN_PYTHON_VERSION: ["311"] + NRN_PYTHON_VERSION: ["312"] NRN_NIGHTLY_UPLOAD: ["false"] nightly: @@ -110,5 +103,5 @@ workflows: - manylinux2014-aarch64: matrix: parameters: - NRN_PYTHON_VERSION: ["38", "39", "310", "311"] + NRN_PYTHON_VERSION: ["38", "39", "310", "311", "312"] NRN_NIGHTLY_UPLOAD: ["true"] diff --git a/.github/workflows/neuron-ci.yml b/.github/workflows/neuron-ci.yml index 9472e60686..0d0d05cc46 100644 --- a/.github/workflows/neuron-ci.yml +++ b/.github/workflows/neuron-ci.yml @@ -40,7 +40,7 @@ jobs: PY_MIN_VERSION: ${{ matrix.config.python_min_version || '3.8' }} PY_MAX_VERSION: ${{ matrix.config.python_max_version || '3.11' }} MUSIC_INSTALL_DIR: /opt/MUSIC - MUSIC_VERSION: 1.2.0 + MUSIC_VERSION: 1.2.1 strategy: matrix: @@ -173,7 +173,8 @@ jobs: curl -L -o MUSIC.zip https://github.com/INCF/MUSIC/archive/refs/tags/${MUSIC_VERSION}.zip unzip MUSIC.zip && mv MUSIC-* MUSIC && cd MUSIC ./autogen.sh - ./configure --with-python-sys-prefix --prefix=$MUSIC_INSTALL_DIR --disable-anysource + # on some systems MPI library detection fails, provide exact flags/compilers + ./configure --with-python-sys-prefix --prefix=$MUSIC_INSTALL_DIR --disable-anysource MPI_CXXFLAGS="-g -O3" MPI_CFLAGS="-g -O3" MPI_LDFLAGS=" " CC=mpicc CXX=mpicxx make -j install deactivate working-directory: ${{runner.temp}} diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index d6baa64fb0..bd944ee463 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -53,6 +53,7 @@ mac_m1_cmake_build: - real_python=$(python3 resolve_shim.py) - echo "python3=$(command -v python3) is really ${real_python}" - PYTHONEXECUTABLE=${real_python} ${real_python} -mvenv venv + - venv/bin/python3 -m ensurepip --upgrade --default-pip - venv/bin/pip install --upgrade pip -r nrn_requirements.txt - git submodule update --init --recursive --force --depth 1 -- external/nmodl - venv/bin/pip install --upgrade -r external/nmodl/requirements.txt diff --git a/README.md b/README.md index 96d5b467f2..f6cba70ead 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -[![Build Status](https://dev.azure.com/neuronsimulator/nrn/_apis/build/status/neuronsimulator.nrn?branchName=master)](https://dev.azure.com/neuronsimulator/nrn/_build/latest?definitionId=1&branchName=master) [![Actions Status](https://github.com/neuronsimulator/nrn/workflows/Windows%20Installer/badge.svg)](https://github.com/neuronsimulator/nrn/actions) [![Actions Status](https://github.com/neuronsimulator/nrn/workflows/NEURON%20CI/badge.svg)](https://github.com/neuronsimulator/nrn/actions) [![codecov](https://codecov.io/gh/neuronsimulator/nrn/branch/master/graph/badge.svg?token=T7PIDw6LrC)](https://codecov.io/gh/neuronsimulator/nrn) [![Documentation Status](https://readthedocs.org/projects/nrn/badge/?version=latest)](http://nrn.readthedocs.io/?badge=latest) +[![Build Status](https://dev.azure.com/neuronsimulator/nrn/_apis/build/status/neuronsimulator.nrn?branchName=master)](https://dev.azure.com/neuronsimulator/nrn/_build/latest?definitionId=1&branchName=master) [![Actions Status](https://github.com/neuronsimulator/nrn/actions/workflows/windows.yml/badge.svg?branch=master)](https://github.com/neuronsimulator/nrn/actions) [![Actions Status](https://github.com/neuronsimulator/nrn/workflows/NEURON%20CI/badge.svg)](https://github.com/neuronsimulator/nrn/actions) [![codecov](https://codecov.io/gh/neuronsimulator/nrn/branch/master/graph/badge.svg?token=T7PIDw6LrC)](https://codecov.io/gh/neuronsimulator/nrn) [![Documentation Status](https://readthedocs.org/projects/nrn/badge/?version=latest)](http://nrn.readthedocs.io/?badge=latest) # NEURON NEURON is a simulator for models of neurons and networks of neuron. See [http://neuron.yale.edu](http://neuron.yale.edu) for installers, source code, documentation, tutorials, announcements of diff --git a/ci/win_build_cmake.sh b/ci/win_build_cmake.sh index 3bf55bc842..c3c84ad37e 100755 --- a/ci/win_build_cmake.sh +++ b/ci/win_build_cmake.sh @@ -32,7 +32,7 @@ cd $BUILD_SOURCESDIRECTORY/build -DNRN_BINARY_DIST_BUILD=ON \ -DPYTHON_EXECUTABLE=/c/Python38/python.exe \ -DNRN_ENABLE_PYTHON_DYNAMIC=ON \ - -DNRN_PYTHON_DYNAMIC='c:/Python38/python.exe;c:/Python39/python.exe;c:/Python310/python.exe;c:/Python311/python.exe' \ + -DNRN_PYTHON_DYNAMIC='c:/Python38/python.exe;c:/Python39/python.exe;c:/Python310/python.exe;c:/Python311/python.exe;c:/Python312/python.exe' \ -DCMAKE_INSTALL_PREFIX='/c/nrn-install' \ -DMPI_CXX_LIB_NAMES:STRING=msmpi \ -DMPI_C_LIB_NAMES:STRING=msmpi \ diff --git a/ci/win_download_deps.cmd b/ci/win_download_deps.cmd index fc58ed4254..a006ffc0d2 100644 --- a/ci/win_download_deps.cmd +++ b/ci/win_download_deps.cmd @@ -7,6 +7,7 @@ pwsh -command Invoke-WebRequest -MaximumRetryCount 4 -OutFile python-3.8.exe htt pwsh -command Invoke-WebRequest -MaximumRetryCount 4 -OutFile python-3.9.exe https://www.python.org/ftp/python/3.9.0/python-3.9.0-amd64.exe || goto :error pwsh -command Invoke-WebRequest -MaximumRetryCount 4 -OutFile python-3.10.exe https://www.python.org/ftp/python/3.10.0/python-3.10.0-amd64.exe || goto :error pwsh -command Invoke-WebRequest -MaximumRetryCount 4 -OutFile python-3.11.exe https://www.python.org/ftp/python/3.11.1/python-3.11.1-amd64.exe || goto :error +pwsh -command Invoke-WebRequest -MaximumRetryCount 4 -OutFile python-3.12.exe https://www.python.org/ftp/python/3.12.1/python-3.12.1-amd64.exe || goto :error :: mpi pwsh -command Invoke-WebRequest -MaximumRetryCount 4 -OutFile msmpisetup.exe https://download.microsoft.com/download/a/5/2/a5207ca5-1203-491a-8fb8-906fd68ae623/msmpisetup.exe || goto :error diff --git a/ci/win_install_deps.cmd b/ci/win_install_deps.cmd index 189287b966..2284b911c0 100644 --- a/ci/win_install_deps.cmd +++ b/ci/win_install_deps.cmd @@ -7,6 +7,7 @@ python-3.8.exe /passive Include_pip=1 Include_test=0 PrependPath=1 DefaultJustFo python-3.9.exe /passive Include_pip=1 Include_test=0 PrependPath=1 DefaultJustForMeTargetDir=C:\Python39 || goto :error python-3.10.exe /passive Include_pip=1 Include_test=0 PrependPath=1 DefaultJustForMeTargetDir=C:\Python310 || goto :error python-3.11.exe /passive Include_pip=1 Include_test=0 PrependPath=1 DefaultJustForMeTargetDir=C:\Python311 || goto :error +python-3.12.exe /passive Include_pip=1 Include_test=0 PrependPath=1 DefaultJustForMeTargetDir=C:\Python312 || goto :error :: fix msvcc version for all python3 pwsh -command "(Get-Content C:\Python38\Lib\distutils\cygwinccompiler.py) -replace 'elif msc_ver == ''1600'':', 'elif msc_ver == ''1916'':' | Out-File C:\Python38\Lib\distutils\cygwinccompiler.py" @@ -26,6 +27,8 @@ C:\Python38\python.exe -m pip install numpy==1.17.5 "cython < 3" || goto :error C:\Python39\python.exe -m pip install numpy==1.19.3 "cython < 3" || goto :error C:\Python310\python.exe -m pip install numpy==1.21.3 "cython < 3" || goto :error C:\Python311\python.exe -m pip install numpy==1.23.5 "cython < 3" || goto :error +C:\Python312\python.exe -m pip install numpy==1.26.3 "cython < 3" || goto :error +C:\Python312\python.exe -m pip install setuptools || goto :error :: install nsis nsis-3.05-setup.exe /S || goto :error diff --git a/cmake/CompilerHelper.cmake b/cmake/CompilerHelper.cmake index 54adf1ce31..e50b120351 100644 --- a/cmake/CompilerHelper.cmake +++ b/cmake/CompilerHelper.cmake @@ -51,7 +51,7 @@ if(CMAKE_C_COMPILER_ID MATCHES "PGI" OR CMAKE_C_COMPILER_ID MATCHES "NVHPC") # Random123 does not play nicely with NVHPC 21.11+'s detection of ABM features, see: # https://github.com/BlueBrain/CoreNeuron/issues/724 and # https://github.com/DEShawResearch/random123/issues/6. - list(APPEND NRN_COMPILE_DEFS R123_USE_INTRIN_H=0) + list(APPEND NRN_R123_COMPILE_DEFS R123_USE_INTRIN_H=0) endif() else() set(NRN_HAVE_NVHPC_COMPILER OFF) diff --git a/cmake/NeuronFileLists.cmake b/cmake/NeuronFileLists.cmake index 8cd9a06553..c1beffb9fd 100644 --- a/cmake/NeuronFileLists.cmake +++ b/cmake/NeuronFileLists.cmake @@ -9,6 +9,9 @@ set(STRUCTURED_HEADER_FILES_TO_INSTALL neuron/container/generic_data_handle.hpp neuron/container/non_owning_soa_identifier.hpp neuron/model_data_fwd.hpp) set(HEADER_FILES_TO_INSTALL + gnu/mcran4.h + gnu/nrnisaac.h + gnu/nrnran123.h nrniv/backtrace_utils.h nrniv/bbsavestate.h nrnmpi/nrnmpidec.h @@ -37,14 +40,12 @@ set(HEADER_FILES_TO_INSTALL oc/hocgetsym.h oc/hoclist.h oc/hocparse.h - oc/mcran4.h oc/mech_api.h + oc/memory.hpp oc/nrnapi.h oc/nrnassrt.h - oc/nrnisaac.h oc/nrnmpi.h oc/nrnrandom.h - oc/nrnran123.h oc/oc_ansi.h oc/ocfunc.h oc/ocmisc.h @@ -68,8 +69,6 @@ set(HEADER_FILES_TO_INSTALL scopmath/sparse_thread.hpp scopmath/ssimplic.hpp scopmath/ssimplic_thread.hpp - sparse13/cspmatrix.h - sparse13/cspredef.h sparse13/spconfig.h sparse13/spmatrix.h) @@ -78,16 +77,10 @@ set(HEADER_FILES_TO_INSTALL # ============================================================================= set(NRN_HEADERS_INCLUDE_LIST) -# ============================================================================= -# Lists of random number related files -# ============================================================================= -set(RAN_FILE_LIST isaac64.cpp mcran4.cpp nrnisaac.cpp nrnran123.cpp) - # ============================================================================= # Files in oc directory # ============================================================================= set(OC_FILE_LIST - ${RAN_FILE_LIST} audit.cpp axis.cpp code.cpp @@ -103,6 +96,8 @@ set(OC_FILE_LIST hoc_oop.cpp list.cpp math.cpp + oc_mcran4.cpp + memory.cpp mswinprt.cpp nonlin.cpp ocerf.cpp @@ -222,7 +217,6 @@ set(NRNIV_FILE_LIST cxprop.cpp datapath.cpp finithnd.cpp - geometry3d.cpp glinerec.cpp hocmech.cpp impedanc.cpp diff --git a/nrn_requirements.txt b/nrn_requirements.txt index 06c1f1920a..77ccb32ddc 100644 --- a/nrn_requirements.txt +++ b/nrn_requirements.txt @@ -13,6 +13,3 @@ pytest-cov mpi4py numpy find_libpython - -# For python 3.12 -setuptools diff --git a/share/lib/python/neuron/__init__.py b/share/lib/python/neuron/__init__.py index df12031798..c48fdd56ef 100644 --- a/share/lib/python/neuron/__init__.py +++ b/share/lib/python/neuron/__init__.py @@ -1217,7 +1217,7 @@ def _get_color(variable, val, cmap, lo, hi, val_range): elif val > hi: col = color_to_hex(cmap(255)) else: - val = color_to_hex(128) + col = color_to_hex(cmap(128)) else: col = color_to_hex( cmap(int(255 * (min(max(val, lo), hi) - lo) / (val_range))) diff --git a/src/coreneuron/mpi/nrnmpiuse.h b/src/coreneuron/mpi/nrnmpiuse.h index 72bcd90092..eec7609094 100644 --- a/src/coreneuron/mpi/nrnmpiuse.h +++ b/src/coreneuron/mpi/nrnmpiuse.h @@ -24,6 +24,3 @@ /* define to the dll path if you want to load automatically */ #undef DLL_DEFAULT_FNAME - -/* Number of times to retry a failed open */ -#undef FILE_OPEN_RETRY diff --git a/src/coreneuron/utils/randoms/nrnran123.cpp b/src/coreneuron/utils/randoms/nrnran123.cpp index a618312a34..2f1b12cb3d 100644 --- a/src/coreneuron/utils/randoms/nrnran123.cpp +++ b/src/coreneuron/utils/randoms/nrnran123.cpp @@ -88,7 +88,7 @@ namespace random123_global { #else #define g_k_qualifiers #endif -g_k_qualifiers philox4x32_key_t g_k{}; +g_k_qualifiers philox4x32_key_t g_k{{0, 0}}; // Cannot refer to g_k directly from a nrn_pragma_acc(routine seq) method like // coreneuron_random123_philox4x32_helper, and cannot have this inlined there at diff --git a/src/gnu/CMakeLists.txt b/src/gnu/CMakeLists.txt index 0e9562297c..e16ae608b9 100644 --- a/src/gnu/CMakeLists.txt +++ b/src/gnu/CMakeLists.txt @@ -6,14 +6,23 @@ add_library( Erlang.cpp Geom.cpp HypGeom.cpp + isaac64.cpp + Isaac64RNG.cpp LogNorm.cpp + MCellRan4RNG.cpp + mcran4.cpp MLCG.cpp NegExp.cpp Normal.cpp + nrnisaac.cpp + nrnran123.cpp Poisson.cpp + Rand.cpp Random.cpp RndInt.cpp RNG.cpp Uniform.cpp Weibull.cpp) set_property(TARGET nrngnu PROPERTY POSITION_INDEPENDENT_CODE ON) +target_include_directories(nrngnu SYSTEM PRIVATE "${PROJECT_SOURCE_DIR}/external/Random123/include") +target_compile_definitions(nrngnu PRIVATE ${NRN_R123_COMPILE_DEFS}) diff --git a/src/gnu/Isaac64RNG.cpp b/src/gnu/Isaac64RNG.cpp new file mode 100644 index 0000000000..559c307d4d --- /dev/null +++ b/src/gnu/Isaac64RNG.cpp @@ -0,0 +1,20 @@ +#include "Isaac64RNG.hpp" + +uint32_t Isaac64::cnt_ = 0; + +Isaac64::Isaac64(std::uint32_t seed) { + if (cnt_ == 0) { + cnt_ = 0xffffffff; + } + --cnt_; + seed_ = seed; + if (seed_ == 0) { + seed_ = cnt_; + } + rng_ = nrnisaac_new(); + reset(); +} + +Isaac64::~Isaac64() { + nrnisaac_delete(rng_); +} diff --git a/src/gnu/Isaac64RNG.hpp b/src/gnu/Isaac64RNG.hpp new file mode 100644 index 0000000000..c825c741ca --- /dev/null +++ b/src/gnu/Isaac64RNG.hpp @@ -0,0 +1,33 @@ +#pragma once + +#include + +#include "RNG.h" +#include "nrnisaac.h" + +class Isaac64: public RNG { + public: + Isaac64(std::uint32_t seed = 0); + ~Isaac64(); + std::uint32_t asLong() { + return nrnisaac_uint32_pick(rng_); + } + void reset() { + nrnisaac_init(rng_, seed_); + } + double asDouble() { + return nrnisaac_dbl_pick(rng_); + } + std::uint32_t seed() { + return seed_; + } + void seed(std::uint32_t s) { + seed_ = s; + reset(); + } + + private: + std::uint32_t seed_; + void* rng_; + static std::uint32_t cnt_; +}; diff --git a/src/gnu/MCellRan4RNG.cpp b/src/gnu/MCellRan4RNG.cpp new file mode 100644 index 0000000000..43b732cbe1 --- /dev/null +++ b/src/gnu/MCellRan4RNG.cpp @@ -0,0 +1,15 @@ +#include "MCellRan4RNG.hpp" + +MCellRan4::MCellRan4(std::uint32_t ihigh, std::uint32_t ilow) { + ++cnt_; + ilow_ = ilow; + ihigh_ = ihigh; + if (ihigh_ == 0) { + ihigh_ = cnt_; + ihigh_ = (std::uint32_t) asLong(); + } + orig_ = ihigh_; +} +MCellRan4::~MCellRan4() {} + +std::uint32_t MCellRan4::cnt_ = 0; diff --git a/src/gnu/MCellRan4RNG.hpp b/src/gnu/MCellRan4RNG.hpp new file mode 100644 index 0000000000..a4356562a9 --- /dev/null +++ b/src/gnu/MCellRan4RNG.hpp @@ -0,0 +1,34 @@ +#pragma once + +#include + +#include "RNG.h" +#include "mcran4.h" + +// The decision that has to be made is whether each generator instance +// should have its own seed or only one seed for all. We choose separate +// seed for each but if arg not present or 0 then seed chosen by system. + +// the addition of ilow > 0 means that value is used for the lowindex +// instead of the mcell_ran4_init global 32 bit lowindex. + +class MCellRan4: public RNG { + public: + MCellRan4(std::uint32_t ihigh = 0, std::uint32_t ilow = 0); + virtual ~MCellRan4(); + virtual std::uint32_t asLong() { + return (std::uint32_t) (ilow_ == 0 ? mcell_iran4(&ihigh_) : nrnRan4int(&ihigh_, ilow_)); + } + virtual void reset() { + ihigh_ = orig_; + } + virtual double asDouble() { + return (ilow_ == 0 ? mcell_ran4a(&ihigh_) : nrnRan4dbl(&ihigh_, ilow_)); + } + std::uint32_t ihigh_; + std::uint32_t orig_; + std::uint32_t ilow_; + + private: + static std::uint32_t cnt_; +}; diff --git a/src/gnu/NrnRandom123RNG.cpp b/src/gnu/NrnRandom123RNG.cpp new file mode 100644 index 0000000000..e69de29bb2 diff --git a/src/gnu/NrnRandom123RNG.hpp b/src/gnu/NrnRandom123RNG.hpp new file mode 100644 index 0000000000..5e790c939c --- /dev/null +++ b/src/gnu/NrnRandom123RNG.hpp @@ -0,0 +1,28 @@ +#pragma once + +#include + +#include "nrnran123.h" +#include "RNG.h" + +class NrnRandom123: public RNG { + public: + NrnRandom123(std::uint32_t id1, std::uint32_t id2, std::uint32_t id3 = 0); + ~NrnRandom123(); + std::uint32_t asLong() { + return nrnran123_ipick(s_); + } + double asDouble() { + return nrnran123_dblpick(s_); + } + void reset() { + nrnran123_setseq(s_, 0, 0); + } + nrnran123_State* s_; +}; +NrnRandom123::NrnRandom123(std::uint32_t id1, std::uint32_t id2, std::uint32_t id3) { + s_ = nrnran123_newstream3(id1, id2, id3); +} +NrnRandom123::~NrnRandom123() { + nrnran123_deletestream(s_); +} diff --git a/src/gnu/Rand.cpp b/src/gnu/Rand.cpp new file mode 100644 index 0000000000..a3b526de5e --- /dev/null +++ b/src/gnu/Rand.cpp @@ -0,0 +1,19 @@ +#include "Rand.hpp" + +#include "ACG.h" +#include "Normal.h" + +Rand::Rand(unsigned long seed, int size, Object* obj) { + // printf("Rand\n"); + gen = new ACG(seed, size); + rand = new Normal(0., 1., gen); + type_ = 0; + obj_ = obj; +} + +Rand::~Rand() { + // printf("~Rand\n"); + delete gen; + delete rand; +} + diff --git a/src/ivoc/random1.h b/src/gnu/Rand.hpp similarity index 63% rename from src/ivoc/random1.h rename to src/gnu/Rand.hpp index 8fedd8f5b0..f5b47a6a61 100644 --- a/src/ivoc/random1.h +++ b/src/gnu/Rand.hpp @@ -1,14 +1,20 @@ -#ifndef random1_h -#define random1_h +#pragma once #include "RNG.h" #include "Random.h" struct Object; +/* type_: + * 0: ACG + * 1: MLCG + * 2: MCellRan4 + * 3: Isaac64 + * 4: Random123 + */ class Rand { public: - Rand(unsigned long seed = 0, int size = 55, Object* obj = NULL); + Rand(unsigned long seed = 0, int size = 55, Object* obj = nullptr); ~Rand(); RNG* gen; Random* rand; @@ -16,5 +22,3 @@ class Rand { // double* looks like random variable that gets changed on every fadvance Object* obj_; }; - -#endif diff --git a/src/oc/isaac64.cpp b/src/gnu/isaac64.cpp similarity index 98% rename from src/oc/isaac64.cpp rename to src/gnu/isaac64.cpp index 579f52d4f0..b0c6c40317 100644 --- a/src/oc/isaac64.cpp +++ b/src/gnu/isaac64.cpp @@ -1,8 +1,3 @@ -#include <../../nrnconf.h> -#if HAVE_SYS_TYPES_H -#include -#endif - /* ------------------------------------------------------------------------------ isaac64.cpp: A Fast cryptographic random number generator diff --git a/src/oc/isaac64.h b/src/gnu/isaac64.h similarity index 97% rename from src/oc/isaac64.h rename to src/gnu/isaac64.h index 6b13dd72f1..4af3b7586e 100644 --- a/src/oc/isaac64.h +++ b/src/gnu/isaac64.h @@ -1,3 +1,5 @@ +#pragma once + /* ------------------------------------------------------------------------------ isaac64.h: Definitions for a fast cryptographic random number generator @@ -12,10 +14,7 @@ Jenkins, R.J. (1996) ISAAC, in Fast Software Encryption, vol. 1039, ------------------------------------------------------------------------------ */ -#ifndef ISAAC64_H -#define ISAAC64_H -#include #include #define RANDSIZL (4) /* I recommend 8 for crypto, 4 for simulations */ @@ -87,5 +86,3 @@ Macros to get individual random numbers rng->randcnt = RANDMAX - 2, \ DBL64 * (*((ub8*) (((ub4*) (rng->randrsl)) + rng->randcnt))))) - -#endif /* ISAAC64_H */ diff --git a/src/oc/mcran4.cpp b/src/gnu/mcran4.cpp similarity index 86% rename from src/oc/mcran4.cpp rename to src/gnu/mcran4.cpp index 1b303164d6..a9553687f6 100644 --- a/src/oc/mcran4.cpp +++ b/src/gnu/mcran4.cpp @@ -39,17 +39,19 @@ contained the header: Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ -#include <../../nrnconf.h> #include #include #include #include #include -#include -#include "hocdec.h" +#include "mcran4.h" static uint32_t lowindex = 0; +double mcell_lowindex() { + return static_cast(lowindex); +} + void mcell_ran4_init(uint32_t low) { lowindex = low; } @@ -70,40 +72,6 @@ uint32_t mcell_iran4(uint32_t* high) { return nrnRan4int(high, lowindex); } -/* Hoc interface */ -extern double chkarg(); -extern int use_mcell_ran4_; - - -void hoc_mcran4() { - uint32_t idx; - double* xidx; - double x; - xidx = hoc_pgetarg(1); - idx = (uint32_t) (*xidx); - x = mcell_ran4a(&idx); - *xidx = idx; - hoc_ret(); - hoc_pushx(x); -} -void hoc_mcran4init() { - double prev = (double) lowindex; - if (ifarg(1)) { - uint32_t idx = (uint32_t) chkarg(1, 0., 4294967295.); - mcell_ran4_init(idx); - } - hoc_ret(); - hoc_pushx(prev); -} -void hoc_usemcran4() { - double prev = (double) use_mcell_ran4_; - if (ifarg(1)) { - use_mcell_ran4_ = (int) chkarg(1, 0., 1.); - } - hoc_ret(); - hoc_pushx(prev); -} - uint32_t nrnRan4int(uint32_t* idx1, uint32_t idx2) { uint32_t u, v, w, m, n; /* 64-bit hash */ diff --git a/src/oc/mcran4.h b/src/gnu/mcran4.h similarity index 98% rename from src/oc/mcran4.h rename to src/gnu/mcran4.h index 7d251b4cc1..9d7ff79dd9 100644 --- a/src/oc/mcran4.h +++ b/src/gnu/mcran4.h @@ -1,6 +1,7 @@ #pragma once #include +double mcell_lowindex(); void mcell_ran4_init(uint32_t); double mcell_ran4(uint32_t* idx1, double* x, unsigned int n, double range); double mcell_ran4a(uint32_t* idx1); diff --git a/src/gnu/nrnisaac.cpp b/src/gnu/nrnisaac.cpp new file mode 100644 index 0000000000..f0549a4d9a --- /dev/null +++ b/src/gnu/nrnisaac.cpp @@ -0,0 +1,25 @@ +#include +#include "nrnisaac.h" +#include "isaac64.h" + +using RNG = struct isaac64_state; + +void* nrnisaac_new() { + return new RNG; +} + +void nrnisaac_delete(void* v) { + delete static_cast(v); +} + +void nrnisaac_init(void* v, unsigned long int seed) { + isaac64_init(static_cast(v), seed); +} + +double nrnisaac_dbl_pick(void* v) { + return isaac64_dbl32(static_cast(v)); +} + +std::uint32_t nrnisaac_uint32_pick(void* v) { + return isaac64_uint32(static_cast(v)); +} diff --git a/src/gnu/nrnisaac.h b/src/gnu/nrnisaac.h new file mode 100644 index 0000000000..5322febafe --- /dev/null +++ b/src/gnu/nrnisaac.h @@ -0,0 +1,9 @@ +#pragma once + +#include + +void* nrnisaac_new(); +void nrnisaac_delete(void* rng); +void nrnisaac_init(void* rng, unsigned long int seed); +double nrnisaac_dbl_pick(void* rng); +std::uint32_t nrnisaac_uint32_pick(void* rng); diff --git a/src/gnu/nrnran123.cpp b/src/gnu/nrnran123.cpp new file mode 100644 index 0000000000..eeb2768115 --- /dev/null +++ b/src/gnu/nrnran123.cpp @@ -0,0 +1,131 @@ +#include +#include "nrnran123.h" +#include +#include +#include + +using RNG = r123::Philox4x32; + +static RNG::key_type k = {{0}}; + +struct nrnran123_State { + RNG::ctr_type c; + RNG::ctr_type r; + char which_; +}; + +void nrnran123_set_globalindex(std::uint32_t gix) { + k[0] = gix; +} + +/* if one sets the global, one should reset all the stream sequences. */ +std::uint32_t nrnran123_get_globalindex() { + return k[0]; +} + +nrnran123_State* nrnran123_newstream(std::uint32_t id1, std::uint32_t id2) { + return nrnran123_newstream3(id1, id2, 0); +} +nrnran123_State* nrnran123_newstream3(std::uint32_t id1, std::uint32_t id2, std::uint32_t id3) { + auto* s = new nrnran123_State; + s->c[1] = id3; + s->c[2] = id1; + s->c[3] = id2; + nrnran123_setseq(s, 0, 0); + return s; +} + +void nrnran123_deletestream(nrnran123_State* s) { + delete s; +} + +void nrnran123_getseq(nrnran123_State* s, std::uint32_t* seq, char* which) { + *seq = s->c[0]; + *which = s->which_; +} + +void nrnran123_setseq(nrnran123_State* s, std::uint32_t seq, char which) { + if (which > 3 || which < 0) { + s->which_ = 0; + } else { + s->which_ = which; + } + s->c[0] = seq; + s->r = philox4x32(s->c, k); +} + +void nrnran123_getids(nrnran123_State* s, std::uint32_t* id1, std::uint32_t* id2) { + *id1 = s->c[2]; + *id2 = s->c[3]; +} + +void nrnran123_getids3(nrnran123_State* s, + std::uint32_t* id1, + std::uint32_t* id2, + std::uint32_t* id3) { + *id3 = s->c[1]; + *id1 = s->c[2]; + *id2 = s->c[3]; +} + +std::uint32_t nrnran123_ipick(nrnran123_State* s) { + char which = s->which_; + std::uint32_t rval = s->r[which++]; + if (which > 3) { + which = 0; + s->c.incr(); + s->r = philox4x32(s->c, k); + } + s->which_ = which; + return rval; +} + +double nrnran123_dblpick(nrnran123_State* s) { + static const double SHIFT32 = 1.0 / 4294967297.0; /* 1/(2^32 + 1) */ + auto u = nrnran123_ipick(s); + return ((double) u + 1.0) * SHIFT32; +} + +double nrnran123_negexp(nrnran123_State* s) { + /* min 2.3283064e-10 to max 22.18071 */ + return -std::log(nrnran123_dblpick(s)); +} + +/* At cost of a cached value we could compute two at a time. */ +/* But that would make it difficult to transfer to coreneuron for t > 0 */ +double nrnran123_normal(nrnran123_State* s) { + double w, x, y; + double u1, u2; + do { + u1 = nrnran123_dblpick(s); + u2 = nrnran123_dblpick(s); + u1 = 2. * u1 - 1.; + u2 = 2. * u2 - 1.; + w = (u1 * u1) + (u2 * u2); + } while (w > 1); + + y = std::sqrt((-2. * log(w)) / w); + x = u1 * y; + return x; +} + +nrnran123_array4x32 nrnran123_iran(std::uint32_t seq, std::uint32_t id1, std::uint32_t id2) { + return nrnran123_iran3(seq, id1, id2, 0); +} +nrnran123_array4x32 nrnran123_iran3(std::uint32_t seq, + std::uint32_t id1, + std::uint32_t id2, + std::uint32_t id3) { + nrnran123_array4x32 a; + RNG::ctr_type c; + c[0] = seq; + c[1] = id3; + c[2] = id1; + c[3] = id2; + RNG::ctr_type r = philox4x32(c, k); + a.v[0] = r[0]; + a.v[1] = r[1]; + a.v[2] = r[2]; + a.v[3] = r[3]; + return a; +} diff --git a/src/gnu/nrnran123.h b/src/gnu/nrnran123.h new file mode 100644 index 0000000000..a4e4844bf6 --- /dev/null +++ b/src/gnu/nrnran123.h @@ -0,0 +1,71 @@ +#ifndef nrnran123_h +#define nrnran123_h + +/* interface to Random123 */ +/* http://www.thesalmons.org/john/random123/papers/random123sc11.pdf */ + +/* +The 4x32 generators utilize a uint32x4 counter and uint32x4 key to transform +into an almost cryptographic quality uint32x4 random result. +There are many possibilites for balancing the sharing of the internal +state instances while reserving a uint32 counter for the stream sequence +and reserving other portions of the counter vector for stream identifiers +and global index used by all streams. + +We currently provide a single instance by default in which the policy is +to use the 0th counter uint32 as the stream sequence, words 2, 3 and 4 as the +stream identifier, and word 0 of the key as the global index. Unused words +are constant uint32 0. + +It is also possible to use Random123 directly without reference to this +interface. See Random123-1.02/docs/html/index.html +of the full distribution available from +http://www.deshawresearch.com/resources_random123.html +*/ + +#include + + +struct nrnran123_State; + +struct nrnran123_array4x32 { + std::uint32_t v[4]; +}; + +/* global index. eg. run number */ +/* all generator instances share this global index */ +extern void nrnran123_set_globalindex(std::uint32_t gix); +extern std::uint32_t nrnran123_get_globalindex(); + +/* minimal data stream */ +extern nrnran123_State* nrnran123_newstream(std::uint32_t id1, std::uint32_t id2); +extern nrnran123_State* nrnran123_newstream3(std::uint32_t id1, + std::uint32_t id2, + std::uint32_t id3); +extern void nrnran123_deletestream(nrnran123_State*); +extern void nrnran123_getseq(nrnran123_State*, std::uint32_t* seq, char* which); +extern void nrnran123_setseq(nrnran123_State*, std::uint32_t seq, char which); +extern void nrnran123_getids(nrnran123_State*, std::uint32_t* id1, std::uint32_t* id2); +extern void nrnran123_getids3(nrnran123_State*, + std::uint32_t* id1, + std::uint32_t* id2, + std::uint32_t* id3); + +// Get a random uint32_t in [0, 2^32-1] +extern std::uint32_t nrnran123_ipick(nrnran123_State*); +// Get a random double on [0, 1] +// nrnran123_dblpick minimum value is 2.3283064e-10 and max value is 1-min +extern double nrnran123_dblpick(nrnran123_State*); + +/* nrnran123_negexp min value is 2.3283064e-10, max is 22.18071 */ +extern double nrnran123_negexp(nrnran123_State*); /* mean 1.0 */ +extern double nrnran123_normal(nrnran123_State*); /* mean 0.0, std 1.0 */ + +/* more fundamental (stateless) (though the global index is still used) */ +extern nrnran123_array4x32 nrnran123_iran(std::uint32_t seq, std::uint32_t id1, std::uint32_t id2); +extern nrnran123_array4x32 nrnran123_iran3(std::uint32_t seq, + std::uint32_t id1, + std::uint32_t id2, + std::uint32_t id3); + +#endif diff --git a/src/ivoc/ivocrand.cpp b/src/ivoc/ivocrand.cpp index d3501913db..f250049299 100644 --- a/src/ivoc/ivocrand.cpp +++ b/src/ivoc/ivocrand.cpp @@ -4,7 +4,7 @@ #include #include -#include "random1.h" +#include "Rand.hpp" #include #include "classreg.h" @@ -17,7 +17,6 @@ #include "ocobserv.h" #include -#include #include #include #include @@ -33,6 +32,9 @@ #include #include #include +#include +#include +#include #if HAVE_IV #include "ivoc.h" @@ -59,118 +61,6 @@ class RandomPlay: public Observer, public Resource { using RandomPlayList = std::vector; static RandomPlayList* random_play_list_; -#include - -class NrnRandom123: public RNG { - public: - NrnRandom123(uint32_t id1, uint32_t id2, uint32_t id3 = 0); - virtual ~NrnRandom123(); - virtual uint32_t asLong() { - return nrnran123_ipick(s_); - } - virtual double asDouble() { - return nrnran123_dblpick(s_); - } - virtual void reset() { - nrnran123_setseq(s_, 0, 0); - } - nrnran123_State* s_; -}; -NrnRandom123::NrnRandom123(uint32_t id1, uint32_t id2, uint32_t id3) { - s_ = nrnran123_newstream3(id1, id2, id3); -} -NrnRandom123::~NrnRandom123() { - nrnran123_deletestream(s_); -} - - -// The decision that has to be made is whether each generator instance -// should have its own seed or only one seed for all. We choose separate -// seed for each but if arg not present or 0 then seed chosen by system. - -// the addition of ilow > 0 means that value is used for the lowindex -// instead of the mcell_ran4_init global 32 bit lowindex. - -class MCellRan4: public RNG { - public: - MCellRan4(uint32_t ihigh = 0, uint32_t ilow = 0); - virtual ~MCellRan4(); - virtual uint32_t asLong() { - return (uint32_t) (ilow_ == 0 ? mcell_iran4(&ihigh_) : nrnRan4int(&ihigh_, ilow_)); - } - virtual void reset() { - ihigh_ = orig_; - } - virtual double asDouble() { - return (ilow_ == 0 ? mcell_ran4a(&ihigh_) : nrnRan4dbl(&ihigh_, ilow_)); - } - uint32_t ihigh_; - uint32_t orig_; - uint32_t ilow_; - - private: - static uint32_t cnt_; -}; - -MCellRan4::MCellRan4(uint32_t ihigh, uint32_t ilow) { - ++cnt_; - ilow_ = ilow; - ihigh_ = ihigh; - if (ihigh_ == 0) { - ihigh_ = cnt_; - ihigh_ = (uint32_t) asLong(); - } - orig_ = ihigh_; -} -MCellRan4::~MCellRan4() {} - -uint32_t MCellRan4::cnt_ = 0; - -class Isaac64: public RNG { - public: - Isaac64(uint32_t seed = 0); - virtual ~Isaac64(); - virtual uint32_t asLong() { - return (uint32_t) nrnisaac_uint32_pick(rng_); - } - virtual void reset() { - nrnisaac_init(rng_, seed_); - } - virtual double asDouble() { - return nrnisaac_dbl_pick(rng_); - } - uint32_t seed() { - return seed_; - } - void seed(uint32_t s) { - seed_ = s; - reset(); - } - - private: - uint32_t seed_; - void* rng_; - static uint32_t cnt_; -}; - -Isaac64::Isaac64(uint32_t seed) { - if (cnt_ == 0) { - cnt_ = 0xffffffff; - } - --cnt_; - seed_ = seed; - if (seed_ == 0) { - seed_ = cnt_; - } - rng_ = nrnisaac_new(); - reset(); -} -Isaac64::~Isaac64() { - nrnisaac_delete(rng_); -} - -uint32_t Isaac64::cnt_ = 0; - RandomPlay::RandomPlay(Rand* r, neuron::container::data_handle px) : r_{r} , px_{std::move(px)} { @@ -197,27 +87,6 @@ void RandomPlay::update(Observable*) { list_remove(); } -Rand::Rand(unsigned long seed, int size, Object* obj) { - // printf("Rand\n"); - gen = new ACG(seed, size); - rand = new Normal(0., 1., gen); - type_ = 0; - obj_ = obj; -} - -Rand::~Rand() { - // printf("~Rand\n"); - delete gen; - delete rand; -} - -// constructor for a random number generator based on the RNG class -// from the gnu c++ class library -// defaults to the ACG generator (see below) - -// syntax: -// a = new Rand([seed],[size]) - static void* r_cons(Object* obj) { unsigned long seed = 0; int size = 55; @@ -337,8 +206,12 @@ static double r_nrnran123(void* r) { id2 = (uint32_t) (chkarg(2, 0., dmaxuint)); if (ifarg(3)) id3 = (uint32_t) (chkarg(3, 0., dmaxuint)); - NrnRandom123* r123 = new NrnRandom123(id1, id2, id3); - x->rand->generator(r123); + try { + NrnRandom123* r123 = new NrnRandom123(id1, id2, id3); + x->rand->generator(r123); + } catch (const std::bad_alloc& e) { + hoc_execerror("Bad allocation for 'NrnRandom123'", e.what()); + } delete x->gen; x->gen = x->rand->generator(); x->type_ = 4; @@ -402,14 +275,22 @@ static double r_Isaac64(void* r) { uint32_t seed1 = 0; - if (ifarg(1)) - seed1 = (uint32_t) (*getarg(1)); - Isaac64* mcr = new Isaac64(seed1); - x->rand->generator(mcr); - delete x->gen; - x->gen = x->rand->generator(); - x->type_ = 3; - return (double) mcr->seed(); + if (ifarg(1)) { + seed1 = static_cast(*getarg(1)); + } + + double seed{}; + try { + Isaac64* mcr = new Isaac64(seed1); + x->rand->generator(mcr); + delete x->gen; + x->gen = x->rand->generator(); + x->type_ = 3; + seed = mcr->seed(); + } catch (const std::bad_alloc& e) { + hoc_execerror("Bad allocation for Isaac64 generator", e.what()); + } + return seed; } // Pick again from the distribution last used diff --git a/src/ivoc/ivocvect.cpp b/src/ivoc/ivocvect.cpp index cc6565bea3..7b9087d328 100644 --- a/src/ivoc/ivocvect.cpp +++ b/src/ivoc/ivocvect.cpp @@ -111,7 +111,7 @@ static double dmaxint_ = 9007199254740992; #include "ivocvect.h" // definition of random numer generator -#include "random1.h" +#include "Rand.hpp" #include #if HAVE_IV diff --git a/src/ivoc/ocfile.cpp b/src/ivoc/ocfile.cpp index ee5b24da7b..61216c9694 100644 --- a/src/ivoc/ocfile.cpp +++ b/src/ivoc/ocfile.cpp @@ -325,14 +325,7 @@ void OcFile::binary_mode() { Use File.seek(0) after opening or use a binary style read/write as first\n\ access to file."); } -#if defined(__MWERKS__) - // printf("can't switch to binary mode. No setmode\n"); - mode_[1] = 'b'; - mode_[2] = '\0'; - file_ = freopen(filename_.c_str(), mode_, file()); -#else setmode(fileno(file()), O_BINARY); -#endif binary_ = true; } } @@ -345,20 +338,6 @@ bool OcFile::open(const char* name, const char* type) { strcpy(mode_, type); #endif file_ = fopen(expand_env_var(name), type); -#if defined(FILE_OPEN_RETRY) && FILE_OPEN_RETRY > 0 - int i; - for (i = 0; !file_ && i < FILE_OPEN_RETRY; ++i) { - // retry occasionally needed on BlueGene - file_ = fopen(expand_env_var(name), type); - } - if (i > 0) { - if (file_) { - printf("%d opened %s after %d retries\n", nrnmpi_myid_world, name, i); - } else { - printf("%d open %s failed after %d retries\n", nrnmpi_myid_world, name, i); - } - } -#endif return is_open(); } diff --git a/src/ivoc/pwman.cpp b/src/ivoc/pwman.cpp index c393eb47a7..898ece83f1 100644 --- a/src/ivoc/pwman.cpp +++ b/src/ivoc/pwman.cpp @@ -3340,13 +3340,6 @@ char* ivoc_get_temp_file() { if (!tdir) { tdir = "/tmp"; } -#if defined(WIN32) && defined(__MWERKS__) - char tname[L_tmpnam + 1]; - tmpnam(tname); - auto const length = strlen(tdir) + 1 + strlen(tname) + 1; - tmpfile = new char[length]; - std::snprintf(tmpfile, length, "%s/%s", tdir, tname); -#else auto const length = strlen(tdir) + 1 + 9 + 1; tmpfile = new char[length]; std::snprintf(tmpfile, length, "%s/nrnXXXXXX", tdir); @@ -3359,7 +3352,6 @@ char* ivoc_get_temp_file() { #else mktemp(tmpfile); #endif -#endif #if defined(WIN32) tmpfile = hoc_back2forward(tmpfile); #endif diff --git a/src/ivoc/rect.h b/src/ivoc/rect.h index 37715f3d0e..74865e1db7 100644 --- a/src/ivoc/rect.h +++ b/src/ivoc/rect.h @@ -36,11 +36,6 @@ class Appear: public Glyph { static const Brush* db_; }; -#if defined(__MWERKS__) -#undef Rect -#define Rect ivoc_Rect -#endif - class Rect: public Appear { public: Rect(Coord left, @@ -63,11 +58,6 @@ class Rect: public Appear { Coord l_, b_, w_, h_; }; -#if defined(__MWERKS__) -#undef Line -#define Line ivoc_Line -#endif - class Line: public Appear { public: Line(Coord dx, Coord dy, const Color* color = NULL, const Brush* brush = NULL); @@ -117,11 +107,6 @@ class Triangle: public Appear { bool filled_; }; -#if defined(__MWERKS__) -#undef Rectangle -#define Rectangle ivoc_Rectangle -#endif - class Rectangle: public Appear { public: Rectangle(float height, diff --git a/src/mswin/mwprefix.h b/src/mswin/mwprefix.h index b28199b663..877b31f9ec 100644 --- a/src/mswin/mwprefix.h +++ b/src/mswin/mwprefix.h @@ -1,9 +1,5 @@ /* auto include file for metrowerks codewarrior for all nrn */ -#if __MWERKS__ >= 7 -#define _MSL_DIRENT_H -#else #include -#endif #pragma once off #ifndef __WIN32__ #define __WIN32__ 1 diff --git a/src/nrniv/CMakeLists.txt b/src/nrniv/CMakeLists.txt index babfe56cc4..48f1c16531 100644 --- a/src/nrniv/CMakeLists.txt +++ b/src/nrniv/CMakeLists.txt @@ -211,7 +211,6 @@ set(NRN_INCLUDE_DIRS ${PROJECT_BINARY_DIR}/src/parallel ${PROJECT_BINARY_DIR}/src/sundials ${PROJECT_BINARY_DIR}/src/sundials/shared - ${PROJECT_SOURCE_DIR}/external/Random123/include ${PROJECT_SOURCE_DIR}/src ${PROJECT_SOURCE_DIR}/src/gnu ${PROJECT_SOURCE_DIR}/src/nrncvode diff --git a/src/nrniv/geometry3d.cpp b/src/nrniv/geometry3d.cpp deleted file mode 100644 index e0ffe7839a..0000000000 --- a/src/nrniv/geometry3d.cpp +++ /dev/null @@ -1,824 +0,0 @@ -/* - This file contains code adapted from p3d.py in - http://code.google.com/p/pythonisosurfaces/source/checkout - which was released under the new BSD license. - accessed 31 July 2012 -*/ - -#include -#include -#include -#include -#include -//#include - -int geometry3d_find_triangles(double value0, - double value1, - double value2, - double value3, - double value4, - double value5, - double value6, - double value7, - double x0, - double x1, - double y0, - double y1, - double z0, - double z1, - double* out, - int offset); - -double geometry3d_llgramarea(double* p0, double* p1, double* p2); -double geometry3d_sum_area_of_triangles(double* tri_vec, int len); - -const int edgeTable[] = { - 0x0, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, - 0xd03, 0xe09, 0xf00, 0x190, 0x99, 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, 0x99c, 0x895, - 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90, 0x230, 0x339, 0x33, 0x13a, 0x636, 0x73f, 0x435, - 0x53c, 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30, 0x3a0, 0x2a9, 0x1a3, 0xaa, - 0x7a6, 0x6af, 0x5a5, 0x4ac, 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0, 0x460, - 0x569, 0x663, 0x76a, 0x66, 0x16f, 0x265, 0x36c, 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, - 0xa69, 0xb60, 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff, 0x3f5, 0x2fc, 0xdfc, 0xcf5, 0xfff, - 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0, 0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55, 0x15c, - 0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950, 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, - 0x2cf, 0x1c5, 0xcc, 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0, 0x8c0, 0x9c9, - 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, 0xcc, 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, - 0x7c0, 0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c, 0x15c, 0x55, 0x35f, 0x256, - 0x55a, 0x453, 0x759, 0x650, 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, 0x2fc, - 0x3f5, 0xff, 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0, 0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, - 0xd65, 0xc6c, 0x36c, 0x265, 0x16f, 0x66, 0x76a, 0x663, 0x569, 0x460, 0xca0, 0xda9, 0xea3, - 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac, 0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa, 0x1a3, 0x2a9, 0x3a0, - 0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c, 0x53c, 0x435, 0x73f, 0x636, 0x13a, - 0x33, 0x339, 0x230, 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c, 0x69c, 0x795, - 0x49f, 0x596, 0x29a, 0x393, 0x99, 0x190, 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, - 0x80c, 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0}; - -/* CTNG:tritable */ -const int triTable[256][16] = {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1}, - {3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1}, - {3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1}, - {3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1}, - {9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1}, - {1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1}, - {9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1}, - {2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1}, - {8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1}, - {9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1}, - {4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1}, - {3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1}, - {1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1}, - {4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1}, - {4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1}, - {9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1}, - {1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1}, - {5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1}, - {2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1}, - {9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1}, - {0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1}, - {2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1}, - {10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1}, - {4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1}, - {5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1}, - {5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1}, - {9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1}, - {0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1}, - {1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1}, - {10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1}, - {8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1}, - {2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1}, - {7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1}, - {9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1}, - {2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1}, - {11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1}, - {9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1}, - {5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1}, - {11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1}, - {11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1}, - {1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1}, - {9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1}, - {5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1}, - {2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1}, - {0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1}, - {5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1}, - {6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1}, - {0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1}, - {3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1}, - {6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1}, - {5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1}, - {1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1}, - {10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1}, - {6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1}, - {1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1}, - {8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1}, - {7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1}, - {3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1}, - {5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1}, - {0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1}, - {9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1}, - {8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1}, - {5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1}, - {0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1}, - {6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1}, - {10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1}, - {10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1}, - {8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1}, - {1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1}, - {3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1}, - {0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1}, - {10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1}, - {0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1}, - {3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1}, - {6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1}, - {9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1}, - {8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1}, - {3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1}, - {6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1}, - {0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1}, - {10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1}, - {10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1}, - {1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1}, - {2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1}, - {7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1}, - {7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1}, - {2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1}, - {1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1}, - {11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1}, - {8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1}, - {0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1}, - {7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1}, - {10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1}, - {2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1}, - {6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1}, - {7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1}, - {2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1}, - {1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1}, - {10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1}, - {10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1}, - {0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1}, - {7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1}, - {6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1}, - {8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1}, - {9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1}, - {6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1}, - {1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1}, - {4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1}, - {10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1}, - {8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1}, - {0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1}, - {1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1}, - {8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1}, - {10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1}, - {4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1}, - {10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1}, - {5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1}, - {11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1}, - {9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1}, - {6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1}, - {7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1}, - {3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1}, - {7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1}, - {9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1}, - {3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1}, - {6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1}, - {9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1}, - {1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1}, - {4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1}, - {7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1}, - {6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1}, - {3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1}, - {0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1}, - {6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1}, - {1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1}, - {0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1}, - {11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1}, - {6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1}, - {5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1}, - {9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1}, - {1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1}, - {1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1}, - {10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1}, - {0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1}, - {5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1}, - {10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1}, - {11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1}, - {0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1}, - {9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1}, - {7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1}, - {2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1}, - {8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1}, - {9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1}, - {9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1}, - {1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1}, - {9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1}, - {9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1}, - {5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1}, - {0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1}, - {10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1}, - {2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1}, - {0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1}, - {0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1}, - {9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1}, - {5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1}, - {3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1}, - {5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1}, - {8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1}, - {0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1}, - {9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1}, - {0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1}, - {1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1}, - {3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1}, - {4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1}, - {9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1}, - {11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1}, - {11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1}, - {2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1}, - {9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1}, - {3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1}, - {1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1}, - {4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1}, - {4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1}, - {0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1}, - {3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1}, - {3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1}, - {0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1}, - {9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1}, - {1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}; - - -/* CTNG:interpedge */ -void geometry3d_vi(double* p1, double* p2, double v1, double v2, double* out) { - if (fabs(v1) < 0.000000000001) { - out[0] = p1[0]; - out[1] = p1[1]; - out[2] = p1[2]; - return; - } - if (fabs(v2) < 0.000000000001) { - out[0] = p2[0]; - out[1] = p2[1]; - out[2] = p2[2]; - return; - } - double delta_v = v1 - v2; - if (fabs(delta_v) < 0.0000000001) { - out[0] = p1[0]; - out[1] = p1[1]; - out[2] = p1[2]; - return; - } - double mu = v1 / delta_v; - if (std::isnan(mu)) { - printf("geometry3d_vi error. delta_v = %g, v1 = %g, v2 = %g\n", delta_v, v1, v2); - } - out[0] = p1[0] + mu * (p2[0] - p1[0]); - out[1] = p1[1] + mu * (p2[1] - p1[1]); - out[2] = p1[2] + mu * (p2[2] - p1[2]); -} - -int geometry3d_find_triangles(double value0, - double value1, - double value2, - double value3, - double value4, - double value5, - double value6, - double value7, - double x0, - double x1, - double y0, - double y1, - double z0, - double z1, - double* out, - int offset) { - out = out + offset; - double position[8][3] = {{x0, y0, z0}, - {x1, y0, z0}, - {x1, y1, z0}, - {x0, y1, z0}, - {x0, y0, z1}, - {x1, y0, z1}, - {x1, y1, z1}, - {x0, y1, z1}}; - - /* CTNG:domarch */ - int cubeIndex = 0; - if (value0 < 0) - cubeIndex |= 1; - if (value1 < 0) - cubeIndex |= 2; - if (value2 < 0) - cubeIndex |= 4; - if (value3 < 0) - cubeIndex |= 8; - if (value4 < 0) - cubeIndex |= 16; - if (value5 < 0) - cubeIndex |= 32; - if (value6 < 0) - cubeIndex |= 64; - if (value7 < 0) - cubeIndex |= 128; - - int et = edgeTable[cubeIndex]; - - if (et == 0) - return 0; - - double vertexList[12][3]; - if (et & 1) - geometry3d_vi(position[0], position[1], value0, value1, vertexList[0]); - if (et & 2) - geometry3d_vi(position[1], position[2], value1, value2, vertexList[1]); - if (et & 4) - geometry3d_vi(position[2], position[3], value2, value3, vertexList[2]); - if (et & 8) - geometry3d_vi(position[3], position[0], value3, value0, vertexList[3]); - if (et & 16) - geometry3d_vi(position[4], position[5], value4, value5, vertexList[4]); - if (et & 32) - geometry3d_vi(position[5], position[6], value5, value6, vertexList[5]); - if (et & 64) - geometry3d_vi(position[6], position[7], value6, value7, vertexList[6]); - if (et & 128) - geometry3d_vi(position[7], position[4], value7, value4, vertexList[7]); - if (et & 256) - geometry3d_vi(position[0], position[4], value0, value4, vertexList[8]); - if (et & 512) - geometry3d_vi(position[1], position[5], value1, value5, vertexList[9]); - if (et & 1024) - geometry3d_vi(position[2], position[6], value2, value6, vertexList[10]); - if (et & 2048) - geometry3d_vi(position[3], position[7], value3, value7, vertexList[11]); - - int const* const tt = triTable[cubeIndex]; - int i, j, k, count; - for (i = 0, count = 0; i < 16; i += 3, count++) { - if (tt[i] == -1) - break; - for (k = 0; k < 3; k++) { - for (j = 0; j < 3; j++) - out[j] = vertexList[tt[i + k]][j]; - out += 3; - } - } - return count; -} - - -double geometry3d_llgramarea(double* p0, double* p1, double* p2) { - /* setup the vectors */ - double a[] = {p0[0] - p1[0], p0[1] - p1[1], p0[2] - p1[2]}; - double b[] = {p0[0] - p2[0], p0[1] - p2[1], p0[2] - p2[2]}; - - /* take the cross-product */ - double cpx = a[1] * b[2] - a[2] * b[1]; - double cpy = a[2] * b[0] - a[0] * b[2]; - double cpz = a[0] * b[1] - a[1] * b[0]; - return std::sqrt(cpx * cpx + cpy * cpy + cpz * cpz); -} - - -double geometry3d_sum_area_of_triangles(double* tri_vec, int len) { - double sum = 0; - for (int i = 0; i < len; i += 9) { - sum += geometry3d_llgramarea(tri_vec + i, tri_vec + i + 3, tri_vec + i + 6); - } - return sum / 2.; -} - - -/***************************************************************************** - * this contains cone, and cylinder code translated and modified - * from Barbier and Galin 2004's example code - * see /u/ramcd/spatial/experiments/one_time_tests/2012-06-28/cone.cpp - * on 7 june 2012, the original code was available online at - * http://jgt.akpeters.com/papers/BarbierGalin04/Cone-Sphere.zip - ****************************************************************************/ - - -class geometry3d_Cylinder { - public: - geometry3d_Cylinder(double x0, double y0, double z0, double x1, double y1, double z1, double r); - //~geometry3d_Cylinder(); - double signed_distance(double px, double py, double pz); - - private: - // double x0, y0, z0, x1, y1, z1, r; - double r, rr, axisx, axisy, axisz, cx, cy, cz, h; -}; - -geometry3d_Cylinder::geometry3d_Cylinder(double x0, - double y0, - double z0, - double x1, - double y1, - double z1, - double r) - : r(r) - , cx((x0 + x1) / 2.) - , cy((y0 + y1) / 2.) - , cz((z0 + z1) / 2.) - , rr(r * r) { - axisx = x1 - x0; - axisy = y1 - y0; - axisz = z1 - z0; - double axislength = std::sqrt(axisx * axisx + axisy * axisy + axisz * axisz); - axisx /= axislength; - axisy /= axislength; - axisz /= axislength; - h = axislength / 2.; -} - -double geometry3d_Cylinder::signed_distance(double px, double py, double pz) { - double const nx{px - cx}; - double const ny{py - cy}; - double const nz{pz - cz}; - double y{std::abs(axisx * nx + axisy * ny + axisz * nz)}; - double yy{y * y}; - double const xx{nx * nx + ny * ny + nz * nz - yy}; - if (y < h) { - return std::max(-std::abs(h - y), std::sqrt(xx) - r); - } else { - y -= h; - yy = y * y; - if (xx < rr) { - return std::abs(y); - } else { - double const x{std::sqrt(xx) - r}; - return std::sqrt(yy + x * x); - } - } -} - -void* geometry3d_new_Cylinder(double x0, - double y0, - double z0, - double x1, - double y1, - double z1, - double r) { - return new geometry3d_Cylinder(x0, y0, z0, x1, y1, z1, r); -} -void geometry3d_delete_Cylinder(void* ptr) { - delete (geometry3d_Cylinder*) ptr; -} -// TODO: add validation for ptr -double geometry3d_Cylinder_signed_distance(void* ptr, double px, double py, double pz) { - return ((geometry3d_Cylinder*) ptr)->signed_distance(px, py, pz); -} - - -class geometry3d_Cone { - public: - geometry3d_Cone(double x0, - double y0, - double z0, - double r0, - double x1, - double y1, - double z1, - double r1); - double signed_distance(double px, double py, double pz); - - private: - double axisx, axisy, axisz, h, rra, rrb, conelength; - double side1, side2, x0, y0, z0, r0, axislength; -}; - -geometry3d_Cone::geometry3d_Cone(double x0, - double y0, - double z0, - double r0, - double x1, - double y1, - double z1, - double r1) - : rra(r0 * r0) - , rrb(r1 * r1) - , x0(x0) - , y0(y0) - , z0(z0) - , r0(r0) { - // TODO: these are preconditions; the python assures them, but could/should - // take care of that here - assert(r1 <= r0); - assert(r1 >= 0); - axisx = x1 - x0; - axisy = y1 - y0; - axisz = z1 - z0; - axislength = std::sqrt(axisx * axisx + axisy * axisy + axisz * axisz); - axisx /= axislength; - axisy /= axislength; - axisz /= axislength; - h = axislength / 2.; - rra = r0 * r0; - rrb = r1 * r1; - conelength = std::sqrt((r1 - r0) * (r1 - r0) + axislength * axislength); - side1 = (r1 - r0) / conelength; - side2 = axislength / conelength; -} - -double geometry3d_Cone::signed_distance(double px, double py, double pz) { - double nx, ny, nz, y, yy, xx, x, rx, ry; - nx = px - x0; - ny = py - y0; - nz = pz - z0; - y = axisx * nx + axisy * ny + axisz * nz; - yy = y * y; - xx = nx * nx + ny * ny + nz * nz - yy; - // in principle, xx >= 0, but roundoff errors may cause trouble - if (xx < 0) - xx = 0; - if (y < 0) { - if (xx < rra) - return -y; - x = std::sqrt(xx) - r0; - return std::sqrt(x * x + yy); - } else if (xx < rrb) { - return y - axislength; - } else { - x = std::sqrt(xx) - r0; - if (y < 0) { - if (x < 0) - return y; - return std::sqrt(x * x + yy); - } else { - ry = x * side1 + y * side2; - if (ry < 0) - return std::sqrt(x * x + yy); - rx = x * side2 - y * side1; - if (ry > conelength) { - ry -= conelength; - return std::sqrt(rx * rx + ry * ry); - } else { - return rx; - } - } - } -} - - -void* geometry3d_new_Cone(double x0, - double y0, - double z0, - double r0, - double x1, - double y1, - double z1, - double r1) { - return new geometry3d_Cone(x0, y0, z0, r0, x1, y1, z1, r1); -} -void geometry3d_delete_Cone(void* ptr) { - delete (geometry3d_Cone*) ptr; -} -// TODO: add validation for ptr -double geometry3d_Cone_signed_distance(void* ptr, double px, double py, double pz) { - return ((geometry3d_Cone*) ptr)->signed_distance(px, py, pz); -} - - -class geometry3d_Sphere { - public: - geometry3d_Sphere(double x, double y, double z, double r); - double signed_distance(double px, double py, double pz); - - private: - double x, y, z, r; -}; - -geometry3d_Sphere::geometry3d_Sphere(double x, double y, double z, double r) - : x(x) - , y(y) - , z(z) - , r(r) {} - -double geometry3d_Sphere::signed_distance(double px, double py, double pz) { - return std::sqrt(std::pow(x - px, 2) + std::pow(y - py, 2) + std::pow(z - pz, 2)) - r; -} - -void* geometry3d_new_Sphere(double x, double y, double z, double r) { - return new geometry3d_Sphere(x, y, z, r); -} -void geometry3d_delete_Sphere(void* ptr) { - delete (geometry3d_Sphere*) ptr; -} -// TODO: add validation for ptr -double geometry3d_Sphere_signed_distance(void* ptr, double px, double py, double pz) { - return ((geometry3d_Sphere*) ptr)->signed_distance(px, py, pz); -} - -class geometry3d_Plane { - public: - geometry3d_Plane(double x, double y, double z, double nx, double ny, double nz); - double signed_distance(double px, double py, double pz); - - private: - double nx, ny, nz; - double d, mul; -}; - -geometry3d_Plane::geometry3d_Plane(double x, double y, double z, double nx, double ny, double nz) - : nx(nx) - , ny(ny) - , nz(nz) - , d(-(nx * x + ny * y + nz * z)) - , mul(1. / std::sqrt(nx * nx + ny * ny + nz * nz)) {} - -double geometry3d_Plane::signed_distance(double px, double py, double pz) { - return (nx * px + ny * py + nz * pz + d) * mul; -} - -void* geometry3d_new_Plane(double x, double y, double z, double nx, double ny, double nz) { - return new geometry3d_Plane(x, y, z, nx, ny, nz); -} -void geometry3d_delete_Plane(void* ptr) { - delete (geometry3d_Plane*) ptr; -} -// TODO: add validation for ptr -double geometry3d_Plane_signed_distance(void* ptr, double px, double py, double pz) { - return ((geometry3d_Plane*) ptr)->signed_distance(px, py, pz); -} - -/* - void print_numbers(PyObject *p) { - for (Py_ssize_t i = 0; i< PyList_Size(p); i++) { - PyObject* obj = PyList_GetItem(p, i); - printf("%g ", PyFloat_AsDouble(obj)); - } - printf("\n"); - } - - // TODO: it would be nice to remove the python dependence, because that - // limits us to mostly single threaded due to the global interpreter - // lock - int geometry3d_process_cell(int i, int j, int k, PyObject* objects, double* xs, double* ys, -double* zs, double* tridata, int start) { double x, y, z, x1, y1, z1, xx, yy, zz; x = xs[i]; y = -ys[j]; z = zs[k]; x1 = xs[i + 1]; y1 = ys[j + 1]; z1 = zs[k + 1]; double value[8], current_value; - PyObject* result; - PyObject* obj; - printf("inside process_cell\n"); - - // march around the cube - for (int m = 0; m < 8; m++) { - // NOTE: only describing changes from previous case - switch(m) { - case 0: xx = x; yy = y; zz = z; break; - case 1: xx = x1; break; - case 2: yy = y1; break; - case 3: xx = x; break; - case 4: yy = y; zz = z1; break; - case 5: xx = x1; break; - case 6: yy = y1; break; - case 7: xx = x; break; - } - printf("phase 0, len(objects) = %ld\n", PyList_Size(objects)); - obj = PyList_GetItem(objects, 0); - printf("phase 0a, obj = %x\n", (void*) obj); - result = PyEval_CallMethod(obj, "distance", "ddd", xx, yy, zz); - //Py_DECREF(obj); - printf("phase 1\n"); - current_value = PyFloat_AsDouble(result); - //Py_DECREF(result); - for (Py_ssize_t n = 1; n < PyList_Size(objects); n++) { - printf("phase 2, n = %ld\n", n); - obj = PyList_GetItem(objects, n); - result = PyEval_CallMethod(obj, "distance", "ddd", xx, yy, zz); - Py_DECREF(obj); - current_value = min(current_value, PyFloat_AsDouble(result)); - //Py_DECREF(result); - } - value[m] = current_value; - } - printf("finishing up; start = %d\n", start); - return start + 9 * geometry3d_find_triangles(value[0], value[1], value[2], value[3], -value[4], value[5], value[6], value[7], x, x1, y, y1, z, z1, tridata, start); - } - - - int geometry3d_test_call_function(PyObject* obj) { - printf("inside\n"); - if (obj == NULL) printf("obj is NULL\n"); - Py_INCREF(obj); - PyEval_CallObject(obj, obj); - return 0; - } - - int geometry3d_test_call_function3(int (*func) (void)) { - return func(); - } - - int geometry3d_test_call_method(PyObject* list, PyObject* obj) { - PyEval_CallMethod(list, "insert", "O", obj); - return 0; - } - - // this works, but is slower than the python version - int geometry3d_contains_surface(int i, int j, int k, double (*objdist)(double, double, double), -double* xs, double* ys, double* zs, double dx, double r_inner, double r_outer) { bool has_neg = -false, has_pos = false; double xbar = xs[i] + dx / 2; double ybar = ys[j] + dx / 2; double zbar = -zs[k] + dx / 2; - - double d = fabs(objdist(xbar, ybar, zbar)); - //if (i == 586 && j == 2169 && k == 83) {printf("at magic point, d=%g\n", d);} - //printf("sphere test: d = %g, r_inner = %g, r_outer = %g\n", d, r_inner, r_outer); - if (d <= r_inner) return 1; - if (d >= r_outer) return 0; -// // spheres alone are indeterminant; check corners -// for (int di = 0; di < 2; di++) { -// for (int dj = 0; dj < 2; dj++) { -// for (int dk = 0; dk < 2; dk++) { -// d = objdist(xs[i + di], ys[j + dj], zs[k + dk]); -// //printf("d = %g\n", d); -// if (d <= 0) has_neg = true; -// if (d >= 0) has_pos = true; -// if (has_neg && has_pos) return 1; -// } -// } -// } - - // spheres alone are indeterminant; check corners - for (double* x = xs + i; x < xs + i + 2; x++) { - for (double* y = ys + j; y < ys + j + 2; y++) { - for (double* z = zs + k; z < zs + k + 2; z++) { - d = objdist(*x, *y, *z); - if (d <= 0) has_neg = true; - if (d >= 0) has_pos = true; - if (has_neg && has_pos) return 1; - } - } - } - - return 0; - } - -} -*/ diff --git a/src/nrniv/have2want.cpp b/src/nrniv/have2want.cpp deleted file mode 100644 index 5ac55709aa..0000000000 --- a/src/nrniv/have2want.cpp +++ /dev/null @@ -1,259 +0,0 @@ -/* -To be included by a file that desires rendezvous rank exchange functionality. -Need to define HAVEWANT_t, HAVEWANT_alltoallv, and HAVEWANT2Int -The latter is a map or unordered_map. -E.g. std::unordered_map -*/ - -#ifdef have2want_cpp -#error "This implementation can only be included once" -// The static function names used to involve a macro name (NrnHash) but now, -// with the use of std::..., it may be the case this could be included -// multiple times or even transformed into a template. -#endif - -#define have2want_cpp - -/* - -A rank owns a set of HAVEWANT_t keys and wants information associated with -a set of HAVEWANT_t keys owned by unknown ranks. Owners do not know which -ranks want their information. Ranks that want info do not know which ranks -own that info. - -The have_to_want function returns two new vectors of keys along with -associated count and displacement vectors of length nhost and nhost+1 -respectively. Note that a send_to_want_displ[i+1] = - send_to_want_cnt[i] + send_to_want_displ[i] . - -send_to_want[send_to_want_displ[i] to send_to_want_displ[i+1]] contains -the keys from this rank for which rank i wants information. - -recv_from_have[recv_from_have_displ[i] to recv_from_have_displ[i+1] contains -the keys from which rank i is sending information to this rank. - -Note that on rank i, the order of keys in the rank j area of send_to_want -is the same order of keys on rank j in the ith area in recv_from_have. - -The rendezvous_rank function is used to parallelize this computation -and minimize memory usage so that no single rank ever needs to know all keys. -*/ - -#ifndef HAVEWANT_t -#define HAVEWANT_t int -#endif - -// round robin default rendezvous rank function -static int default_rendezvous(HAVEWANT_t key) { - return key % nrnmpi_numprocs; -} - -static int* cnt2displ(int* cnt) { - int* displ = new int[nrnmpi_numprocs + 1]; - displ[0] = 0; - for (int i = 0; i < nrnmpi_numprocs; ++i) { - displ[i + 1] = displ[i] + cnt[i]; - } - return displ; -} - -static int* srccnt2destcnt(int* srccnt) { - int* destcnt = new int[nrnmpi_numprocs]; - nrnmpi_int_alltoall(srccnt, destcnt, 1); - return destcnt; -} - -static void rendezvous_rank_get(HAVEWANT_t* data, - int size, - HAVEWANT_t*& sdata, - int*& scnt, - int*& sdispl, - HAVEWANT_t*& rdata, - int*& rcnt, - int*& rdispl, - int (*rendezvous_rank)(HAVEWANT_t)) { - int nhost = nrnmpi_numprocs; - - // count what gets sent - scnt = new int[nhost]; - for (int i = 0; i < nhost; ++i) { - scnt[i] = 0; - } - for (int i = 0; i < size; ++i) { - int r = (*rendezvous_rank)(data[i]); - ++scnt[r]; - } - - sdispl = cnt2displ(scnt); - rcnt = srccnt2destcnt(scnt); - rdispl = cnt2displ(rcnt); - sdata = new HAVEWANT_t[sdispl[nhost] + 1]; // ensure not 0 size - rdata = new HAVEWANT_t[rdispl[nhost] + 1]; // ensure not 0 size - // scatter data into sdata by recalculating scnt. - for (int i = 0; i < nhost; ++i) { - scnt[i] = 0; - } - for (int i = 0; i < size; ++i) { - int r = (*rendezvous_rank)(data[i]); - sdata[sdispl[r] + scnt[r]] = data[i]; - ++scnt[r]; - } - HAVEWANT_alltoallv(sdata, scnt, sdispl, rdata, rcnt, rdispl); -} - -static void have_to_want(HAVEWANT_t* have, - int have_size, - HAVEWANT_t* want, - int want_size, - HAVEWANT_t*& send_to_want, - int*& send_to_want_cnt, - int*& send_to_want_displ, - HAVEWANT_t*& recv_from_have, - int*& recv_from_have_cnt, - int*& recv_from_have_displ, - int (*rendezvous_rank)(HAVEWANT_t)) { - // 1) Send have and want to the rendezvous ranks. - // 2) Rendezvous rank matches have and want. - // 3) Rendezvous ranks tell the want ranks which ranks own the keys - // 4) Ranks that want tell owner ranks where to send. - - int nhost = nrnmpi_numprocs; - - // 1) Send have and want to the rendezvous ranks. - HAVEWANT_t *have_s_data, *have_r_data; - int *have_s_cnt, *have_s_displ, *have_r_cnt, *have_r_displ; - rendezvous_rank_get(have, - have_size, - have_s_data, - have_s_cnt, - have_s_displ, - have_r_data, - have_r_cnt, - have_r_displ, - rendezvous_rank); - delete[] have_s_cnt; - delete[] have_s_displ; - delete[] have_s_data; - // assume it is an error if two ranks have the same key so create - // hash table of key2rank. Will also need it for matching have and want - HAVEWANT2Int havekey2rank = HAVEWANT2Int(have_r_displ[nhost] + 1); // ensure not empty. - for (int r = 0; r < nhost; ++r) { - for (int i = 0; i < have_r_cnt[r]; ++i) { - HAVEWANT_t key = have_r_data[have_r_displ[r] + i]; - if (havekey2rank.find(key) != havekey2rank.end()) { - hoc_execerr_ext( - "internal error in have_to_want: key %lld owned by multiple ranks\n", - (long long) key); - } - havekey2rank[key] = r; - } - } - delete[] have_r_data; - delete[] have_r_cnt; - delete[] have_r_displ; - - HAVEWANT_t *want_s_data, *want_r_data; - int *want_s_cnt, *want_s_displ, *want_r_cnt, *want_r_displ; - rendezvous_rank_get(want, - want_size, - want_s_data, - want_s_cnt, - want_s_displ, - want_r_data, - want_r_cnt, - want_r_displ, - rendezvous_rank); - - // 2) Rendezvous rank matches have and want. - // we already have made the havekey2rank map. - // Create an array parallel to want_r_data which contains the ranks that - // have that data. - int n = want_r_displ[nhost]; - int* want_r_ownerranks = new int[n]; - for (int r = 0; r < nhost; ++r) { - for (int i = 0; i < want_r_cnt[r]; ++i) { - int ix = want_r_displ[r] + i; - HAVEWANT_t key = want_r_data[ix]; - auto search = havekey2rank.find(key); - if (search == havekey2rank.end()) { - hoc_execerr_ext( - "internal error in have_to_want: key = %lld is wanted but does not exist\n", - (long long) key); - } - want_r_ownerranks[ix] = search->second; - } - } - delete[] want_r_data; - - // 3) Rendezvous ranks tell the want ranks which ranks own the keys - // The ranks that want keys need to know the ranks that own those keys. - // The want_s_ownerranks will be parallel to the want_s_data. - // That is, each item defines the rank from which information associated - // with that key is coming from - int* want_s_ownerranks = new int[want_s_displ[nhost]]; - if (nrn_sparse_partrans > 0) { - nrnmpi_int_alltoallv_sparse(want_r_ownerranks, - want_r_cnt, - want_r_displ, - want_s_ownerranks, - want_s_cnt, - want_s_displ); - } else { - nrnmpi_int_alltoallv(want_r_ownerranks, - want_r_cnt, - want_r_displ, - want_s_ownerranks, - want_s_cnt, - want_s_displ); - } - - delete[] want_r_ownerranks; - delete[] want_r_cnt; - delete[] want_r_displ; - - // 4) Ranks that want tell owner ranks where to send. - // Finished with the rendezvous ranks. The ranks that want keys know the - // owner ranks for those keys. The next step is for the want ranks to - // tell the owner ranks where to send. - // The parallel want_s_ownerranks and want_s_data are now uselessly ordered - // by rendezvous rank. Reorganize so that want ranks can tell owner ranks - // what they want. - n = want_s_displ[nhost]; - delete[] want_s_displ; - for (int i = 0; i < nhost; ++i) { - want_s_cnt[i] = 0; - } - HAVEWANT_t* old_want_s_data = want_s_data; - want_s_data = new HAVEWANT_t[n]; - // compute the counts - for (int i = 0; i < n; ++i) { - int r = want_s_ownerranks[i]; - ++want_s_cnt[r]; - } - want_s_displ = cnt2displ(want_s_cnt); - for (int i = 0; i < nhost; ++i) { - want_s_cnt[i] = 0; - } // recount while filling - for (int i = 0; i < n; ++i) { - int r = want_s_ownerranks[i]; - HAVEWANT_t key = old_want_s_data[i]; - want_s_data[want_s_displ[r] + want_s_cnt[r]] = key; - ++want_s_cnt[r]; - } - delete[] want_s_ownerranks; - delete[] old_want_s_data; - want_r_cnt = srccnt2destcnt(want_s_cnt); - want_r_displ = cnt2displ(want_r_cnt); - want_r_data = new HAVEWANT_t[want_r_displ[nhost]]; - HAVEWANT_alltoallv( - want_s_data, want_s_cnt, want_s_displ, want_r_data, want_r_cnt, want_r_displ); - // now the want_r_data on the have_ranks are grouped according to the ranks - // that want those keys. - - send_to_want = want_r_data; - send_to_want_cnt = want_r_cnt; - send_to_want_displ = want_r_displ; - recv_from_have = want_s_data; - recv_from_have_cnt = want_s_cnt; - recv_from_have_displ = want_s_displ; -} diff --git a/src/nrniv/have2want.hpp b/src/nrniv/have2want.hpp new file mode 100644 index 0000000000..11b22d28ed --- /dev/null +++ b/src/nrniv/have2want.hpp @@ -0,0 +1,202 @@ +#pragma once + +#include +#include +#include + +/* +A rank owns a set of T keys and wants information associated with +a set of T keys owned by unknown ranks. Owners do not know which +ranks want their information. Ranks that want info do not know which ranks +own that info. + +The have_to_want function returns two new vectors of keys along with +associated count and displacement vectors of length nhost and nhost+1 +respectively. Note that a send_to_want_displ[i+1] = + send_to_want_cnt[i] + send_to_want_displ[i] . + +send_to_want[send_to_want_displ[i] to send_to_want_displ[i+1]] contains +the keys from this rank for which rank i wants information. + +recv_from_have[recv_from_have_displ[i] to recv_from_have_displ[i+1] contains +the keys from which rank i is sending information to this rank. + +Note that on rank i, the order of keys in the rank j area of send_to_want +is the same order of keys on rank j in the ith area in recv_from_have. + +The rendezvous_rank function is used to parallelize this computation +and minimize memory usage so that no single rank ever needs to know all keys. +*/ + +// round robin rendezvous rank function +template +int rendezvous_rank(const T& key) { + return key % nrnmpi_numprocs; +} + +template +struct Data { + std::vector data{}; + std::vector cnt{}; + std::vector displ{}; +}; + +static std::vector cnt2displ(const std::vector& cnt) { + std::vector displ(nrnmpi_numprocs + 1); + std::partial_sum(cnt.cbegin(), cnt.cend(), displ.begin() + 1); + return displ; +} + +static std::vector srccnt2destcnt(std::vector srccnt) { + std::vector destcnt(nrnmpi_numprocs); + nrnmpi_int_alltoall(srccnt.data(), destcnt.data(), 1); + return destcnt; +} + +template +static std::tuple, Data> rendezvous_rank_get(const std::vector& data, + F alltoall_function) { + int nhost = nrnmpi_numprocs; + + Data s; + // count what gets sent + s.cnt = std::vector(nhost); + + for (const auto& e: data) { + int r = rendezvous_rank(e); + s.cnt[r] += 1; + } + + s.displ = cnt2displ(s.cnt); + s.data.resize(s.displ[nhost] + 1); + + Data r; + r.cnt = srccnt2destcnt(s.cnt); + r.displ = cnt2displ(r.cnt); + r.data.resize(r.displ[nhost]); + // scatter data into sdata by recalculating s.cnt. + std::fill(s.cnt.begin(), s.cnt.end(), 0); + for (const auto& e: data) { + int rank = rendezvous_rank(e); + s.data[s.displ[rank] + s.cnt[rank]] = e; + s.cnt[rank] += 1; + } + alltoall_function(s, r); + return {s, r}; +} + +template +std::pair, Data> have_to_want(const std::vector& have, + const std::vector& want, + F alltoall_function) { + // 1) Send have and want to the rendezvous ranks. + // 2) Rendezvous rank matches have and want. + // 3) Rendezvous ranks tell the want ranks which ranks own the keys + // 4) Ranks that want tell owner ranks where to send. + + int nhost = nrnmpi_numprocs; + + // 1) Send have and want to the rendezvous ranks. + + // hash table of key2rank. Will also need it for matching have and want + std::unordered_map havekey2rank{}; + { + auto [_, have_r] = rendezvous_rank_get(have, alltoall_function); + // assume it is an error if two ranks have the same key so create + havekey2rank.reserve(have_r.displ[nhost] + 1); + for (int r = 0; r < nhost; ++r) { + for (int i = 0; i < have_r.cnt[r]; ++i) { + T key = have_r.data[have_r.displ[r] + i]; + if (havekey2rank.find(key) != havekey2rank.end()) { + hoc_execerr_ext( + "internal error in have_to_want: key %lld owned by multiple ranks\n", + (long long) key); + } + havekey2rank[key] = r; + } + } + } + + auto [want_s, want_r] = rendezvous_rank_get(want, alltoall_function); + + // 2) Rendezvous rank matches have and want. + // we already have made the havekey2rank map. + // Create an array parallel to want_r_data which contains the ranks that + // have that data. + int n = want_r.displ[nhost]; + std::vector want_r_ownerranks(n); + for (int r = 0; r < nhost; ++r) { + for (int i = 0; i < want_r.cnt[r]; ++i) { + int ix = want_r.displ[r] + i; + T key = want_r.data[ix]; + auto search = havekey2rank.find(key); + if (search == havekey2rank.end()) { + hoc_execerr_ext( + "internal error in have_to_want: key = %lld is wanted but does not exist\n", + (long long) key); + } + want_r_ownerranks[ix] = search->second; + } + } + + // 3) Rendezvous ranks tell the want ranks which ranks own the keys + // The ranks that want keys need to know the ranks that own those keys. + // The want_s_ownerranks will be parallel to the want_s_data. + // That is, each item defines the rank from which information associated + // with that key is coming from + std::vector want_s_ownerranks(want_s.displ[nhost]); + if (nrn_sparse_partrans > 0) { + nrnmpi_int_alltoallv_sparse(want_r_ownerranks.data(), + want_r.cnt.data(), + want_r.displ.data(), + want_s_ownerranks.data(), + want_s.cnt.data(), + want_s.displ.data()); + } else { + nrnmpi_int_alltoallv(want_r_ownerranks.data(), + want_r.cnt.data(), + want_r.displ.data(), + want_s_ownerranks.data(), + want_s.cnt.data(), + want_s.displ.data()); + } + + // 4) Ranks that want tell owner ranks where to send. + // Finished with the rendezvous ranks. The ranks that want keys know the + // owner ranks for those keys. The next step is for the want ranks to + // tell the owner ranks where to send. + // The parallel want_s_ownerranks and want_s_data are now uselessly ordered + // by rendezvous rank. Reorganize so that want ranks can tell owner ranks + // what they want. + n = want_s.displ[nhost]; + for (int i = 0; i < nhost; ++i) { + want_s.cnt[i] = 0; + } + std::vector old_want_s_data(n); + std::swap(old_want_s_data, want_s.data); + // compute the counts + for (int i = 0; i < n; ++i) { + int r = want_s_ownerranks[i]; + ++want_s.cnt[r]; + } + want_s.displ = cnt2displ(want_s.cnt); + for (int i = 0; i < nhost; ++i) { + want_s.cnt[i] = 0; + } // recount while filling + for (int i = 0; i < n; ++i) { + int r = want_s_ownerranks[i]; + T key = old_want_s_data[i]; + want_s.data[want_s.displ[r] + want_s.cnt[r]] = key; + ++want_s.cnt[r]; + } + + Data new_want_r{}; + new_want_r.cnt = srccnt2destcnt(want_s.cnt); + new_want_r.displ = cnt2displ(new_want_r.cnt); + new_want_r.data.resize(new_want_r.displ[nhost]); + alltoall_function(want_s, new_want_r); + // now the want_r_data on the have_ranks are grouped according to the ranks + // that want those keys. + + return {new_want_r, want_s}; +} diff --git a/src/nrniv/hocmech.cpp b/src/nrniv/hocmech.cpp index 0235647100..e72027c07e 100644 --- a/src/nrniv/hocmech.cpp +++ b/src/nrniv/hocmech.cpp @@ -134,7 +134,7 @@ static void alloc_pnt(Prop* p) { p->ob = nrn_point_prop_->ob; // printf("p->ob comes from nrn_point_prop_ %s\n", hoc_object_name(p->ob)); } else { - p->dparam = (Datum*) hoc_Ecalloc(2, sizeof(Datum)); + nrn_prop_datum_alloc(p->_type, 2, p); if (last_created_pp_ob_) { p->ob = last_created_pp_ob_; // printf("p->ob comes from last_created %s\n", hoc_object_name(p->ob)); diff --git a/src/nrniv/kschan.cpp b/src/nrniv/kschan.cpp index f5e2fe0ac4..b3ab59ddc1 100644 --- a/src/nrniv/kschan.cpp +++ b/src/nrniv/kschan.cpp @@ -11,10 +11,6 @@ #include "nrniv_mf.h" #define NSingleIndex 0 -#if defined(__MWERKS__) && !defined(_MSC_VER) -#include -#define strdup _strdup -#endif using KSChanList = std::vector; static KSChanList* channels; diff --git a/src/nrniv/nonlinz.cpp b/src/nrniv/nonlinz.cpp index 65bb00dd5f..f0a24eac8d 100644 --- a/src/nrniv/nonlinz.cpp +++ b/src/nrniv/nonlinz.cpp @@ -6,56 +6,58 @@ #include "nrniv_mf.h" #include "nrnoc2iv.h" #include "nrnmpi.h" -#include "cspmatrix.h" #include "membfunc.h" +#include +#include + +using namespace std::complex_literals; + extern void v_setup_vectors(); extern int nrndae_extra_eqn_count(); extern Symlist* hoc_built_in_symlist; extern void (*nrnthread_v_transfer_)(NrnThread*); -extern spREAL* spGetElement(char*, int, int); -extern void pargap_jacobi_rhs(double*, double*); +extern void pargap_jacobi_rhs(std::vector>&, + const std::vector>&); extern void pargap_jacobi_setup(int mode); class NonLinImpRep { public: NonLinImpRep(); - virtual ~NonLinImpRep(); void delta(double); + + // Functions to fill the matrix void didv(); void dids(); void dsdv(); void dsds(); + int gapsolve(); - char* m_; + // Matrix containing the non linear system to solve. + Eigen::SparseMatrix> m_{}; + // The solver of the matrix using the LU decomposition method. + Eigen::SparseLU>> lu_{}; int scnt_; // structure_change int n_v_, n_ext_, n_lin_, n_ode_, neq_v_, neq_; std::vector> pv_, pvdot_; - int* v_index_; - double* rv_; - double* jv_; - double** diag_; - double* deltavec_; // just like cvode.atol*cvode.atolscale for ode's - double delta_; // slightly more efficient and easier for v. + std::vector> v_; + std::vector deltavec_; // just like cvode.atol*cvode.atolscale for ode's + double delta_; // slightly more efficient and easier for v. void current(int, Memb_list*, int); void ode(int, Memb_list*); double omega_; int iloc_; // current injection site of last solve - float* vsymtol_; - int maxiter_; + float* vsymtol_{}; + int maxiter_{500}; }; -NonLinImp::NonLinImp() { - rep_ = NULL; -} NonLinImp::~NonLinImp() { - if (rep_) { - delete rep_; - } + delete rep_; } + double NonLinImp::transfer_amp(int curloc, int vloc) { if (nrnmpi_numprocs > 1 && nrnthread_v_transfer_ && curloc != rep_->iloc_) { hoc_execerror( @@ -64,9 +66,7 @@ double NonLinImp::transfer_amp(int curloc, int vloc) { if (curloc != rep_->iloc_) { solve(curloc); } - double x = rep_->rv_[vloc]; - double y = rep_->jv_[vloc]; - return sqrt(x * x + y * y); + return std::abs(rep_->v_[vloc]); } double NonLinImp::input_amp(int curloc) { if (nrnmpi_numprocs > 1 && nrnthread_v_transfer_) { @@ -78,9 +78,7 @@ double NonLinImp::input_amp(int curloc) { if (curloc < 0) { return 0.0; } - double x = rep_->rv_[curloc]; - double y = rep_->jv_[curloc]; - return sqrt(x * x + y * y); + return std::abs(rep_->v_[curloc]); } double NonLinImp::transfer_phase(int curloc, int vloc) { if (nrnmpi_numprocs > 1 && nrnthread_v_transfer_ && curloc != rep_->iloc_) { @@ -90,9 +88,7 @@ double NonLinImp::transfer_phase(int curloc, int vloc) { if (curloc != rep_->iloc_) { solve(curloc); } - double x = rep_->rv_[vloc]; - double y = rep_->jv_[vloc]; - return atan2(y, x); + return std::arg(rep_->v_[vloc]); } double NonLinImp::input_phase(int curloc) { if (nrnmpi_numprocs > 1 && nrnthread_v_transfer_) { @@ -104,9 +100,7 @@ double NonLinImp::input_phase(int curloc) { if (curloc < 0) { return 0.0; } - double x = rep_->rv_[curloc]; - double y = rep_->jv_[curloc]; - return atan2(y, x); + return std::arg(rep_->v_[curloc]); } double NonLinImp::ratio_amp(int clmploc, int vloc) { if (nrnmpi_numprocs > 1 && nrnthread_v_transfer_) { @@ -118,22 +112,14 @@ double NonLinImp::ratio_amp(int clmploc, int vloc) { if (clmploc != rep_->iloc_) { solve(clmploc); } - double ax, bx, cx, ay, by, cy, bb; - ax = rep_->rv_[vloc]; - ay = rep_->jv_[vloc]; - bx = rep_->rv_[clmploc]; - by = rep_->jv_[clmploc]; - bb = bx * bx + by * by; - cx = (ax * bx + ay * by) / bb; - cy = (ay * bx - ax * by) / bb; - return sqrt(cx * cx + cy * cy); + return std::abs(rep_->v_[vloc] * std::conj(rep_->v_[clmploc]) / std::norm(rep_->v_[clmploc])); } void NonLinImp::compute(double omega, double deltafac, int maxiter) { v_setup_vectors(); nrn_rhs(nrn_ensure_model_data_are_sorted(), nrn_threads[0]); if (rep_ && rep_->scnt_ != structure_change_cnt) { delete rep_; - rep_ = NULL; + rep_ = nullptr; } if (!rep_) { rep_ = new NonLinImpRep(); @@ -151,8 +137,10 @@ void NonLinImp::compute(double omega, double deltafac, int maxiter) { rep_->omega_ = 1000. * omega; rep_->delta(deltafac); + + rep_->m_.setZero(); + // fill matrix - cmplx_spClear(rep_->m_); rep_->didv(); rep_->dsds(); #if 1 // when 0 equivalent to standard method @@ -160,18 +148,32 @@ void NonLinImp::compute(double omega, double deltafac, int maxiter) { rep_->dsdv(); #endif - // cmplx_spPrint(rep_->m_, 0, 1, 1); - // for (int i=0; i < rep_->neq_; ++i) { - // printf("i=%d %g %g\n", i, rep_->diag_[i][0], rep_->diag_[i][1]); - // } - int e = cmplx_spFactor(rep_->m_); - switch (e) { - case spZERO_DIAG: - hoc_execerror("cmplx_spFactor error:", "Zero Diagonal"); - case spNO_MEMORY: - hoc_execerror("cmplx_spFactor error:", "No Memory"); - case spSINGULAR: - hoc_execerror("cmplx_spFactor error:", "Singular"); + // Now that the matrix is filled we can compress it (mandatory for SparseLU) + rep_->m_.makeCompressed(); + + // Factorize the matrix so this is ready to solve + rep_->lu_.compute(rep_->m_); + switch (rep_->lu_.info()) { + case Eigen::Success: + // Everything fine + break; + case Eigen::NumericalIssue: + hoc_execerror( + "Eigen Sparse LU factorization failed with Eigen::NumericalIssue, please check the " + "input matrix:", + rep_->lu_.lastErrorMessage().c_str()); + break; + case Eigen::NoConvergence: + hoc_execerror( + "Eigen Sparse LU factorization reports Eigen::NonConvergence after calling compute():", + rep_->lu_.lastErrorMessage().c_str()); + break; + case Eigen::InvalidInput: + hoc_execerror( + "Eigen Sparse LU factorization failed with Eigen::InvalidInput, the input matrix seems " + "invalid:", + rep_->lu_.lastErrorMessage().c_str()); + break; } rep_->iloc_ = -2; @@ -184,20 +186,18 @@ int NonLinImp::solve(int curloc) { hoc_execerror("Must call Impedance.compute first", 0); } if (rep_->iloc_ != curloc) { - int i; rep_->iloc_ = curloc; - for (i = 0; i < rep_->neq_; ++i) { - rep_->rv_[i] = 0; - rep_->jv_[i] = 0; - } + rep_->v_ = std::vector>(rep_->neq_); if (curloc >= 0) { - rep_->rv_[curloc] = 1.e2 / NODEAREA(_nt->_v_node[curloc]); + rep_->v_[curloc] = 1.e2 / NODEAREA(_nt->_v_node[curloc]); } if (nrnthread_v_transfer_) { rval = rep_->gapsolve(); } else { - assert(rep_->m_); - cmplx_spSolve(rep_->m_, rep_->rv_ - 1, rep_->rv_ - 1, rep_->jv_ - 1, rep_->jv_ - 1); + auto v = + Eigen::Map, Eigen::Dynamic>>(rep_->v_.data(), + rep_->v_.size()); + v = rep_->lu_.solve(v); } } return rval; @@ -207,13 +207,8 @@ int NonLinImp::solve(int curloc) { // mapping is already done there. NonLinImpRep::NonLinImpRep() { - int err; - int i, j, ieq, cnt; NrnThread* _nt = nrn_threads; - maxiter_ = 500; - m_ = NULL; - vsymtol_ = NULL; Symbol* vsym = hoc_table_lookup("v", hoc_built_in_symlist); if (vsym->extra) { vsymtol_ = &vsym->extra->tolerance; @@ -233,10 +228,10 @@ NonLinImpRep::NonLinImpRep() { n_ode_ = 0; for (NrnThreadMembList* tml = _nt->tml; tml; tml = tml->next) { Memb_list* ml = tml->ml; - i = tml->index; + int i = tml->index; nrn_ode_count_t s = memb_func[i].ode_count; if (s) { - cnt = (*s)(i); + int cnt = (*s)(i); n_ode_ += cnt * ml->nodecount; } } @@ -245,46 +240,21 @@ NonLinImpRep::NonLinImpRep() { if (neq_ == 0) { return; } - m_ = cmplx_spCreate(neq_, 1, &err); - assert(err == spOKAY); + m_.resize(neq_, neq_); pv_.resize(neq_); pvdot_.resize(neq_); - v_index_ = new int[n_v_]; - rv_ = new double[neq_ + 1]; - rv_ += 1; - jv_ = new double[neq_ + 1]; - jv_ += 1; - diag_ = new double*[neq_]; - deltavec_ = new double[neq_]; - - for (i = 0; i < n_v_; ++i) { + v_.resize(neq_); + deltavec_.resize(neq_); + + for (int i = 0; i < n_v_; ++i) { // utilize nd->eqn_index in case of use_sparse13 later Node* nd = _nt->_v_node[i]; pv_[i] = nd->v_handle(); pvdot_[i] = nd->rhs_handle(); - v_index_[i] = i + 1; - } - for (i = 0; i < n_v_; ++i) { - diag_[i] = cmplx_spGetElement(m_, v_index_[i], v_index_[i]); - } - for (i = neq_v_; i < neq_; ++i) { - diag_[i] = cmplx_spGetElement(m_, i + 1, i + 1); } scnt_ = structure_change_cnt; } -NonLinImpRep::~NonLinImpRep() { - if (!m_) { - return; - } - cmplx_spDestroy(m_); - delete[] v_index_; - delete[](rv_ - 1); - delete[](jv_ - 1); - delete[] diag_; - delete[] deltavec_; -} - void NonLinImpRep::delta(double deltafac) { // also defines pv_,pvdot_ map for ode int i, nc, cnt, ieq; NrnThread* nt = nrn_threads; @@ -299,8 +269,12 @@ void NonLinImpRep::delta(double deltafac) { // also defines pv_,pvdot_ map for if (nrn_ode_count_t s = memb_func[i].ode_count; s && (cnt = s(i)) > 0) { nrn_ode_map_t ode_map = memb_func[i].ode_map; for (auto j = 0; j < nc; ++j) { - ode_map( - ml->prop[j], ieq, pv_.data() + ieq, pvdot_.data() + ieq, deltavec_ + ieq, i); + ode_map(ml->prop[j], + ieq, + pv_.data() + ieq, + pvdot_.data() + ieq, + deltavec_.data() + ieq, + i); ieq += cnt; } } @@ -317,19 +291,17 @@ void NonLinImpRep::didv() { for (i = _nt->ncell; i < n_v_; ++i) { nd = _nt->_v_node[i]; ip = _nt->_v_parent[i]->v_node_index; - double* a = cmplx_spGetElement(m_, v_index_[ip], v_index_[i]); - double* b = cmplx_spGetElement(m_, v_index_[i], v_index_[ip]); - *a += NODEA(nd); - *b += NODEB(nd); - *diag_[i] -= NODEB(nd); - *diag_[ip] -= NODEA(nd); + m_.coeffRef(ip, i) += NODEA(nd); + m_.coeffRef(i, ip) += NODEB(nd); + m_.coeffRef(i, i) -= NODEB(nd); + m_.coeffRef(ip, ip) -= NODEA(nd); } // jwC term Memb_list* mlc = _nt->tml->ml; int n = mlc->nodecount; for (i = 0; i < n; ++i) { j = mlc->nodelist[i]->v_node_index; - diag_[v_index_[j] - 1][1] += .001 * mlc->data(i, 0) * omega_; + m_.coeffRef(j, j) += .001 * mlc->data(i, 0) * omega_ * 1i; } // di/dv terms // because there may be several point processes of the same type @@ -370,7 +342,7 @@ void NonLinImpRep::didv() { current(i, ml, j); // conductance // add to matrix - *diag_[v_index_[nd->v_node_index] - 1] -= (x2 - NODERHS(nd)) / delta_; + m_.coeffRef(nd->v_node_index, nd->v_node_index) -= (x2 - NODERHS(nd)) / delta_; } } } @@ -390,8 +362,6 @@ void NonLinImpRep::dids() { nrn_ode_count_t s = memb_func[i].ode_count; int cnt = (*s)(i); if (memb_func[i].current) { - double* x1 = rv_; // use as temporary storage - double* x2 = jv_; for (in = 0; in < ml->nodecount; ++in) { Node* nd = ml->nodelist[in]; // zero rhs @@ -399,22 +369,20 @@ void NonLinImpRep::dids() { // compute rhs current(i, ml, in); // save rhs - x2[in] = NODERHS(nd); + v_[in].imag(NODERHS(nd)); // each state incremented separately and restored for (iis = 0; iis < cnt; ++iis) { is = ieq + in * cnt + iis; // save s - x1[is] = *pv_[is]; + v_[is].real(*pv_[is]); // increment s and zero rhs *pv_[is] += deltavec_[is]; NODERHS(nd) = 0; current(i, ml, in); - *pv_[is] = x1[is]; // restore s - double g = (NODERHS(nd) - x2[in]) / deltavec_[is]; + *pv_[is] = v_[is].real(); // restore s + double g = (NODERHS(nd) - v_[in].imag()) / deltavec_[is]; if (g != 0.) { - double* elm = - cmplx_spGetElement(m_, v_index_[nd->v_node_index], is + 1); - elm[0] = -g; + m_.coeffRef(nd->v_node_index, is) = -g; } } // don't know if this is necessary but make sure last @@ -439,22 +407,20 @@ void NonLinImpRep::dsdv() { nrn_ode_count_t s = memb_func[i].ode_count; int cnt = (*s)(i); if (memb_func[i].current) { - double* x1 = rv_; // use as temporary storage - double* x2 = jv_; // zero rhs, save v for (in = 0; in < ml->nodecount; ++in) { Node* nd = ml->nodelist[in]; for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) { *pvdot_[is] = 0.; } - x1[in] = NODEV(nd); + v_[in].real(NODEV(nd)); } // increment v only once in case there are multiple // point processes at the same location for (in = 0; in < ml->nodecount; ++in) { Node* nd = ml->nodelist[in]; auto const v = nd->v(); - if (x1[in] == v) { + if (v_[in].real() == v) { nd->v() = v + delta_; } } @@ -464,10 +430,10 @@ void NonLinImpRep::dsdv() { for (in = 0; in < ml->nodecount; ++in) { Node* nd = ml->nodelist[in]; for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) { - x2[is] = *pvdot_[is]; + v_[is].imag(*pvdot_[is]); *pvdot_[is] = 0; } - nd->v() = x1[in]; + nd->v() = v_[in].real(); } // compute the rhs(v) ode(i, ml); @@ -475,11 +441,9 @@ void NonLinImpRep::dsdv() { for (in = 0; in < ml->nodecount; ++in) { Node* nd = ml->nodelist[in]; for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) { - double ds = (x2[is] - *pvdot_[is]) / delta_; + double ds = (v_[is].imag() - *pvdot_[is]) / delta_; if (ds != 0.) { - double* elm = - cmplx_spGetElement(m_, is + 1, v_index_[nd->v_node_index]); - elm[0] = -ds; + m_.coeffRef(is, nd->v_node_index) = -ds; } } } @@ -494,7 +458,7 @@ void NonLinImpRep::dsds() { NrnThread* nt = nrn_threads; // jw term for (i = neq_v_; i < neq_; ++i) { - diag_[i][1] += omega_; + m_.coeffRef(i, i) += omega_ * 1i; } ieq = neq_v_; for (NrnThreadMembList* tml = nt->tml; tml; tml = tml->next) { @@ -504,13 +468,11 @@ void NonLinImpRep::dsds() { int nc = ml->nodecount; nrn_ode_count_t s = memb_func[i].ode_count; int cnt = (*s)(i); - double* x1 = rv_; // use as temporary storage - double* x2 = jv_; // zero rhs, save s for (in = 0; in < ml->nodecount; ++in) { for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) { *pvdot_[is] = 0.; - x1[is] = *pv_[is]; + v_[is].real(*pv_[is]); } } // compute rhs. this is the rhs(s) @@ -518,7 +480,7 @@ void NonLinImpRep::dsds() { // save rhs for (in = 0; in < ml->nodecount; ++in) { for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) { - x2[is] = *pvdot_[is]; + v_[is].imag(*pvdot_[is]); } } // iterate over the states @@ -538,12 +500,11 @@ void NonLinImpRep::dsds() { Node* nd = ml->nodelist[in]; ks = ieq + in * cnt + kks; for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) { - double ds = (*pvdot_[is] - x2[is]) / deltavec_[is]; + double ds = (*pvdot_[is] - v_[is].imag()) / deltavec_[is]; if (ds != 0.) { - double* elm = cmplx_spGetElement(m_, is + 1, ks + 1); - elm[0] = -ds; + m_.coeffRef(is, ks) = -ds; } - *pv_[ks] = x1[ks]; + *pv_[ks] = v_[ks].real(); } } // perhaps not necessary but ensures the last computation is with @@ -573,9 +534,11 @@ void NonLinImpRep::ode(int im, Memb_list* ml) { // assume there is in fact an o memb_func[im].ode_spec(nrn_ensure_model_data_are_sorted(), nrn_threads, ml, im); } +// This function compute a solution of a converging system by iteration. +// The value returned is the number of iterations to reach a precision of "tol" (1e-9). int NonLinImpRep::gapsolve() { - // On entry, rv_ and jv_ contain the complex b for A*x = b. - // On return rv_ and jv_ contain complex solution, x. + // On entry, v_ contains the complex b for A*x = b. + // On return v_ contains complex solution, x. // m_ is the factored matrix for the trees without gap junctions // Jacobi method (easy for parallel) // A = D + R @@ -592,48 +555,40 @@ int NonLinImpRep::gapsolve() { } #endif - pargap_jacobi_setup(0); + pargap_jacobi_setup(0); // 0 means 'setup' - double *rx, *jx, *rx1, *jx1, *rb, *jb; - if (neq_) { - rx = new double[neq_]; - jx = new double[neq_]; - rx1 = new double[neq_]; - jx1 = new double[neq_]; - rb = new double[neq_]; - jb = new double[neq_]; - } - - // initialize for first iteration - for (int i = 0; i < neq_; ++i) { - rx[i] = jx[i] = 0.0; - rb[i] = rv_[i]; - jb[i] = jv_[i]; - } + std::vector> x_old(neq_); + std::vector> x(neq_); + std::vector> b(v_); // iterate till change in x is small double tol = 1e-9; - double delta; + double delta{}; int success = 0; int iter; for (iter = 1; iter <= maxiter_; ++iter) { if (neq_) { - cmplx_spSolve(m_, rb - 1, rx1 - 1, jb - 1, jx1 - 1); + auto b_ = Eigen::Map, Eigen::Dynamic>>(b.data(), + b.size()); + auto x_ = Eigen::Map, Eigen::Dynamic>>(x.data(), + x.size()); + x_ = lu_.solve(b_); } // if any change in x > tol, then do another iteration. success = 1; delta = 0.0; + // Do the substraction of the previous result (x_old) and current result (x). + // If all differences are < tol stop the loop, otherwise continue to iterate for (int i = 0; i < neq_; ++i) { - double err = fabs(rx1[i] - rx[i]) + fabs(jx1[i] - jx[i]); + auto diff = x[i] - x_old[i]; + double err = std::abs(diff.real()) + std::abs(diff.imag()); if (err > tol) { success = 0; } - if (delta < err) { - delta = err; - } + delta = std::max(err, delta); } #if NRNMPI if (nrnmpi_numprocs > 1) { @@ -641,34 +596,17 @@ int NonLinImpRep::gapsolve() { } #endif if (success) { - for (int i = 0; i < neq_; ++i) { - rv_[i] = rx1[i]; - jv_[i] = jx1[i]; - } + v_ = x; break; } // setup for next iteration - for (int i = 0; i < neq_; ++i) { - rx[i] = rx1[i]; - jx[i] = jx1[i]; - rb[i] = rv_[i]; - jb[i] = jv_[i]; - } - pargap_jacobi_rhs(rb, rx); - pargap_jacobi_rhs(jb, jx); + x_old = x; + b = v_; + pargap_jacobi_rhs(b, x_old); } - pargap_jacobi_setup(1); // tear down - - if (neq_) { - delete[] rx; - delete[] jx; - delete[] rx1; - delete[] jx1; - delete[] rb; - delete[] jb; - } + pargap_jacobi_setup(1); // 1 means 'tear down' if (!success) { char buf[256]; @@ -678,7 +616,7 @@ int NonLinImpRep::gapsolve() { maxiter_, delta, tol); - execerror(buf, 0); + hoc_execerror(buf, nullptr); } return iter; } diff --git a/src/nrniv/nonlinz.h b/src/nrniv/nonlinz.h index 115a191b27..b477bf783a 100644 --- a/src/nrniv/nonlinz.h +++ b/src/nrniv/nonlinz.h @@ -3,20 +3,24 @@ class NonLinImpRep; +// A solver for non linear equation of complex numbers. +// Matrix should be squared. class NonLinImp { public: - NonLinImp(); - virtual ~NonLinImp(); + ~NonLinImp(); + // Prepare the matrix before solving it. void compute(double omega, double deltafac, int maxiter); - double transfer_amp(int curloc, int vloc); // v_node[arg] is the node + + double transfer_amp(int curloc, int vloc); double transfer_phase(int curloc, int vloc); double input_amp(int curloc); double input_phase(int curloc); double ratio_amp(int clmploc, int vloc); + int solve(int curloc); private: - NonLinImpRep* rep_; + NonLinImpRep* rep_{}; }; #endif diff --git a/src/nrniv/partrans.cpp b/src/nrniv/partrans.cpp index 1f34f9bbf2..a6779bad9e 100644 --- a/src/nrniv/partrans.cpp +++ b/src/nrniv/partrans.cpp @@ -13,27 +13,53 @@ #include #include +#include #include // Replaces NrnHash for MapSgid2Int #include #include +#if NRNMPI +#include "have2want.hpp" +#endif + + #if NRNLONGSGID #if NRNMPI -static void sgid_alltoallv(sgid_t* s, int* scnt, int* sdispl, sgid_t* r, int* rcnt, int* rdispl) { +static void sgid_alltoallv(Data& s, Data& r) { if (nrn_sparse_partrans > 0) { - nrnmpi_long_alltoallv_sparse(s, scnt, sdispl, r, rcnt, rdispl); + nrnmpi_long_alltoallv_sparse(s.data.data(), + s.cnt.data(), + s.displ.data(), + r.data.data(), + r.cnt.data(), + r.displ.data()); } else { - nrnmpi_long_alltoallv(s, scnt, sdispl, r, rcnt, rdispl); + nrnmpi_long_alltoallv(s.data.data(), + s.cnt.data(), + s.displ.data(), + r.data.data(), + r.cnt.data(), + r.displ.data()); } } #endif // NRNMPI #else // not NRNLONGSGID #if NRNMPI -static void sgid_alltoallv(sgid_t* s, int* scnt, int* sdispl, sgid_t* r, int* rcnt, int* rdispl) { +static void sgid_alltoallv(Data& s, Data& r) { if (nrn_sparse_partrans > 0) { - nrnmpi_int_alltoallv_sparse(s, scnt, sdispl, r, rcnt, rdispl); + nrnmpi_int_alltoallv_sparse(s.data.data(), + s.cnt.data(), + s.displ.data(), + r.data.data(), + r.cnt.data(), + r.displ.data()); } else { - nrnmpi_int_alltoallv(s, scnt, sdispl, r, rcnt, rdispl); + nrnmpi_int_alltoallv(s.data.data(), + s.cnt.data(), + s.displ.data(), + r.data.data(), + r.cnt.data(), + r.displ.data()); } } #endif // NRNMPI @@ -152,8 +178,12 @@ static std::vector> poutsrc_; // prior t // to proper place in // outsrc_buf_ static int* poutsrc_indices_; // for recalc pointers -static int insrc_buf_size_, *insrccnt_, *insrcdspl_; -static int outsrc_buf_size_, *outsrccnt_, *outsrcdspl_; +static int insrc_buf_size_; +static std::vector insrccnt_; +static std::vector insrcdspl_; +static int outsrc_buf_size_; +static std::vector outsrccnt_; +static std::vector outsrcdspl_; static MapSgid2Int sid2insrc_; // received interprocessor sid data is // associated with which insrc_buf index. Created by nrnmpi_setup_transfer // and used by mk_ttd @@ -522,11 +552,19 @@ static void mpi_transfer() { if (nrnmpi_numprocs > 1) { double wt = nrnmpi_wtime(); if (nrn_sparse_partrans > 0) { - nrnmpi_dbl_alltoallv_sparse( - outsrc_buf_, outsrccnt_, outsrcdspl_, insrc_buf_, insrccnt_, insrcdspl_); + nrnmpi_dbl_alltoallv_sparse(outsrc_buf_, + outsrccnt_.data(), + outsrcdspl_.data(), + insrc_buf_, + insrccnt_.data(), + insrcdspl_.data()); } else { - nrnmpi_dbl_alltoallv( - outsrc_buf_, outsrccnt_, outsrcdspl_, insrc_buf_, insrccnt_, insrcdspl_); + nrnmpi_dbl_alltoallv(outsrc_buf_, + outsrccnt_.data(), + outsrcdspl_.data(), + insrc_buf_, + insrccnt_.data(), + insrcdspl_.data()); } nrnmpi_transfer_wait_ += nrnmpi_wtime() - wt; errno = 0; @@ -577,15 +615,6 @@ static void thread_transfer(NrnThread* _nt) { // " But this was a mistake as many mpi implementations do not allow overlap // of send and receive buffers. -// 22-08-2014 For setup of the All2allv pattern, use the rendezvous rank -// idiom. -#define HAVEWANT_t sgid_t -#define HAVEWANT_alltoallv sgid_alltoallv -#define HAVEWANT2Int MapSgid2Int -#if NRNMPI -#include "have2want.cpp" -#endif - void nrnmpi_setup_transfer() { #if !NRNMPI if (nrnmpi_numprocs > 1) { @@ -611,10 +640,14 @@ void nrnmpi_setup_transfer() { return; } if (nrnmpi_numprocs > 1) { - delete[] std::exchange(insrccnt_, nullptr); - delete[] std::exchange(insrcdspl_, nullptr); - delete[] std::exchange(outsrccnt_, nullptr); - delete[] std::exchange(outsrcdspl_, nullptr); + insrccnt_.clear(); + insrccnt_.shrink_to_fit(); + insrcdspl_.clear(); + insrcdspl_.shrink_to_fit(); + outsrccnt_.clear(); + outsrccnt_.shrink_to_fit(); + outsrcdspl_.clear(); + outsrcdspl_.shrink_to_fit(); // This is an old comment prior to using the want_to_have rendezvous // rank function in want2have.cpp. The old method did not scale // to more sgids than could fit on a single rank, because @@ -649,78 +682,54 @@ void nrnmpi_setup_transfer() { // sids needed by this machine. The 'seen' table values are unused // but the keys are all the (unique) sgid needed by this process. // At the end seen is in fact what we want for sid2insrc_. - int needsrc_cnt = 0; int szalloc = targets_.size(); szalloc = szalloc ? szalloc : 1; // At the moment sid2insrc_ is serving as 'seen' sid2insrc_.clear(); - sid2insrc_.reserve(szalloc); // for single counting - sgid_t* needsrc = new sgid_t[szalloc]; // more than we need + sid2insrc_.reserve(szalloc); // for single counting + std::vector needsrc{}; for (size_t i = 0; i < sgid2targets_.size(); ++i) { sgid_t sid = sgid2targets_[i]; auto search = sid2insrc_.find(sid); if (search == sid2insrc_.end()) { sid2insrc_[sid] = 0; // at the moment, value does not matter - needsrc[needsrc_cnt++] = sid; + needsrc.push_back(sid); } } // 1 continued) Create an array of sources this rank owns. // This already exists as a vector in the SgidList sgids_ but // that is private so go ahead and copy. - sgid_t* ownsrc = new sgid_t[sgids_.size() + 1]; // not 0 length if count is 0 - for (size_t i = 0; i < sgids_.size(); ++i) { - ownsrc[i] = sgids_[i]; - } + std::vector ownsrc = sgids_; // 2) Call the have_to_want function. - sgid_t* send_to_want; - int *send_to_want_cnt, *send_to_want_displ; - sgid_t* recv_from_have; - int *recv_from_have_cnt, *recv_from_have_displ; - - have_to_want(ownsrc, - sgids_.size(), - needsrc, - needsrc_cnt, - send_to_want, - send_to_want_cnt, - send_to_want_displ, - recv_from_have, - recv_from_have_cnt, - recv_from_have_displ, - default_rendezvous); + auto [send_to_want, recv_from_have] = have_to_want(ownsrc, needsrc, sgid_alltoallv); // sanity check. all the sgids we are asked to send, we actually have - int n = send_to_want_displ[nhost]; + int n = send_to_want.displ[nhost]; // sanity check. all the sgids we receive, we actually need. // also set the sid2insrc_ value to the proper recv_from_have index. - n = recv_from_have_displ[nhost]; + n = recv_from_have.displ[nhost]; for (int i = 0; i < n; ++i) { - sgid_t sgid = recv_from_have[i]; + sgid_t sgid = recv_from_have.data[i]; nrn_assert(sid2insrc_.find(sgid) != sid2insrc_.end()); sid2insrc_[sgid] = i; } - // clean up a little - delete[] std::exchange(ownsrc, nullptr); - delete[] std::exchange(needsrc, nullptr); - delete[] std::exchange(recv_from_have, nullptr); - // 3) First return triple creates the proper outsrc_buf_. // Now that we know what machines are interested in our sids... // construct outsrc_buf, outsrc_buf_size, outsrccnt_, outsrcdspl_ // and poutsrc_; - outsrccnt_ = send_to_want_cnt; - outsrcdspl_ = send_to_want_displ; + std::swap(outsrccnt_, send_to_want.cnt); + std::swap(outsrcdspl_, send_to_want.displ); outsrc_buf_size_ = outsrcdspl_[nrnmpi_numprocs]; szalloc = std::max(1, outsrc_buf_size_); outsrc_buf_ = new double[szalloc]; poutsrc_.resize(szalloc); poutsrc_indices_ = new int[szalloc]; for (int i = 0; i < outsrc_buf_size_; ++i) { - sgid_t sid = send_to_want[i]; + sgid_t sid = send_to_want.data[i]; auto search = sgid2srcindex_.find(sid); nrn_assert(search != sgid2srcindex_.end()); Node* nd = visources_[search->second]; @@ -735,12 +744,11 @@ void nrnmpi_setup_transfer() { poutsrc_indices_[i] = search->second; outsrc_buf_[i] = double(sid); // see step 5 } - delete[] send_to_want; // 4) The second triple is creates the insrc_buf_. // From the recv_from_have and sid2insrc_ table, construct the insrc... - insrccnt_ = recv_from_have_cnt; - insrcdspl_ = recv_from_have_displ; + std::swap(insrccnt_, recv_from_have.cnt); + std::swap(insrcdspl_, recv_from_have.displ); insrc_buf_size_ = insrcdspl_[nrnmpi_numprocs]; szalloc = insrc_buf_size_ ? insrc_buf_size_ : 1; insrc_buf_ = new double[szalloc]; @@ -873,47 +881,59 @@ void pargap_jacobi_setup(int mode) { } } -void pargap_jacobi_rhs(double* b, double* x) { - // helper for complex impedance with parallel gap junctions - // b = b - R*x R are the off diagonal gap elements of the jacobian. - // we presume 1 thread. First nrn_thread[0].end equations are in node order. - if (!nrnthread_v_transfer_) { - return; - } +void pargap_jacobi_rhs(std::vector>& b, + const std::vector>& x) { + // First loop for real, second for imag + for (int real_imag = 0; real_imag < 2; ++real_imag) { + // helper for complex impedance with parallel gap junctions + // b = b - R*x R are the off diagonal gap elements of the jacobian. + // we presume 1 thread. First nrn_thread[0].end equations are in node order. + if (!nrnthread_v_transfer_) { + return; + } - NrnThread* _nt = nrn_threads; + NrnThread* _nt = nrn_threads; - // transfer gap node voltages to gap vpre - for (size_t i = 0; i < visources_.size(); ++i) { - Node* nd = visources_[i]; - nd->v() = x[nd->v_node_index]; - } - mpi_transfer(); - thread_transfer(_nt); + // transfer gap node voltages to gap vpre + for (size_t i = 0; i < visources_.size(); ++i) { + Node* nd = visources_[i]; + if (real_imag == 0) { + nd->v() = x[nd->v_node_index].real(); + } else { + nd->v() = x[nd->v_node_index].imag(); + } + } + mpi_transfer(); + thread_transfer(_nt); - // set gap node voltages to 0 so we can use nrn_cur to set rhs - for (size_t i = 0; i < visources_.size(); ++i) { - Node* nd = visources_[i]; - nd->v() = 0.0; - } - auto const sorted_token = nrn_ensure_model_data_are_sorted(); - auto* const vec_rhs = _nt->node_rhs_storage(); - // Initialize rhs to 0. - for (int i = 0; i < _nt->end; ++i) { - vec_rhs[i] = 0.0; - } - for (int k = 0; k < imped_current_type_count_; ++k) { - int type = imped_current_type_[k]; - Memb_list* ml = imped_current_ml_[k]; - memb_func[type].current(sorted_token, _nt, ml, type); - } + // set gap node voltages to 0 so we can use nrn_cur to set rhs + for (size_t i = 0; i < visources_.size(); ++i) { + Node* nd = visources_[i]; + nd->v() = 0.0; + } + auto const sorted_token = nrn_ensure_model_data_are_sorted(); + auto* const vec_rhs = _nt->node_rhs_storage(); + // Initialize rhs to 0. + for (int i = 0; i < _nt->end; ++i) { + vec_rhs[i] = 0.0; + } + for (int k = 0; k < imped_current_type_count_; ++k) { + int type = imped_current_type_[k]; + Memb_list* ml = imped_current_ml_[k]; + memb_func[type].current(sorted_token, _nt, ml, type); + } - // possibly many gap junctions in same node (and possible even different - // types) but rhs is the accumulation of all those instances at each node - // so ... The only thing that can go wrong is if there are intances of - // gap junctions that are not being used (not in the target list). - for (int i = 0; i < _nt->end; ++i) { - b[i] += vec_rhs[i]; + // possibly many gap junctions in same node (and possible even different + // types) but rhs is the accumulation of all those instances at each node + // so ... The only thing that can go wrong is if there are intances of + // gap junctions that are not being used (not in the target list). + for (int i = 0; i < _nt->end; ++i) { + if (real_imag == 0) { + b[i] += vec_rhs[i]; + } else { + b[i] += std::complex(0, vec_rhs[i]); + } + } } } diff --git a/src/oc/hoc_init.cpp b/src/oc/hoc_init.cpp index 6e022f7047..17675d7753 100644 --- a/src/oc/hoc_init.cpp +++ b/src/oc/hoc_init.cpp @@ -10,6 +10,7 @@ #include "nrn_ansi.h" #include "ocfunc.h" +#include "oc_mcran4.hpp" extern void hoc_nrnmpi_init(); @@ -231,7 +232,6 @@ double hoc_default_dll_loaded_; char* neuron_home; const char* nrn_mech_dll; /* but actually only for NEURON mswin and linux */ int nrn_noauto_dlopen_nrnmech; /* 0 except when binary special. */ -int use_mcell_ran4_; int nrn_xopen_broadcast_; void hoc_init(void) /* install constants and built-ins table */ @@ -255,7 +255,7 @@ void hoc_init(void) /* install constants and built-ins table */ } } - use_mcell_ran4_ = 0; + set_use_mcran4(false); nrn_xopen_broadcast_ = 255; extern void hoc_init_space(void); hoc_init_space(); diff --git a/src/oc/hocdec.h b/src/oc/hocdec.h index 113defaf53..944d4e9d72 100644 --- a/src/oc/hocdec.h +++ b/src/oc/hocdec.h @@ -224,10 +224,6 @@ struct HocParmUnits { /* units for symbol values */ #include "oc_ansi.h" -void* emalloc(size_t n); -void* ecalloc(size_t n, size_t size); -void* erealloc(void* ptr, size_t n); - extern Inst *hoc_progp, *hoc_progbase, *hoc_prog, *hoc_prog_parse_recover; extern Inst* hoc_pc; diff --git a/src/oc/memory.cpp b/src/oc/memory.cpp new file mode 100644 index 0000000000..44a6b93e51 --- /dev/null +++ b/src/oc/memory.cpp @@ -0,0 +1,112 @@ +#include "memory.hpp" + +#include +#include + +// For hoc_warning and hoc_execerror +#include "oc_ansi.h" + +#if HAVE_POSIX_MEMALIGN +#define HAVE_MEMALIGN 1 +#endif +#if defined(DARWIN) /* posix_memalign seems not to work on Darwin 10.6.2 */ +#undef HAVE_MEMALIGN +#endif +#if HAVE_MEMALIGN +#undef _XOPEN_SOURCE /* avoid warnings about redefining this */ +#define _XOPEN_SOURCE 600 +#endif + +static bool emalloc_error = false; + +void* hoc_Emalloc(std::size_t n) { /* check return from malloc */ + void* p = std::malloc(n); + if (p == nullptr) { + emalloc_error = true; + } + return p; +} + +void* hoc_Ecalloc(std::size_t n, std::size_t size) { /* check return from calloc */ + if (n == 0) { + return nullptr; + } + void* p = std::calloc(n, size); + if (p == nullptr) { + emalloc_error = true; + } + return p; +} + +void* hoc_Erealloc(void* ptr, std::size_t size) { /* check return from realloc */ + if (!ptr) { + return hoc_Emalloc(size); + } + void* p = std::realloc(ptr, size); + if (p == nullptr) { + std::free(ptr); + emalloc_error = true; + } + return p; +} + +void hoc_malchk(void) { + if (emalloc_error) { + emalloc_error = false; + hoc_execerror("out of memory", nullptr); + } +} + +void* emalloc(std::size_t n) { + void* p = hoc_Emalloc(n); + if (emalloc_error) { + hoc_malchk(); + } + return p; +} + +void* ecalloc(std::size_t n, std::size_t size) { + void* p = hoc_Ecalloc(n, size); + if (emalloc_error) { + hoc_malchk(); + } + return p; +} + +void* erealloc(void* ptr, std::size_t size) { + void* p = hoc_Erealloc(ptr, size); + if (emalloc_error) { + hoc_malchk(); + } + return p; +} + +void* nrn_cacheline_alloc(void** memptr, std::size_t size) { +#if HAVE_MEMALIGN + static bool memalign_is_working = true; + if (memalign_is_working) { + if (posix_memalign(memptr, 64, size) != 0) { + hoc_warning("posix_memalign not working, falling back to using malloc\n", nullptr); + memalign_is_working = false; + *memptr = hoc_Emalloc(size); + hoc_malchk(); + } + } else +#endif + { + *memptr = hoc_Emalloc(size); + hoc_malchk(); + } + return *memptr; +} + +void* nrn_cacheline_calloc(void** memptr, std::size_t nmemb, std::size_t size) { +#if HAVE_MEMALIGN + nrn_cacheline_alloc(memptr, nmemb * size); + std::memset(*memptr, 0, nmemb * size); +#else + *memptr = hoc_Ecalloc(nmemb, size); + hoc_malchk(); +#endif + return *memptr; +} diff --git a/src/oc/memory.hpp b/src/oc/memory.hpp new file mode 100644 index 0000000000..38960e0eba --- /dev/null +++ b/src/oc/memory.hpp @@ -0,0 +1,18 @@ +#pragma once + +// Some functions here are prepend with 'hoc_' but they are unrelated to hoc itself. +#include + +/* check return from malloc */ +void* hoc_Emalloc(std::size_t n); +void* hoc_Ecalloc(std::size_t n, std::size_t size); +void* hoc_Erealloc(void* ptr, std::size_t size); + +void hoc_malchk(void); + +void* emalloc(std::size_t n); +void* ecalloc(std::size_t n, std::size_t size); +void* erealloc(void* ptr, std::size_t size); + +void* nrn_cacheline_alloc(void** memptr, std::size_t size); +void* nrn_cacheline_calloc(void** memptr, std::size_t nmemb, std::size_t size); diff --git a/src/oc/mswinprt.cpp b/src/oc/mswinprt.cpp index 202a86db82..905aa101b2 100644 --- a/src/oc/mswinprt.cpp +++ b/src/oc/mswinprt.cpp @@ -148,11 +148,9 @@ void hoc_win_exec(void) { void hoc_winio_show(int b) {} -#if !defined(__MWERKS__) int getpid() { return 1; } -#endif void hoc_Plt() { TRY_GUI_REDIRECT_DOUBLE("plt", NULL); diff --git a/src/oc/nrnisaac.cpp b/src/oc/nrnisaac.cpp deleted file mode 100644 index 1aa0050339..0000000000 --- a/src/oc/nrnisaac.cpp +++ /dev/null @@ -1,39 +0,0 @@ -#include <../../nrnconf.h> -#include -#if HAVE_SYS_TYPES_H -#include -#endif -#include -#include -#include "hocdec.h" - -typedef struct isaac64_state Rng; - -void* nrnisaac_new(void) { - Rng* rng; - rng = (Rng*) hoc_Emalloc(sizeof(Rng)); - hoc_malchk(); - return (void*) rng; -} - -void nrnisaac_delete(void* v) { - free(v); -} - -void nrnisaac_init(void* v, unsigned long int seed) { - isaac64_init((Rng*) v, seed); -} - -double nrnisaac_dbl_pick(void* v) { - Rng* rng = (Rng*) v; - double x = isaac64_dbl32(rng); - /*printf("dbl %d %d %d %d %g\n", sizeof(ub8), sizeof(ub4), sizeof(ub2), sizeof(ub1), x);*/ - return x; -} - -uint32_t nrnisaac_uint32_pick(void* v) { - Rng* rng = (Rng*) v; - double x = isaac64_uint32(rng); - /*printf("uint32 %g\n", x);*/ - return x; -} diff --git a/src/oc/nrnisaac.h b/src/oc/nrnisaac.h deleted file mode 100644 index 053366684c..0000000000 --- a/src/oc/nrnisaac.h +++ /dev/null @@ -1,15 +0,0 @@ -#ifndef nrnisaac_h -#define nrnisaac_h - -#include <../../nrnconf.h> -#include - - -void* nrnisaac_new(void); -void nrnisaac_delete(void* rng); -void nrnisaac_init(void* rng, unsigned long int seed); -double nrnisaac_dbl_pick(void* rng); -uint32_t nrnisaac_uint32_pick(void* rng); - - -#endif diff --git a/src/oc/nrnmpiuse.h.in b/src/oc/nrnmpiuse.h.in index 7f17e5c0b2..4db1c75155 100755 --- a/src/oc/nrnmpiuse.h.in +++ b/src/oc/nrnmpiuse.h.in @@ -16,7 +16,4 @@ /* define if needed */ #undef ALWAYS_CALL_MPI_INIT -/* Number of times to retry a failed open */ -#undef FILE_OPEN_RETRY - #endif diff --git a/src/oc/nrnran123.cpp b/src/oc/nrnran123.cpp deleted file mode 100644 index d9db49351c..0000000000 --- a/src/oc/nrnran123.cpp +++ /dev/null @@ -1,134 +0,0 @@ -#include <../../nrnconf.h> -#include -#include -#include -#include -#include -#include - -static const double SHIFT32 = 1.0 / 4294967297.0; /* 1/(2^32 + 1) */ - -static philox4x32_key_t k = {{0}}; - -struct nrnran123_State { - philox4x32_ctr_t c; - philox4x32_ctr_t r; - char which_; -}; - -void nrnran123_set_globalindex(uint32_t gix) { - k.v[0] = gix; -} - -/* if one sets the global, one should reset all the stream sequences. */ -uint32_t nrnran123_get_globalindex() { - return k.v[0]; -} - -nrnran123_State* nrnran123_newstream(uint32_t id1, uint32_t id2) { - return nrnran123_newstream3(id1, id2, 0); -} -nrnran123_State* nrnran123_newstream3(uint32_t id1, uint32_t id2, uint32_t id3) { - nrnran123_State* s; - s = (nrnran123_State*) ecalloc(sizeof(nrnran123_State), 1); - s->c.v[1] = id3; - s->c.v[2] = id1; - s->c.v[3] = id2; - nrnran123_setseq(s, 0, 0); - return s; -} - -void nrnran123_deletestream(nrnran123_State* s) { - free(s); -} - -void nrnran123_getseq(nrnran123_State* s, uint32_t* seq, char* which) { - *seq = s->c.v[0]; - *which = s->which_; -} - -void nrnran123_setseq(nrnran123_State* s, uint32_t seq, char which) { - if (which > 3 || which < 0) { - s->which_ = 0; - } else { - s->which_ = which; - } - s->c.v[0] = seq; - s->r = philox4x32(s->c, k); -} - -void nrnran123_getids(nrnran123_State* s, uint32_t* id1, uint32_t* id2) { - *id1 = s->c.v[2]; - *id2 = s->c.v[3]; -} - -void nrnran123_getids3(nrnran123_State* s, uint32_t* id1, uint32_t* id2, uint32_t* id3) { - *id3 = s->c.v[1]; - *id1 = s->c.v[2]; - *id2 = s->c.v[3]; -} - -uint32_t nrnran123_ipick(nrnran123_State* s) { - uint32_t rval; - char which = s->which_; - assert(which < 4); - rval = s->r.v[which++]; - if (which > 3) { - which = 0; - s->c.v[0]++; - s->r = philox4x32(s->c, k); - } - s->which_ = which; - return rval; -} - -double nrnran123_dblpick(nrnran123_State* s) { - return nrnran123_uint2dbl(nrnran123_ipick(s)); -} - -double nrnran123_negexp(nrnran123_State* s) { - /* min 2.3283064e-10 to max 22.18071 */ - return -log(nrnran123_dblpick(s)); -} - -/* At cost of a cached value we could compute two at a time. */ -/* But that would make it difficult to transfer to coreneuron for t > 0 */ -double nrnran123_normal(nrnran123_State* s) { - double w, x, y; - double u1, u2; - do { - u1 = nrnran123_dblpick(s); - u2 = nrnran123_dblpick(s); - u1 = 2. * u1 - 1.; - u2 = 2. * u2 - 1.; - w = (u1 * u1) + (u2 * u2); - } while (w > 1); - - y = sqrt((-2. * log(w)) / w); - x = u1 * y; - return x; -} - -nrnran123_array4x32 nrnran123_iran(uint32_t seq, uint32_t id1, uint32_t id2) { - return nrnran123_iran3(seq, id1, id2, 0); -} -nrnran123_array4x32 nrnran123_iran3(uint32_t seq, uint32_t id1, uint32_t id2, uint32_t id3) { - nrnran123_array4x32 a; - philox4x32_ctr_t c; - c.v[0] = seq; - c.v[1] = id3; - c.v[2] = id1; - c.v[3] = id2; - philox4x32_ctr_t r = philox4x32(c, k); - a.v[0] = r.v[0]; - a.v[1] = r.v[1]; - a.v[2] = r.v[2]; - a.v[3] = r.v[3]; - return a; -} - -double nrnran123_uint2dbl(uint32_t u) { - /* 0 to 2^32-1 transforms to double value in open (0,1) interval */ - /* min 2.3283064e-10 to max (1 - 2.3283064e-10) */ - return ((double) u + 1.0) * SHIFT32; -} diff --git a/src/oc/nrnran123.h b/src/oc/nrnran123.h deleted file mode 100644 index 0a0acb6a71..0000000000 --- a/src/oc/nrnran123.h +++ /dev/null @@ -1,62 +0,0 @@ -#ifndef nrnran123_h -#define nrnran123_h - -/* interface to Random123 */ -/* http://www.thesalmons.org/john/random123/papers/random123sc11.pdf */ - -/* -The 4x32 generators utilize a uint32x4 counter and uint32x4 key to transform -into an almost cryptographic quality uint32x4 random result. -There are many possibilites for balancing the sharing of the internal -state instances while reserving a uint32 counter for the stream sequence -and reserving other portions of the counter vector for stream identifiers -and global index used by all streams. - -We currently provide a single instance by default in which the policy is -to use the 0th counter uint32 as the stream sequence, words 2, 3 and 4 as the -stream identifier, and word 0 of the key as the global index. Unused words -are constant uint32 0. - -It is also possible to use Random123 directly without reference to this -interface. See Random123-1.02/docs/html/index.html -of the full distribution available from -http://www.deshawresearch.com/resources_random123.html -*/ - -#include - - -typedef struct nrnran123_State nrnran123_State; - -typedef struct nrnran123_array4x32 { - uint32_t v[4]; -} nrnran123_array4x32; - -/* global index. eg. run number */ -/* all generator instances share this global index */ -extern void nrnran123_set_globalindex(uint32_t gix); -extern uint32_t nrnran123_get_globalindex(); - -/* minimal data stream */ -extern nrnran123_State* nrnran123_newstream(uint32_t id1, uint32_t id2); -extern nrnran123_State* nrnran123_newstream3(uint32_t id1, uint32_t id2, uint32_t id3); -extern void nrnran123_deletestream(nrnran123_State*); -extern void nrnran123_getseq(nrnran123_State*, uint32_t* seq, char* which); -extern void nrnran123_setseq(nrnran123_State*, uint32_t seq, char which); -extern void nrnran123_getids(nrnran123_State*, uint32_t* id1, uint32_t* id2); -extern void nrnran123_getids3(nrnran123_State*, uint32_t* id1, uint32_t* id2, uint32_t* id3); - -extern double nrnran123_negexp(nrnran123_State*); /* mean 1.0 */ -extern uint32_t nrnran123_ipick(nrnran123_State*); /* uniform 0 to 2^32-1 */ -extern double nrnran123_dblpick(nrnran123_State*); /* uniform open interval (0,1)*/ -/* nrnran123_dblpick minimum value is 2.3283064e-10 and max value is 1-min */ - -/* nrnran123_negexp min value is 2.3283064e-10, max is 22.18071 */ -extern double nrnran123_normal(nrnran123_State*); /* mean 0.0, std 1.0 */ - -/* more fundamental (stateless) (though the global index is still used) */ -extern nrnran123_array4x32 nrnran123_iran(uint32_t seq, uint32_t id1, uint32_t id2); -extern nrnran123_array4x32 nrnran123_iran3(uint32_t seq, uint32_t id1, uint32_t id2, uint32_t id3); -extern double nrnran123_uint2dbl(uint32_t); - -#endif diff --git a/src/oc/oc_ansi.h b/src/oc/oc_ansi.h index 399a86e408..34cc37062f 100644 --- a/src/oc/oc_ansi.h +++ b/src/oc/oc_ansi.h @@ -8,6 +8,8 @@ #include #include #include + +#include "memory.hpp" /** * \dir * \brief HOC Interpreter @@ -49,9 +51,6 @@ void ivoc_help(const char*); Symbol* hoc_lookup(const char*); -void* hoc_Ecalloc(std::size_t nmemb, std::size_t size); -void* hoc_Emalloc(size_t size); -void hoc_malchk(); [[noreturn]] void hoc_execerror(const char*, const char*); [[noreturn]] void hoc_execerr_ext(const char* fmt, ...); char* hoc_object_name(Object*); @@ -320,10 +319,6 @@ void hoc_obj_set(int i, Object*); void nrn_hoc_lock(); void nrn_hoc_unlock(); -void* hoc_Erealloc(void* ptr, std::size_t size); - -void* nrn_cacheline_alloc(void** memptr, std::size_t size); -void* nrn_cacheline_calloc(void** memptr, std::size_t nmemb, std::size_t size); [[noreturn]] void nrn_exit(int); void hoc_free_list(Symlist**); int hoc_errno_check(); diff --git a/src/oc/oc_mcran4.cpp b/src/oc/oc_mcran4.cpp new file mode 100644 index 0000000000..128ccd7d04 --- /dev/null +++ b/src/oc/oc_mcran4.cpp @@ -0,0 +1,41 @@ +#include "hocdec.h" +#include "mcran4.h" + +int use_mcell_ran4_; + +void set_use_mcran4(bool value) { + use_mcell_ran4_ = value ? 1 : 0; +} + +bool use_mcran4() { + return use_mcell_ran4_ != 0; +} + +void hoc_mcran4() { + uint32_t idx; + double* xidx; + double x; + xidx = hoc_pgetarg(1); + idx = (uint32_t) (*xidx); + x = mcell_ran4a(&idx); + *xidx = idx; + hoc_ret(); + hoc_pushx(x); +} +void hoc_mcran4init() { + double prev = mcell_lowindex(); + if (ifarg(1)) { + uint32_t idx = (uint32_t) chkarg(1, 0., 4294967295.); + mcell_ran4_init(idx); + } + hoc_ret(); + hoc_pushx(prev); +} +void hoc_usemcran4() { + double prev = (double) use_mcell_ran4_; + if (ifarg(1)) { + use_mcell_ran4_ = (int) chkarg(1, 0., 1.); + } + hoc_ret(); + hoc_pushx(prev); +} diff --git a/src/oc/oc_mcran4.hpp b/src/oc/oc_mcran4.hpp new file mode 100644 index 0000000000..d275e3e9bf --- /dev/null +++ b/src/oc/oc_mcran4.hpp @@ -0,0 +1,5 @@ +void set_use_mcran4(bool value); +bool use_mcran4(); +void hoc_mcran4(); +void hoc_mcran4init(); +void hoc_usemcran4(); diff --git a/src/oc/ocfunc.h b/src/oc/ocfunc.h index fca598c0fa..61a620aa44 100644 --- a/src/oc/ocfunc.h +++ b/src/oc/ocfunc.h @@ -32,7 +32,6 @@ extern void hoc_neuronhome(void), hoc_Execerror(void); extern void hoc_sscanf(void), hoc_save_session(void), hoc_print_session(void); extern void hoc_Chdir(void), hoc_getcwd(void), hoc_Symbol_units(void), hoc_stdout(void); extern void hoc_name_declared(void), hoc_unix_mac_pc(void), hoc_show_winio(void); -extern void hoc_usemcran4(void), hoc_mcran4(void), hoc_mcran4init(void); extern void hoc_nrn_load_dll(void), hoc_nrnversion(void), hoc_object_pushed(void); extern void hoc_mallinfo(void); extern void hoc_Setcolor(void); diff --git a/src/oc/scoprand.cpp b/src/oc/scoprand.cpp index 7f08ae8814..70651f5608 100644 --- a/src/oc/scoprand.cpp +++ b/src/oc/scoprand.cpp @@ -23,7 +23,8 @@ static char RCSid[] = "random.cpp,v 1.4 1999/01/04 12:46:49 hines Exp"; #endif #include -#include +#include "oc_mcran4.hpp" +#include "mcran4.h" #include "scoplib.h" static uint32_t value = 1; @@ -58,8 +59,7 @@ static uint32_t value = 1; *--------------------------------------------------------------------------- */ double scop_random(void) { - extern int use_mcell_ran4_; - if (use_mcell_ran4_) { + if (use_mcran4()) { /*perhaps 4 times slower but much higher quality*/ return mcell_ran4a(&value); } else { diff --git a/src/oc/symbol.cpp b/src/oc/symbol.cpp index 7257a91aa3..605a7ba2a6 100644 --- a/src/oc/symbol.cpp +++ b/src/oc/symbol.cpp @@ -2,17 +2,6 @@ /* /local/src/master/nrn/src/oc/symbol.cpp,v 1.9 1999/02/25 18:01:58 hines Exp */ /* version 7.2.1 2-jan-89 */ -#if HAVE_POSIX_MEMALIGN -#define HAVE_MEMALIGN 1 -#endif -#if defined(DARWIN) /* posix_memalign seems not to work on Darwin 10.6.2 */ -#undef HAVE_MEMALIGN -#endif -#if HAVE_MEMALIGN -#undef _XOPEN_SOURCE /* avoid warnings about redefining this */ -#define _XOPEN_SOURCE 600 -#endif - #include "hoc.h" #include "hocdec.h" #include "hoclist.h" @@ -173,96 +162,6 @@ void hoc_link_symbol(Symbol* sp, Symlist* list) { sp->next = nullptr; } -static int emalloc_error = 0; - -void hoc_malchk(void) { - if (emalloc_error) { - emalloc_error = 0; - execerror("out of memory", nullptr); - } -} - -void* hoc_Emalloc(size_t n) { /* check return from malloc */ - void* p = malloc(n); - if (p == nullptr) - emalloc_error = 1; - return p; -} - -void* emalloc(size_t n) { - void* p = hoc_Emalloc(n); - if (emalloc_error) { - hoc_malchk(); - } - return p; -} - -void* hoc_Ecalloc(size_t n, size_t size) { /* check return from calloc */ - if (n == 0) { - return nullptr; - } - void* p = calloc(n, size); - if (p == nullptr) - emalloc_error = 1; - return p; -} - -void* ecalloc(size_t n, size_t size) { - void* p = hoc_Ecalloc(n, size); - if (emalloc_error) { - hoc_malchk(); - } - return p; -} - -void* nrn_cacheline_alloc(void** memptr, size_t size) { -#if HAVE_MEMALIGN - static int memalign_is_working = 1; - if (memalign_is_working) { - if (posix_memalign(memptr, 64, size) != 0) { - fprintf(stderr, "posix_memalign not working, falling back to using malloc\n"); - memalign_is_working = 0; - *memptr = hoc_Emalloc(size); - hoc_malchk(); - } - } else -#endif - *memptr = hoc_Emalloc(size); - hoc_malchk(); - return *memptr; -} - -void* nrn_cacheline_calloc(void** memptr, size_t nmemb, size_t size) { -#if HAVE_MEMALIGN - nrn_cacheline_alloc(memptr, nmemb * size); - memset(*memptr, 0, nmemb * size); -#else - *memptr = hoc_Ecalloc(nmemb, size); - hoc_malchk(); -#endif - return *memptr; -} - -void* hoc_Erealloc(void* ptr, size_t size) { /* check return from realloc */ - if (!ptr) { - return hoc_Emalloc(size); - } - void* p = realloc(ptr, size); - if (p == nullptr) { - free(ptr); - emalloc_error = 1; - } - return p; -} - -void* erealloc(void* ptr, size_t size) { - void* p = hoc_Erealloc(ptr, size); - if (emalloc_error) { - hoc_malchk(); - } - return p; -} - void hoc_free_symspace(Symbol* s1) { /* frees symbol space. Marks it UNDEF */ if (s1 && s1->cpublic != 2) { switch (s1->type) { diff --git a/src/sparse13/CMakeLists.txt b/src/sparse13/CMakeLists.txt index 69f9e1a955..b2adc31dc2 100644 --- a/src/sparse13/CMakeLists.txt +++ b/src/sparse13/CMakeLists.txt @@ -1,15 +1,3 @@ -add_library( - sparse13 STATIC - spalloc.cpp - spbuild.cpp - spfactor.cpp - spoutput.cpp - spsolve.cpp - sputils.cpp - cspalloc.cpp - cspbuild.cpp - cspfactor.cpp - cspoutput.cpp - cspsolve.cpp - csputils.cpp) +add_library(sparse13 STATIC spalloc.cpp spbuild.cpp spfactor.cpp spoutput.cpp spsolve.cpp + sputils.cpp) set_property(TARGET sparse13 PROPERTY POSITION_INDEPENDENT_CODE ON) diff --git a/src/sparse13/cspalloc.cpp b/src/sparse13/cspalloc.cpp deleted file mode 100644 index 29b7ee44c8..0000000000 --- a/src/sparse13/cspalloc.cpp +++ /dev/null @@ -1,2 +0,0 @@ -#define cmplx_spPrefix -#include "spalloc.cpp" diff --git a/src/sparse13/cspbuild.cpp b/src/sparse13/cspbuild.cpp deleted file mode 100644 index b002c1e5c2..0000000000 --- a/src/sparse13/cspbuild.cpp +++ /dev/null @@ -1,2 +0,0 @@ -#define cmplx_spPrefix -#include "spbuild.cpp" diff --git a/src/sparse13/cspfactor.cpp b/src/sparse13/cspfactor.cpp deleted file mode 100644 index 2a9d91aaaa..0000000000 --- a/src/sparse13/cspfactor.cpp +++ /dev/null @@ -1,2 +0,0 @@ -#define cmplx_spPrefix -#include "spfactor.cpp" diff --git a/src/sparse13/cspmatrix.h b/src/sparse13/cspmatrix.h deleted file mode 100644 index 7200d07b55..0000000000 --- a/src/sparse13/cspmatrix.h +++ /dev/null @@ -1,2 +0,0 @@ -#define cmplx_spPrefix -#include "spmatrix.h" diff --git a/src/sparse13/cspoutput.cpp b/src/sparse13/cspoutput.cpp deleted file mode 100644 index 59325c189a..0000000000 --- a/src/sparse13/cspoutput.cpp +++ /dev/null @@ -1,2 +0,0 @@ -#define cmplx_spPrefix -#include "spoutput.cpp" diff --git a/src/sparse13/cspredef.h b/src/sparse13/cspredef.h deleted file mode 100644 index 2df7687583..0000000000 --- a/src/sparse13/cspredef.h +++ /dev/null @@ -1,49 +0,0 @@ -/* mostly generated from -cat temp | /usr/bin/tr -cs "[:alpha:]" "[\n*]" | sort | uniq | grep '^sp' > sp_redef.h -where temp is the last part of spmatrix.h -*/ -#define spClear cmplx_spClear -#define spCondition cmplx_spCondition -#define spCreate cmplx_spCreate -#define spDeleteRowAndCol cmplx_spDeleteRowAndCol -#define spDestroy cmplx_spDestroy -#define spDeterminant cmplx_spDeterminant -#define spElementCount cmplx_spElementCount -#define spError cmplx_spError -#define spFactor cmplx_spFactor -#define spFileMatrix cmplx_spFileMatrix -#define spFileStats cmplx_spFileStats -#define spFileVector cmplx_spFileVector -#define spFillinCount cmplx_spFillinCount -#define spGetAdmittance cmplx_spGetAdmittance -#define spGetElement cmplx_spGetElement -#define spGetInitInfo cmplx_spGetInitInfo -#define spGetOnes cmplx_spGetOnes -#define spGetQuad cmplx_spGetQuad -#define spGetSize cmplx_spGetSize -#define spInitialize cmplx_spInitialize -#define spInstallInitInfo cmplx_spInstallInitInfo -#define spLargestElement cmplx_spLargestElement -#define spMNA_Preorder cmplx_spMNA_Preorder -#define spMultTransposed cmplx_spMultTransposed -#define spMultiply cmplx_spMultiply -#define spNorm cmplx_spNorm -#define spOrderAndFactor cmplx_spOrderAndFactor -#define spPartition cmplx_spPartition -#define spPrint cmplx_spPrint -#define spPseudoCondition cmplx_spPseudoCondition -#define spRoundoff cmplx_spRoundoff -#define spScale cmplx_spScale -#define spSetComplex cmplx_spSetComplex -#define spSetReal cmplx_spSetReal -#define spSolve cmplx_spSolve -#define spSolveTransposed cmplx_spSolveTransposed -#define spStripFills cmplx_spStripFills -#define spWhereSingular cmplx_spWhereSingular -#define spcGetFillin cmplx_spcGetFillin -#define spcGetElement cmplx_spcGetElement -#define spcLinkRows cmplx_spcLinkRows -#define spcFindElementInCol cmplx_spcFindElementInCol -#define spcCreateElement cmplx_spcCreateElement -#define spcRowExchange cmplx_spcRowExchange -#define spcColExchange cmplx_spcColExchange diff --git a/src/sparse13/cspsolve.cpp b/src/sparse13/cspsolve.cpp deleted file mode 100644 index f1aeba7d67..0000000000 --- a/src/sparse13/cspsolve.cpp +++ /dev/null @@ -1,2 +0,0 @@ -#define cmplx_spPrefix -#include "spsolve.cpp" diff --git a/src/sparse13/csputils.cpp b/src/sparse13/csputils.cpp deleted file mode 100644 index 0f0363b8ae..0000000000 --- a/src/sparse13/csputils.cpp +++ /dev/null @@ -1,2 +0,0 @@ -#define cmplx_spPrefix -#include "sputils.cpp"