summaryrefslogtreecommitdiffstatshomepage
path: root/src/lib/netlist/plib/gmres.h
diff options
context:
space:
mode:
Diffstat (limited to 'src/lib/netlist/plib/gmres.h')
-rw-r--r--src/lib/netlist/plib/gmres.h499
1 files changed, 287 insertions, 212 deletions
diff --git a/src/lib/netlist/plib/gmres.h b/src/lib/netlist/plib/gmres.h
index 2c357e97624..6ce96ed2ae6 100644
--- a/src/lib/netlist/plib/gmres.h
+++ b/src/lib/netlist/plib/gmres.h
@@ -1,36 +1,50 @@
-// license:GPL-2.0+
+// license:BSD-3-Clause
// copyright-holders:Couriersud
-/*
- * gmres.h
- *
- */
#ifndef PLIB_GMRES_H_
#define PLIB_GMRES_H_
-#include "mat_cr.h"
+// Names
+// spell-checker: words Burkardt, Saad, Yousef, Givens
+//
+
+///
+/// \file gmres.h
+///
+
#include "parray.h"
#include "pconfig.h"
+#include "pmatrix_cr.h"
#include "vector_ops.h"
#include <algorithm>
-#include <cmath>
-
namespace plib
{
- template <typename FT, int SIZE>
+ template <int k>
+ struct do_k_helper
+ {
+ static constexpr bool value = true;
+ };
+
+ template <>
+ struct do_k_helper<-1>
+ {
+ static constexpr float value = 0.0;
+ };
+
+ template <typename ARENA, typename FT, int SIZE>
struct mat_precondition_ILU
{
- using mat_type = plib::matrix_compressed_rows_t<FT, SIZE>;
-
- mat_precondition_ILU(std::size_t size, int ilu_scale = 4
- , std::size_t bw = plib::matrix_compressed_rows_t<FT, SIZE>::FILL_INFINITY)
- : m_mat(static_cast<typename mat_type::index_type>(size))
- , m_LU(static_cast<typename mat_type::index_type>(size))
- , m_use_iLU_preconditioning(ilu_scale >= 0)
- , m_ILU_scale(static_cast<std::size_t>(ilu_scale))
+ using mat_type = plib::pmatrix_cr<ARENA, FT, SIZE>;
+ using matLU_type = plib::pLUmatrix_cr<mat_type>;
+
+ mat_precondition_ILU(ARENA &arena, std::size_t size, std::size_t ilu_scale = 4
+ , std::size_t bw = plib::pmatrix_cr<ARENA, FT, SIZE>::FILL_INFINITY)
+ : m_mat(arena, narrow_cast<typename mat_type::index_type>(size))
+ , m_LU(arena, narrow_cast<typename mat_type::index_type>(size))
+ , m_ILU_scale(narrow_cast<std::size_t>(ilu_scale))
, m_band_width(bw)
{
}
@@ -39,16 +53,11 @@ namespace plib
void build(M &fill)
{
m_mat.build_from_fill_mat(fill, 0);
- if (m_use_iLU_preconditioning)
- {
- m_LU.gaussian_extend_fill_mat(fill);
- m_LU.build_from_fill_mat(fill, m_ILU_scale, m_band_width); // ILU(2)
- //m_LU.build_from_fill_mat(fill, 9999, 20); // Band matrix width 20
- }
+ m_LU.build(fill, m_ILU_scale);
}
- template<typename R, typename V>
+ template <typename R, typename V>
void calc_rhs(R &rhs, const V &v)
{
m_mat.mult_vec(rhs, v);
@@ -56,41 +65,30 @@ namespace plib
void precondition()
{
- if (m_use_iLU_preconditioning)
- {
- if (m_ILU_scale < 1)
- m_LU.raw_copy_from(m_mat);
- else
- m_LU.reduction_copy_from(m_mat);
- m_LU.incomplete_LU_factorization();
- }
+ m_LU.incomplete_LU_factorization(m_mat);
}
- template<typename V>
- void solve_LU_inplace(V &v)
+ template <typename V>
+ void solve_inplace(V &v)
{
- if (m_use_iLU_preconditioning)
- {
- m_LU.solveLUx(v);
- }
+ m_LU.solveLU(v);
}
PALIGNAS_VECTOROPT()
mat_type m_mat;
PALIGNAS_VECTOROPT()
- mat_type m_LU;
- bool m_use_iLU_preconditioning;
+ matLU_type m_LU;
std::size_t m_ILU_scale;
std::size_t m_band_width;
};
- template <typename FT, int SIZE>
- struct mat_precondition_diag
+ template <typename ARENA, typename FT, int SIZE>
+ struct mat_precondition_diagonal
{
- mat_precondition_diag(std::size_t size)
- : m_mat(size)
- , m_diag(size)
- , m_use_iLU_preconditioning(true)
+ mat_precondition_diagonal(ARENA &arena, std::size_t size, [[maybe_unused]] int dummy = 0)
+ : m_mat(arena, size)
+ , m_diagonal(size)
+ , nz_col(size)
{
}
@@ -98,9 +96,21 @@ namespace plib
void build(M &fill)
{
m_mat.build_from_fill_mat(fill, 0);
+ for (std::size_t i = 0; i< m_diagonal.size(); i++)
+ {
+ for (std::size_t j = 0; j< m_diagonal.size(); j++)
+ {
+ std::size_t k=m_mat.row_idx[j];
+ while (m_mat.col_idx[k] < i && k < m_mat.row_idx[j+1])
+ k++;
+ if (m_mat.col_idx[k] == i && k < m_mat.row_idx[j+1])
+ nz_col[i].push_back(k);
+ }
+ nz_col[i].push_back(narrow_cast<std::size_t>(-1));
+ }
}
- template<typename R, typename V>
+ template <typename R, typename V>
void calc_rhs(R &rhs, const V &v)
{
m_mat.mult_vec(rhs, v);
@@ -108,87 +118,148 @@ namespace plib
void precondition()
{
- if (m_use_iLU_preconditioning)
+ for (std::size_t i = 0; i< m_diagonal.size(); i++)
{
- for (std::size_t i = 0; i< m_diag.size(); i++)
+ // ILUT: 265%
+ FT v(0.0);
+#if 0
+ // doesn't works, Mame performance drops significantly%
+ // 136%
+ for (std::size_t j = m_mat.row_idx[i]; j< m_mat.row_idx[i+1]; j++)
+ v += m_mat.A[j] * m_mat.A[j];
+ m_diagonal[i] = reciprocal(std::sqrt(v));
+#elif 0
+ // works halfway, i.e. Mame performance 50%
+ // 147% - lowest average solution time with 7.094
+ for (std::size_t j = m_mat.row_idx[i]; j< m_mat.row_idx[i+1]; j++)
+ v += m_mat.A[j] * m_mat.A[j];
+ m_diagonal[i] = m_mat.A[m_mat.diagonal[i]] / v;
+#elif 0
+ // works halfway, i.e. Mame performance 50%
+ // sum over column i
+ // 344% - lowest average solution time with 3.06
+ std::size_t nzcolp = 0;
+ const auto &nz = nz_col[i];
+ std::size_t j;
+
+ while ((j = nz[nzcolp++])!=narrow_cast<std::size_t>(-1)) // NOLINT(bugprone-infinite-loop)
{
- m_diag[i] = 1.0 / m_mat.A[m_mat.diag[i]];
+ v += m_mat.A[j] * m_mat.A[j];
}
+ m_diagonal[i] = m_mat.A[m_mat.diagonal[i]] / v;
+#elif 0
+ // works halfway, i.e. Mame performance 50%
+ // 151%
+ for (std::size_t j = m_mat.row_idx[i]; j< m_mat.row_idx[i+1]; j++)
+ v += plib::abs(m_mat.A[j]);
+ m_diagonal[i] = reciprocal(v);
+#else
+ // 124%
+ for (std::size_t j = m_mat.row_idx[i]; j< m_mat.row_idx[i+1]; j++)
+ v = std::max(v, plib::abs(m_mat.A[j]));
+ m_diagonal[i] = reciprocal(v);
+#endif
+ //m_diagonal[i] = reciprocal(m_mat.A[m_mat.diagonal[i]]);
}
}
- template<typename V>
- void solve_LU_inplace(V &v)
+ template <typename V>
+ void solve_inplace(V &v)
{
- if (m_use_iLU_preconditioning)
- {
- for (std::size_t i = 0; i< m_diag.size(); i++)
- v[i] = v[i] * m_diag[i];
- }
+ for (std::size_t i = 0; i< m_diagonal.size(); i++)
+ v[i] = v[i] * m_diagonal[i];
}
- plib::matrix_compressed_rows_t<FT, SIZE> m_mat;
- plib::parray<FT, SIZE> m_diag;
- bool m_use_iLU_preconditioning;
+ plib::pmatrix_cr<ARENA, FT, SIZE> m_mat;
+ plib::parray<FT, SIZE> m_diagonal;
+ plib::parray<std::vector<std::size_t>, SIZE > nz_col;
};
- /* FIXME: hardcoding RESTART to 20 becomes an issue on very large
- * systems.
- */
- template <typename FT, int SIZE, int RESTART = 20>
+ template <typename ARENA, typename FT, int SIZE>
+ struct mat_precondition_none
+ {
+ mat_precondition_none(std::size_t size, [[maybe_unused]] int dummy = 0)
+ : m_mat(size)
+ {
+ }
+
+ template <typename M>
+ void build(M &fill)
+ {
+ m_mat.build_from_fill_mat(fill, 0);
+ }
+
+ template <typename R, typename V>
+ void calc_rhs(R &rhs, const V &v)
+ {
+ m_mat.mult_vec(rhs, v);
+ }
+
+ void precondition()
+ {
+ }
+
+ template <typename V>
+ void solve_inplace([[maybe_unused]] V &v)
+ {
+ }
+
+ plib::pmatrix_cr<ARENA, FT, SIZE> m_mat;
+ };
+
+ // FIXME: hard coding RESTART to 20 becomes an issue on very large
+ // systems.
+
+ template <typename FT, int SIZE, int RESTARTMAX = 16>
struct gmres_t
{
public:
using float_type = FT;
- // FIXME: dirty hack to make this compile
- static constexpr const std::size_t storage_N = plib::sizeabs<FT, SIZE>::ABS();
- gmres_t(std::size_t size)
+ //constexpr static int RESTART = RESTARTMAX;
+ constexpr static const int RESTART = (SIZE > 0) ? ((SIZE < RESTARTMAX) ? SIZE : RESTARTMAX)
+ : ((SIZE < 0) ? ((-SIZE < RESTARTMAX) ? -SIZE : RESTARTMAX) : RESTARTMAX);
+
+ explicit gmres_t(std::size_t size)
: residual(size)
, Ax(size)
+ , m_ht(RESTART +1, RESTART)
+ , m_v(RESTART + 1, size)
, m_size(size)
, m_use_more_precise_stop_condition(false)
{
}
- void givens_mult( const FT c, const FT s, FT & g0, FT & g1 )
- {
- const FT g0_last(g0);
-
- g0 = c * g0 - s * g1;
- g1 = s * g0_last + c * g1;
- }
-
- std::size_t size() const { return (SIZE<=0) ? m_size : static_cast<std::size_t>(SIZE); }
+ std::size_t size() const { return (SIZE<=0) ? m_size : narrow_cast<std::size_t>(SIZE); }
template <typename OPS, typename VT, typename VRHS>
std::size_t solve(OPS &ops, VT &x, const VRHS & rhs, const std::size_t itr_max, float_type accuracy)
{
- /*-------------------------------------------------------------------------
- * The code below was inspired by code published by John Burkardt under
- * the LPGL here:
- *
- * http://people.sc.fsu.edu/~jburkardt/cpp_src/mgmres/mgmres.html
- *
- * The code below was completely written from scratch based on the pseudo code
- * found here:
- *
- * http://de.wikipedia.org/wiki/GMRES-Verfahren
- *
- * The Algorithm itself is described in
- *
- * Yousef Saad,
- * Iterative Methods for Sparse Linear Systems,
- * Second Edition,
- * SIAM, 20003,
- * ISBN: 0898715342,
- * LC: QA188.S17.
- *
- *------------------------------------------------------------------------*/
+ // -------------------------------------------------------------------------
+ // The code below was inspired by code published by John Burkardt under
+ // the LPGL here:
+ //
+ // http://people.sc.fsu.edu/~jburkardt/cpp_src/mgmres/mgmres.html
+ //
+ // The code below was completely written from scratch based on the pseudo code
+ // found here:
+ //
+ // http://de.wikipedia.org/wiki/GMRES-Verfahren
+ //
+ // The Algorithm itself is described in
+ //
+ // Yousef Saad,
+ // Iterative Methods for Sparse Linear Systems,
+ // Second Edition,
+ // SIAM, 20003,
+ // ISBN: 0898715342,
+ // LC: QA188.S17.
+ //
+ //------------------------------------------------------------------------
std::size_t itr_used = 0;
- double rho_delta = 0.0;
+ float_type rho_delta(plib::constants<float_type>::zero());
const std::size_t n = size();
@@ -196,165 +267,174 @@ namespace plib
if (m_use_more_precise_stop_condition)
{
- /* derive residual for a given delta x
- *
- * LU y = A dx
- *
- * ==> rho / accuracy = sqrt(y * y)
- *
- * This approach will approximate the iterative stop condition
- * based |xnew - xold| pretty precisely. But it is slow, or expressed
- * differently: The invest doesn't pay off.
- */
-
- vec_set_scalar(n, residual, accuracy);
+ // derive residual for a given delta x
+ //
+ // LU y = A dx
+ //
+ // ==> rho / accuracy = sqrt(y * y)
+ //
+ // This approach will approximate the iterative stop condition
+ // based `|xnew - xold|` pretty precisely. But it is slow, or expressed
+ // differently: The invest doesn't pay off.
+ //
+
+ vec_set_scalar(residual, accuracy);
ops.calc_rhs(Ax, residual);
- ops.solve_LU_inplace(Ax);
+ ops.solve_inplace(Ax);
- const float_type rho_to_accuracy = std::sqrt(vec_mult2<FT>(n, Ax)) / accuracy;
+ const float_type rho_to_accuracy = plib::sqrt(vec_mult2<FT>(Ax)) / accuracy;
rho_delta = accuracy * rho_to_accuracy;
}
else
- rho_delta = accuracy * std::sqrt(static_cast<FT>(n));
-
- /*
- * Using
- *
- * vec_set(n, x, rhs);
- * ops.solve_LU_inplace(x);
- *
- * to get a starting point for x degrades convergence speed compared
- * to using the last solution for x.
- *
- * LU x = b; solve for x;
- *
- */
+ //rho_delta = accuracy * plib::sqrt(vec_mult2<FT>(n, rhs))
+ // + 1e-4 * std::sqrt(n);
+ rho_delta = accuracy * plib::sqrt(narrow_cast<FT>(n));
+
+ //
+ // LU x = b; solve for x;
+ //
+ // Using
+ //
+ // vec_set(n, x, rhs);
+ // ops.solve_inplace(x);
+ //
+ // to get a starting point for x degrades convergence speed compared
+ // to using the last solution for x.
while (itr_used < itr_max)
{
- std::size_t last_k = RESTART;
float_type rho;
ops.calc_rhs(Ax, x);
- vec_sub(n, residual, rhs, Ax);
+ vec_sub(residual, rhs, Ax);
- ops.solve_LU_inplace(residual);
+ ops.solve_inplace(residual);
- rho = std::sqrt(vec_mult2<FT>(n, residual));
+ rho = plib::sqrt(vec_mult2<FT>(residual));
if (rho < rho_delta)
return itr_used + 1;
- /* FIXME: The "+" is necessary to avoid link issues
- * on some systems / compiler versions. Issue reported by
- * AJR, no details known yet.
- */
- vec_set_scalar(RESTART+1, m_g, +constants<FT>::zero());
+ // FIXME: The "+" is necessary to avoid link issues
+ // on some systems / compiler versions. Issue reported by
+ // AJR, no details known yet.
+
+ vec_set_scalar(m_g, +constants<FT>::zero());
m_g[0] = rho;
- //for (std::size_t i = 0; i < mr + 1; i++)
- // vec_set_scalar(mr, m_ht[i], NL_FCONST(0.0));
+ vec_mult_scalar(m_v[0], residual, plib::reciprocal(rho));
- vec_mult_scalar(n, m_v[0], residual, constants<FT>::one() / rho);
+ if (do_k<RESTART-1>(ops, x, itr_used, rho_delta, true))
+ // converged
+ break;
+ }
+ return itr_used;
+ }
- for (std::size_t k = 0; k < RESTART; k++)
- {
- const std::size_t kp1 = k + 1;
+ private:
- ops.calc_rhs(m_v[kp1], m_v[k]);
- ops.solve_LU_inplace(m_v[kp1]);
+ static void givens_mult(FT c, FT s, FT & g0, FT & g1 )
+ {
+ const FT g0_last(g0);
+
+ g0 = c * g0 - s * g1;
+ g1 = s * g0_last + c * g1;
+ }
- for (std::size_t j = 0; j <= k; j++)
- {
- m_ht[j][k] = vec_mult<FT>(n, m_v[kp1], m_v[j]);
- vec_add_mult_scalar(n, m_v[kp1], m_v[j], -m_ht[j][k]);
- }
- m_ht[kp1][k] = std::sqrt(vec_mult2<FT>(n, m_v[kp1]));
+ template <int k, typename OPS, typename VT>
+ bool do_k(OPS &ops, VT &x, std::size_t &itr_used, FT rho_delta, [[maybe_unused]] bool dummy)
+ {
+ if (do_k<k-1, OPS>(ops, x, itr_used, rho_delta, do_k_helper<k-1>::value))
+ return true;
+
+ constexpr const std::size_t kp1 = k + 1;
+ //const std::size_t n = size();
- if (m_ht[kp1][k] != 0.0)
- vec_scale(n, m_v[kp1], constants<FT>::one() / m_ht[kp1][k]);
+ ops.calc_rhs(m_v[kp1], m_v[k]);
+ ops.solve_inplace(m_v[kp1]);
- for (std::size_t j = 0; j < k; j++)
- givens_mult(m_c[j], m_s[j], m_ht[j][k], m_ht[j+1][k]);
+ for (std::size_t j = 0; j <= k; j++)
+ {
+ m_ht[j][k] = vec_mult<FT>(m_v[kp1], m_v[j]);
+ vec_add_mult_scalar(m_v[kp1], m_v[j], -m_ht[j][k]);
+ }
+ m_ht[kp1][k] = plib::sqrt(vec_mult2<FT>(m_v[kp1]));
- const float_type mu = 1.0 / std::hypot(m_ht[k][k], m_ht[kp1][k]);
+ // FIXME: comparison to zero
+ if (m_ht[kp1][k] != plib::constants<FT>::zero())
+ vec_scale(m_v[kp1], reciprocal(m_ht[kp1][k]));
- m_c[k] = m_ht[k][k] * mu;
- m_s[k] = -m_ht[kp1][k] * mu;
- m_ht[k][k] = m_c[k] * m_ht[k][k] - m_s[k] * m_ht[kp1][k];
- m_ht[kp1][k] = 0.0;
+ for (std::size_t j = 0; j < k; j++)
+ givens_mult(m_c[j], m_s[j], m_ht[j][k], m_ht[j+1][k]);
- givens_mult(m_c[k], m_s[k], m_g[k], m_g[kp1]);
+ const float_type mu = reciprocal(plib::hypot(m_ht[k][k], m_ht[kp1][k]));
- rho = std::abs(m_g[kp1]);
+ m_c[k] = m_ht[k][k] * mu;
+ m_s[k] = -m_ht[kp1][k] * mu;
+ m_ht[k][k] = m_c[k] * m_ht[k][k] - m_s[k] * m_ht[kp1][k];
+ m_ht[kp1][k] = plib::constants<FT>::zero();
- itr_used = itr_used + 1;
+ givens_mult(m_c[k], m_s[k], m_g[k], m_g[kp1]);
- if (rho <= rho_delta)
- {
- last_k = k;
- break;
- }
- }
+ const float_type rho = plib::abs(m_g[kp1]);
- if (last_k >= RESTART)
- /* didn't converge within accuracy */
- last_k = RESTART - 1;
+ // FIXME ..
+ itr_used = itr_used + 1;
- /* Solve the system H * y = g */
- /* x += m_v[j] * m_y[j] */
- for (std::size_t i = last_k + 1; i-- > 0;)
+ if (rho <= rho_delta || k == RESTART-1)
+ {
+ // Solve the system H * y = g
+ // x += m_v[j] * m_y[j]
+ for (std::size_t i = k + 1; i-- > 0;)
{
- double tmp = m_g[i];
- for (std::size_t j = i + 1; j <= last_k; j++)
+ auto tmp = m_g[i];
+ const auto ht_i_i = plib::reciprocal(m_ht[i][i]);
+ for (std::size_t j = i + 1; j <= k; j++)
tmp -= m_ht[i][j] * m_y[j];
- m_y[i] = tmp / m_ht[i][i];
+ m_y[i] = tmp * ht_i_i;
+ vec_add_mult_scalar(x, m_v[i], m_y[i]);
}
- for (std::size_t i = 0; i <= last_k; i++)
- vec_add_mult_scalar(n, x, m_v[i], m_y[i]);
-
- if (rho <= rho_delta)
- break;
-
+ //for (std::size_t i = 0; i <= k; i++)
+ // vec_add_mult_scalar(n, x, m_v[i], m_y[i]);
+ return true;
}
- return itr_used;
+ return false;
}
- private:
-
- //typedef typename plib::mat_cr_t<FT, SIZE>::index_type mattype;
+ template <int k, typename OPS, typename VT>
+ bool do_k(OPS &ops, VT &x, std::size_t &itr_used, FT rho_delta, float dummy)
+ {
+ plib::unused_var(ops, x, itr_used, rho_delta, dummy);
+ return false;
+ }
plib::parray<float_type, SIZE> residual;
plib::parray<float_type, SIZE> Ax;
- plib::parray<float_type, RESTART + 1> m_c; /* mr + 1 */
- plib::parray<float_type, RESTART + 1> m_g; /* mr + 1 */
- plib::parray<plib::parray<float_type, RESTART>, RESTART + 1> m_ht; /* (mr + 1), mr */
- plib::parray<float_type, RESTART + 1> m_s; /* mr + 1 */
- plib::parray<float_type, RESTART + 1> m_y; /* mr + 1 */
+ plib::parray<float_type, RESTART + 1> m_c; // mr + 1
+ plib::parray<float_type, RESTART + 1> m_g; // mr + 1
+ plib::parray2D<float_type, RESTART + 1, RESTART> m_ht; // (mr + 1), mr
+ plib::parray<float_type, RESTART + 1> m_s; // mr + 1
+ plib::parray<float_type, RESTART + 1> m_y; // mr + 1
- //plib::parray<float_type, SIZE> m_v[RESTART + 1]; /* mr + 1, n */
- plib::parray<plib::parray<float_type, storage_N>, RESTART + 1> m_v; /* mr + 1, n */
+ plib::parray2D<float_type, RESTART + 1, SIZE> m_v; // mr + 1, n
std::size_t m_size;
bool m_use_more_precise_stop_condition;
-
-
};
#if 0
- /* Example of a Chebyshev iteration solver. This one doesn't work yet,
- * it needs to be extended for non-symmetric matrix operation and
- * depends on spectral radius estimates - which we don't have.
- *
- * Left here as another example.
- */
+ // Example of a Chebyshev iteration solver. This one doesn't work yet,
+ // it needs to be extended for non-symmetric matrix operation and
+ // depends on spectral radius estimates - which we don't have.
+ //
+ // Left here as another example.
template <typename FT, int SIZE>
struct ch_t
@@ -375,16 +455,11 @@ namespace plib
{
}
- std::size_t size() const { return (SIZE<=0) ? m_size : static_cast<std::size_t>(SIZE); }
+ std::size_t size() const { return (SIZE<=0) ? m_size : narrow_cast<std::size_t>(SIZE); }
template <typename OPS, typename VT, typename VRHS>
std::size_t solve(OPS &ops, VT &x0, const VRHS & rhs, const std::size_t iter_max, float_type accuracy)
{
- /*-------------------------------------------------------------------------
- *
- *
- *------------------------------------------------------------------------*/
-
ops.precondition();
const FT lmax = 20.0;
@@ -404,13 +479,13 @@ namespace plib
ops.calc_rhs(Ax, x);
vec_sub(size(), rhs, Ax, residual);
- FT rho_delta = accuracy * std::sqrt(static_cast<FT>(size()));
+ FT rho_delta = accuracy * std::sqrt(narrow_cast<FT>(size()));
rho_delta = 1e-9;
for (int i = 0; i < iter_max; i++)
{
- ops.solve_LU_inplace(residual);
+ ops.solve_inplace(residual);
if (i==0)
{
vec_set(size(), p, residual);
@@ -419,7 +494,7 @@ namespace plib
else
{
beta = alpha * ( c / 2.0)*( c / 2.0);
- alpha = 1.0 / (d - beta);
+ alpha = reciprocal(d - beta);
for (std::size_t k = 0; k < size(); k++)
p[k] = residual[k] + beta * p[k];
}
@@ -435,7 +510,7 @@ namespace plib
}
private:
- //typedef typename plib::mat_cr_t<FT, SIZE>::index_type mattype;
+ //#typedef typename plib::mat_cr_t<FT, SIZE>::index_type mattype;
plib::parray<float_type, SIZE> residual;
plib::parray<float_type, SIZE> Ax;
@@ -447,4 +522,4 @@ namespace plib
} // namespace plib
-#endif /* PLIB_GMRES_H_ */
+#endif // PLIB_GMRES_H_