summaryrefslogtreecommitdiffstatshomepage
path: root/src/lib/netlist/solver/mat_cr.h
diff options
context:
space:
mode:
Diffstat (limited to 'src/lib/netlist/solver/mat_cr.h')
-rw-r--r--src/lib/netlist/solver/mat_cr.h114
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_ */