References for 'BlueBrain/nmodl#1344'.

at_time/coreneuron/solver/crout.hpp

@@ -0,0 +1,207 @@
+/*
+ * Copyright 2023 Blue Brain Project, EPFL.
+ * See the top-level LICENSE file for details.
+ *
+ * SPDX-License-Identifier: Apache-2.0
+ */
+
+#pragma once
+
+/**
+ * \dir
+ * \brief Solver for a system of linear equations : Crout matrix decomposition
+ *
+ * \file
+ * \brief Implementation of Crout matrix decomposition (LU decomposition) followed by
+ * Forward/Backward substitution: Implementation details : (Legacy code) nrn / scopmath / crout.c
+ */
+
+#include <Eigen/Core>
+#include <cmath>
+
+#if defined(CORENEURON_ENABLE_GPU) && !defined(DISABLE_OPENACC)
+#include "coreneuron/utils/offload.hpp"
+#endif
+
+namespace nmodl {
+namespace crout {
+
+/**
+ * \brief Crout matrix decomposition : in-place LU Decomposition of matrix a.
+ *
+ * Implementation details : (Legacy code) nrn / scopmath / crout.c
+ *
+ * Returns: 0 if no error; -1 if matrix is singular or ill-conditioned
+ */
+#if defined(CORENEURON_ENABLE_GPU) && !defined(DISABLE_OPENACC)
+nrn_pragma_acc(routine seq)
+nrn_pragma_omp(declare target)
+#endif
+template <typename T>
+EIGEN_DEVICE_FUNC inline int Crout(int n, T* const a, int* const perm, double* const rowmax) {
+    // roundoff is the minimal value for a pivot element without its being considered too close to
+    // zero
+    double roundoff = 1.e-20;
+    int i, j, k, r, pivot, irow, save_i = 0, krow;
+    T sum, equil_1, equil_2;
+
+    /* Initialize permutation and rowmax vectors */
+
+    for (i = 0; i < n; i++) {
+        perm[i] = i;
+        k = 0;
+        for (j = 1; j < n; j++)
+            if (std::fabs(a[i * n + j]) > std::fabs(a[i * n + k]))
+                k = j;
+        rowmax[i] = a[i * n + k];
+    }
+
+    /* Loop over rows and columns r */
+
+    for (r = 0; r < n; r++) {
+        /*
+         * Operate on rth column.  This produces the lower triangular matrix
+         * of terms needed to transform the constant vector.
+         */
+
+        for (i = r; i < n; i++) {
+            sum = 0.0;
+            irow = perm[i];
+            for (k = 0; k < r; k++) {
+                krow = perm[k];
+                sum += a[irow * n + k] * a[krow * n + r];
+            }
+            a[irow * n + r] -= sum;
+        }
+
+        /* Find row containing the pivot in the rth column */
+
+        pivot = perm[r];
+        equil_1 = std::fabs(a[pivot * n + r] / rowmax[pivot]);
+        for (i = r + 1; i < n; i++) {
+            irow = perm[i];
+            equil_2 = std::fabs(a[irow * n + r] / rowmax[irow]);
+            if (equil_2 > equil_1) {
+                /* make irow the new pivot row */
+
+                pivot = irow;
+                save_i = i;
+                equil_1 = equil_2;
+            }
+        }
+
+        /* Interchange entries in permutation vector if necessary */
+
+        if (pivot != perm[r]) {
+            perm[save_i] = perm[r];
+            perm[r] = pivot;
+        }
+
+        /* Check that pivot element is not too small */
+
+        if (std::fabs(a[pivot * n + r]) < roundoff)
+            return -1;
+
+        /*
+         * Operate on row in rth position.  This produces the upper
+         * triangular matrix whose diagonal elements are assumed to be unity.
+         * This matrix is used in the back substitution algorithm.
+         */
+
+        for (j = r + 1; j < n; j++) {
+            sum = 0.0;
+            for (k = 0; k < r; k++) {
+                krow = perm[k];
+                sum += a[pivot * n + k] * a[krow * n + j];
+            }
+            a[pivot * n + j] = (a[pivot * n + j] - sum) / a[pivot * n + r];
+        }
+    }
+    return 0;
+}
+#if defined(CORENEURON_ENABLE_GPU) && !defined(DISABLE_OPENACC)
+nrn_pragma_omp(end declare target)
+#endif
+
+/**
+ * \brief Crout matrix decomposition : Forward/Backward substitution.
+ *
+ * Implementation details : (Legacy code) nrn / scopmath / crout.c
+ *
+ * Returns: no return variable
+ */
+#define y_(arg) p[y[arg]]
+#define b_(arg) b[arg]
+#if defined(CORENEURON_ENABLE_GPU) && !defined(DISABLE_OPENACC)
+nrn_pragma_acc(routine seq)
+nrn_pragma_omp(declare target)
+#endif
+template <typename T>
+EIGEN_DEVICE_FUNC inline int solveCrout(int n,
+                                        T const* const a,
+                                        T const* const b,
+                                        T* const p,
+                                        int const* const perm,
+                                        int const* const y = nullptr) {
+    int i, j, pivot;
+    T sum;
+
+    /* Perform forward substitution with pivoting */
+    if (y) {
+        for (i = 0; i < n; i++) {
+            pivot = perm[i];
+            sum = 0.0;
+            for (j = 0; j < i; j++)
+                sum += a[pivot * n + j] * (y_(j));
+            y_(i) = (b_(pivot) - sum) / a[pivot * n + i];
+        }
+
+        /*
+         * Note that the y vector is already in the correct order for back
+         * substitution.  Perform back substitution, pivoting the matrix but not
+         * the y vector.  There is no need to divide by the diagonal element as
+         * this is assumed to be unity.
+         */
+
+        for (i = n - 1; i >= 0; i--) {
+            pivot = perm[i];
+            sum = 0.0;
+            for (j = i + 1; j < n; j++)
+                sum += a[pivot * n + j] * (y_(j));
+            y_(i) -= sum;
+        }
+    } else {
+        for (i = 0; i < n; i++) {
+            pivot = perm[i];
+            sum = 0.0;
+            for (j = 0; j < i; j++)
+                sum += a[pivot * n + j] * (p[j]);
+            p[i] = (b_(pivot) - sum) / a[pivot * n + i];
+        }
+
+        /*
+         * Note that the y vector is already in the correct order for back
+         * substitution.  Perform back substitution, pivoting the matrix but not
+         * the y vector.  There is no need to divide by the diagonal element as
+         * this is assumed to be unity.
+         */
+
+        for (i = n - 1; i >= 0; i--) {
+            pivot = perm[i];
+            sum = 0.0;
+            for (j = i + 1; j < n; j++)
+                sum += a[pivot * n + j] * (p[j]);
+            p[i] -= sum;
+        }
+    }
+    return 0;
+}
+#if defined(CORENEURON_ENABLE_GPU) && !defined(DISABLE_OPENACC)
+nrn_pragma_omp(end declare target)
+#endif
+
+#undef y_
+#undef b_
+
+}  // namespace crout
+}  // namespace nmodl