diff options
Diffstat (limited to 'src/lib/netlist/solver/mat_cr.h')
-rw-r--r-- | src/lib/netlist/solver/mat_cr.h | 114 |
1 files changed, 75 insertions, 39 deletions
diff --git a/src/lib/netlist/solver/mat_cr.h b/src/lib/netlist/solver/mat_cr.h index 340c662752f..c39af1b2f5f 100644 --- a/src/lib/netlist/solver/mat_cr.h +++ b/src/lib/netlist/solver/mat_cr.h @@ -12,38 +12,79 @@ #include <algorithm> #include "plib/pconfig.h" +#include "plib/palloc.h" -template<int storage_N> +template<std::size_t N, typename C = uint16_t, typename T = double> struct mat_cr_t { - unsigned nz_num; - unsigned ia[storage_N + 1]; - unsigned ja[storage_N * storage_N]; - unsigned diag[storage_N]; /* n */ + typedef C index_type; + typedef T value_type; - template<typename T> - void mult_vec(const T * RESTRICT A, const T * RESTRICT x, T * RESTRICT res) + C diag[N]; // diagonal index pointer n + C ia[N+1]; // row index pointer n + 1 + C ja[N*N]; // column index array nz_num, initially (n * n) + T A[N*N]; // Matrix elements nz_num, initially (n * n) + + std::size_t size; + std::size_t nz_num; + + mat_cr_t(const std::size_t n) + : size(n) + , nz_num(0) + { +#if 0 +#if 0 + ia = plib::palloc_array<C>(n + 1); + ja = plib::palloc_array<C>(n * n); + diag = plib::palloc_array<C>(n); +#else + diag = plib::palloc_array<C>(n + (n + 1) + n * n); + ia = diag + n; + ja = ia + (n+1); + A = plib::palloc_array<T>(n * n); +#endif +#endif + } + + ~mat_cr_t() + { +#if 0 + plib::pfree_array(diag); +#if 0 + plib::pfree_array(ia); + plib::pfree_array(ja); +#endif + plib::pfree_array(A); +#endif + } + + void set_scalar(const T scalar) + { + for (std::size_t i=0, e=nz_num; i<e; i++) + A[i] = scalar; + } + + void mult_vec(const T * RESTRICT x, T * RESTRICT res) { /* * res = A * x */ - unsigned i = 0; - unsigned k = 0; - const unsigned oe = nz_num; + std::size_t i = 0; + std::size_t k = 0; + const std::size_t oe = nz_num; while (k < oe) { T tmp = 0.0; - const unsigned e = ia[i+1]; + const std::size_t e = ia[i+1]; for (; k < e; k++) tmp += A[k] * x[ja[k]]; res[i++] = tmp; } } - template<typename T> - void incomplete_LU_factorization(const T * RESTRICT A, T * RESTRICT LU) + void incomplete_LU_factorization(T * RESTRICT LU) { /* * incomplete LU Factorization according to http://de.wikipedia.org/wiki/ILU-Zerlegung @@ -52,28 +93,28 @@ struct mat_cr_t * */ - const unsigned lnz = nz_num; + const std::size_t lnz = nz_num; - for (unsigned k = 0; k < lnz; k++) + for (std::size_t k = 0; k < lnz; k++) LU[k] = A[k]; - for (unsigned i = 1; ia[i] < lnz; i++) // row i + for (std::size_t i = 1; ia[i] < lnz; i++) // row i { - const unsigned iai1 = ia[i + 1]; - const unsigned pke = diag[i]; - for (unsigned pk = ia[i]; pk < pke; pk++) // all columns left of diag in row i + const std::size_t iai1 = ia[i + 1]; + const std::size_t pke = diag[i]; + for (std::size_t pk = ia[i]; pk < pke; pk++) // all columns left of diag in row i { // pk == (i, k) - const unsigned k = ja[pk]; - const unsigned iak1 = ia[k + 1]; + const std::size_t k = ja[pk]; + const std::size_t iak1 = ia[k + 1]; const T LUpk = LU[pk] = LU[pk] / LU[diag[k]]; - unsigned pt = ia[k]; + std::size_t pt = ia[k]; - for (unsigned pj = pk + 1; pj < iai1; pj++) // pj = (i, j) + for (std::size_t pj = pk + 1; pj < iai1; pj++) // pj = (i, j) { // we can assume that within a row ja increases continuously */ - const unsigned ej = ja[pj]; + const std::size_t ej = ja[pj]; while (ja[pt] < ej && pt < iak1) pt++; if (pt < iak1 && ja[pt] == ej) @@ -83,7 +124,6 @@ struct mat_cr_t } } - template<typename T> void solveLUx (const T * RESTRICT LU, T * RESTRICT r) { /* @@ -96,7 +136,7 @@ struct mat_cr_t * * ==> LUx = r * - * ==> Ux = L?????r = w + * ==> Ux = L⁻¹ r = w * * ==> r = Lw * @@ -108,32 +148,28 @@ struct mat_cr_t * */ - unsigned i; - - for (i = 1; ia[i] < nz_num; i++ ) + for (std::size_t i = 1; ia[i] < nz_num; ++i ) { T tmp = 0.0; - const unsigned j1 = ia[i]; - const unsigned j2 = diag[i]; + const std::size_t j1 = ia[i]; + const std::size_t j2 = diag[i]; - for (unsigned j = j1; j < j2; j++ ) + for (std::size_t j = j1; j < j2; ++j ) tmp += LU[j] * r[ja[j]]; r[i] -= tmp; } // i now is equal to n; - for (; 0 < i; i-- ) + for (std::size_t i = size; i-- > 0; ) { - const unsigned im1 = i - 1; T tmp = 0.0; - const unsigned j1 = diag[im1] + 1; - const unsigned j2 = ia[im1+1]; - for (unsigned j = j1; j < j2; j++ ) + const std::size_t di = diag[i]; + const std::size_t j2 = ia[i+1]; + for (std::size_t j = di + 1; j < j2; j++ ) tmp += LU[j] * r[ja[j]]; - r[im1] = (r[im1] - tmp) / LU[diag[im1]]; + r[i] = (r[i] - tmp) / LU[di]; } } - }; #endif /* MAT_CR_H_ */ |