summaryrefslogtreecommitdiffstats
path: root/src/lib/netlist/plib/mat_cr.h
diff options
context:
space:
mode:
Diffstat (limited to 'src/lib/netlist/plib/mat_cr.h')
-rw-r--r--src/lib/netlist/plib/mat_cr.h121
1 files changed, 61 insertions, 60 deletions
diff --git a/src/lib/netlist/plib/mat_cr.h b/src/lib/netlist/plib/mat_cr.h
index d87faef8225..9d885ffb884 100644
--- a/src/lib/netlist/plib/mat_cr.h
+++ b/src/lib/netlist/plib/mat_cr.h
@@ -1,15 +1,15 @@
// license:GPL-2.0+
// copyright-holders:Couriersud
-/*
- * mat_cr.h
- *
- * Compressed row format matrices
- *
- */
#ifndef MAT_CR_H_
#define MAT_CR_H_
+///
+/// \file mat_cr.h
+///
+/// Compressed row format matrices
+///
+
#include "palloc.h"
#include "parray.h"
#include "pconfig.h"
@@ -132,7 +132,7 @@ namespace plib
row_idx[size()] = nz;
nz_num = nz;
- /* build nzbd */
+ // build nzbd
for (std::size_t k=0; k < size(); k++)
{
@@ -147,9 +147,8 @@ namespace plib
template <typename VTV, typename VTR>
void mult_vec(VTR & res, const VTV & x) const noexcept
{
- /*
- * res = A * x
- */
+
+ // res = A * x
#if 0
//plib::omp::set_num_threads(4);
plib::omp::for_static(0, constants<std::size_t>::zero(), m_size, [this, &res, &x](std::size_t row)
@@ -176,7 +175,7 @@ namespace plib
#endif
}
- /* throws error if P(source)>P(destination) */
+ // throws error if P(source)>P(destination)
template <typename LUMAT>
void slim_copy_from(LUMAT & src)
{
@@ -185,20 +184,20 @@ namespace plib
C dp = row_idx[r];
for (C sp = src.row_idx[r]; sp < src.row_idx[r+1]; sp++)
{
- /* advance dp to source column and fill 0s if necessary */
+ // advance dp to source column and fill 0s if necessary
while (col_idx[dp] < src.col_idx[sp])
A[dp++] = 0;
if (row_idx[r+1] <= dp || col_idx[dp] != src.col_idx[sp])
pthrow<pexception>("slim_copy_from error");
A[dp++] = src.A[sp];
}
- /* fill remaining elements in row */
+ // fill remaining elements in row
while (dp < row_idx[r+1])
A[dp++] = 0;
}
}
- /* only copies common elements */
+ // only copies common elements
template <typename LUMAT>
void reduction_copy_from(LUMAT & src) noexcept
{
@@ -208,7 +207,7 @@ namespace plib
C dp(row_idx[r]);
while(sp < src.row_idx[r+1])
{
- /* advance dp to source column and fill 0s if necessary */
+ // advance dp to source column and fill 0s if necessary
if (col_idx[dp] < src.col_idx[sp])
A[dp++] = 0;
else if (col_idx[dp] == src.col_idx[sp])
@@ -216,13 +215,13 @@ namespace plib
else
sp++;
}
- /* fill remaining elements in row */
+ // fill remaining elements in row
while (dp < row_idx[r+1])
A[dp++] = 0;
}
}
- /* no checks at all - may crash */
+ // no checks at all - may crash
template <typename LUMAT>
void raw_copy_from(LUMAT & src) noexcept
{
@@ -379,7 +378,7 @@ namespace plib
void gaussian_back_substitution(V1 &V, const V2 &RHS)
{
const std::size_t iN = base::size();
- /* row n-1 */
+ // row n-1
V[iN - 1] = RHS[iN - 1] / base::A[base::diag[iN - 1]];
for (std::size_t j = iN - 1; j-- > 0;)
@@ -397,7 +396,7 @@ namespace plib
void gaussian_back_substitution(V1 &V)
{
const std::size_t iN = base::size();
- /* row n-1 */
+ // row n-1
V[iN - 1] = V[iN - 1] / base::A[base::diag[iN - 1]];
for (std::size_t j = iN - 1; j-- > 0;)
@@ -492,7 +491,7 @@ namespace plib
void build(M &fill, std::size_t ilup)
{
std::size_t p(0);
- /* build ilu_rows */
+ // build ilu_rows
for (decltype(fill.size()) i=1; i < fill.size(); i++)
{
bool found(false);
@@ -526,25 +525,25 @@ namespace plib
m_ILUp = ilup;
}
+ /// \brief incomplete LU Factorization.
+ ///
+ /// We are following http://de.wikipedia.org/wiki/ILU-Zerlegung here.
+ ///
+ /// The result is stored in matrix LU
+ ///
+ /// For i = 1,...,N-1
+ /// For k = 0, ... , i - 1
+ /// If a[i,k] != 0
+ /// a[i,k] = a[i,k] / a[k,k]
+ /// For j = k + 1, ... , N - 1
+ /// If a[i,j] != 0
+ /// a[i,j] = a[i,j] - a[i,k] * a[k,j]
+ /// j=j+1
+ /// k=k+1
+ /// i=i+1
+ ///
void incomplete_LU_factorization(const base &mat)
{
- /*
- * incomplete LU Factorization according to http://de.wikipedia.org/wiki/ILU-Zerlegung
- *
- * Result is stored in matrix LU
- *
- * For i = 1,...,N-1
- * For k = 0, ... , i - 1
- * If a[i,k] != 0
- * a[i,k] = a[i,k] / a[k,k]
- * For j = k + 1, ... , N - 1
- * If a[i,j] != 0
- * a[i,j] = a[i,j] - a[i,k] * a[k,j]
- * j=j+1
- * k=k+1
- * i=i+1
- *
- */
if (m_ILUp < 1)
this->raw_copy_from(mat);
else
@@ -568,7 +567,7 @@ namespace plib
while (i_j < p_i_end && k_j < p_k_end ) // pj = (i, j)
{
- // we can assume that within a row ja increases continuously */
+ // we can assume that within a row ja increases continuously
const std::size_t c_i_j(base::col_idx[i_j]); // row i, column j
const auto c_k_j(base::col_idx[k_j]); // row k, column j
@@ -582,30 +581,32 @@ namespace plib
}
}
+
+ /// \brief Solve a linear equation
+ ///
+ /// 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.
+ ///
template <typename R>
void solveLU (R &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; i < base::size(); ++i )
{
typename base::value_type tmp(0);
@@ -634,4 +635,4 @@ namespace plib
} // namespace plib
-#endif /* MAT_CR_H_ */
+#endif // MAT_CR_H_