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.h175
1 files changed, 0 insertions, 175 deletions
diff --git a/src/lib/netlist/solver/mat_cr.h b/src/lib/netlist/solver/mat_cr.h
deleted file mode 100644
index 8693c49c3e1..00000000000
--- a/src/lib/netlist/solver/mat_cr.h
+++ /dev/null
@@ -1,175 +0,0 @@
-// 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"
-#include "../plib/palloc.h"
-
-template<std::size_t N, typename C = uint16_t, typename T = double>
-struct mat_cr_t
-{
- typedef C index_type;
- typedef T value_type;
-
- 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;
-
- explicit 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
- */
-
- 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 std::size_t e = ia[i+1];
- for (; k < e; k++)
- tmp += A[k] * x[ja[k]];
- res[i++] = tmp;
- }
- }
-
- void incomplete_LU_factorization(T * RESTRICT LU)
- {
- /*
- * incomplete LU Factorization according to http://de.wikipedia.org/wiki/ILU-Zerlegung
- *
- * Result is stored in matrix LU
- *
- */
-
- const std::size_t lnz = nz_num;
-
- for (std::size_t k = 0; k < lnz; k++)
- LU[k] = A[k];
-
- for (std::size_t i = 1; ia[i] < lnz; i++) // 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 std::size_t k = ja[pk];
- const std::size_t iak1 = ia[k + 1];
- const T LUpk = LU[pk] = LU[pk] / LU[diag[k]];
-
- std::size_t pt = ia[k];
-
- for (std::size_t pj = pk + 1; pj < iai1; pj++) // pj = (i, j)
- {
- // we can assume that within a row ja increases continuously */
- const std::size_t 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 T * RESTRICT LU, T * 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.
- *
- */
-
- for (std::size_t i = 1; ia[i] < nz_num; ++i )
- {
- T tmp = 0.0;
- const std::size_t j1 = ia[i];
- const std::size_t j2 = diag[i];
-
- 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 (std::size_t i = size; i-- > 0; )
- {
- T tmp = 0.0;
- 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[i] = (r[i] - tmp) / LU[di];
- }
- }
-};
-
-#endif /* MAT_CR_H_ */