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.h136
1 files changed, 136 insertions, 0 deletions
diff --git a/src/lib/netlist/solver/mat_cr.h b/src/lib/netlist/solver/mat_cr.h
new file mode 100644
index 00000000000..ea72884c1de
--- /dev/null
+++ b/src/lib/netlist/solver/mat_cr.h
@@ -0,0 +1,136 @@
+// license:GPL-2.0+
+// copyright-holders:Couriersud
+/*
+ * mat_cr.h
+ *
+ * Compressed row format matrices
+ *
+ */
+
+#ifndef MAT_CR_H_
+#define MAT_CR_H_
+
+#include <algorithm>
+#include "plib/pconfig.h"
+
+template<int _storage_N>
+struct mat_cr_t
+{
+ unsigned nz_num;
+ unsigned ia[_storage_N + 1];
+ unsigned ja[_storage_N * _storage_N];
+ unsigned diag[_storage_N]; /* n */
+
+ template<typename T>
+ void mult_vec(const T * RESTRICT A, const T * RESTRICT x, T * RESTRICT res)
+ {
+ /*
+ * res = A * x
+ */
+
+ unsigned i = 0;
+ unsigned k = 0;
+ const unsigned oe = nz_num;
+
+ while (k < oe)
+ {
+ double tmp = 0.0;
+ const unsigned e = ia[i+1];
+ for (; k < e; k++)
+ tmp += A[k] * x[ja[k]];
+ res[i++] = tmp;
+ }
+ }
+
+ void incomplete_LU_factorization(const nl_double * RESTRICT A, nl_double * RESTRICT LU)
+ {
+ /*
+ * incomplete LU Factorization according to http://de.wikipedia.org/wiki/ILU-Zerlegung
+ *
+ * Result is stored in matrix LU
+ *
+ */
+
+ const unsigned lnz = nz_num;
+
+ for (unsigned k = 0; k < lnz; k++)
+ LU[k] = A[k];
+
+ for (unsigned i = 1; ia[i] < lnz; i++) // row i
+ {
+ const unsigned iai1 = ia[i + 1];
+ for (unsigned pk = ia[i]; pk < diag[i]; pk++) // all columns left of diag in row i
+ {
+ // pk == (i, k)
+ const unsigned k = ja[pk];
+ const unsigned iak1 = ia[k + 1];
+ const double LUpk = LU[pk] = LU[pk] / LU[diag[k]];
+
+ unsigned pt = ia[k];
+
+ for (unsigned pj = pk + 1; pj < iai1; pj++) // pj = (i, j)
+ {
+ // we can assume that within a row ja increases continuously */
+ const unsigned ej = ja[pj];
+ while (ja[pt] < ej && pt < iak1)
+ pt++;
+ if (pt < iak1 && ja[pt] == ej)
+ LU[pj] = LU[pj] - LUpk * LU[pt];
+ }
+ }
+ }
+ }
+
+ void solveLUx (const nl_double * RESTRICT LU, nl_double * RESTRICT r)
+ {
+ /*
+ * Solve a linear equation Ax = r
+ * where
+ * A = L*U
+ *
+ * L unit lower triangular
+ * U upper triangular
+ *
+ * ==> LUx = r
+ *
+ * ==> Ux = L?????r = w
+ *
+ * ==> r = Lw
+ *
+ * This can be solved for w using backwards elimination in L.
+ *
+ * Now Ux = w
+ *
+ * This can be solved for x using backwards elimination in U.
+ *
+ */
+
+ unsigned i;
+
+ for (i = 1; ia[i] < nz_num; i++ )
+ {
+ double tmp = 0.0;
+ const unsigned j1 = ia[i];
+ const unsigned j2 = diag[i];
+
+ for (unsigned j = j1; j < j2; j++ )
+ tmp += LU[j] * r[ja[j]];
+
+ r[i] -= tmp;
+ }
+ // i now is equal to n;
+ for (; 0 < i; i-- )
+ {
+ const unsigned im1 = i - 1;
+ double tmp = 0.0;
+ const unsigned j1 = diag[im1] + 1;
+ const unsigned j2 = ia[im1+1];
+ for (unsigned j = j1; j < j2; j++ )
+ tmp += LU[j] * r[ja[j]];
+ r[im1] = (r[im1] - tmp) / LU[diag[im1]];
+ }
+ }
+
+};
+
+#endif /* MAT_CR_H_ */