2#include <cusp/dia_matrix.h>
3#include <cusp/coo_matrix.h>
66template<
class ContainerType >
71 using HCooMatrix = cusp::coo_matrix<int, value_type, cusp::host_memory>;
72 using HDiaMatrix = cusp::dia_matrix<int, value_type, cusp::host_memory>;
84 m_v = m_vp = m_vm = copyable;
85 m_max_iter = max_iterations;
86 m_iter = max_iterations;
88 set_iter( max_iterations);
91 template<
class ...Params>
101 m_max_iter = new_max;
144 template <
class MatrixType,
class ContainerType0,
class ContainerType1,
145 class ContainerType2,
class FuncTe1>
146 unsigned solve(ContainerType0& x, FuncTe1 f,
147 MatrixType&& A,
const ContainerType1& b,
150 std::string error_norm =
"universal",
155 nrmb_correction, error_norm, res_fac, q);
156 if(
"residual" == error_norm)
159 normMbVy(std::forward<MatrixType>(A), m_TH, m_yH,
x, b,
183 template<
class MatrixType,
class ContainerType0,
class ContainerType1>
187 std::string error_norm =
"universal",
value_type res_fac = 1.,
193 nrmb_correction, error_norm, res_fac, q);
215 template<
class MatrixType,
class DiaMatrixType,
class ContainerType0,
216 class ContainerType1,
class ContainerType2>
217 void normMbVy( MatrixType&& A,
218 const DiaMatrixType& T,
219 const ContainerType0& y,
231 unsigned less_iter = 0;
232 for(
unsigned i=0; i<
y.size(); i++)
236 for (
unsigned i=0; i<less_iter-1; i++)
240 -T.values(i,0)/T.values(i,2), m_vm,
241 -T.values(i,1)/T.values(i,2), m_v,
242 1.0/T.values(i,2), m_vp);
249 template <
class MatrixType,
class ContainerType1,
250 class ContainerType2,
class UnaryOp>
252 MatrixType&& A,
const ContainerType1& b,
255 std::string error_norm =
"residual",
261 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
266 DG_RANK0 std::cout <<
"# Norm of b "<<m_bnorm <<
"\n";
267 DG_RANK0 std::cout <<
"# Res factor "<<res_fac <<
"\n";
268 DG_RANK0 std::cout <<
"# Residual errors: \n";
279 for(
unsigned i=0; i<m_max_iter; i++)
281 m_TH.values(i,0) = betaip;
285 m_TH.values(i,1) = alphai;
291 DG_RANK0 std::cout <<
"beta["<<i+1 <<
"]=0 encountered\n";
295 m_TH.values(i,2) = betaip;
298 if(
"residual" == error_norm)
300 residual = compute_residual_error( m_TH, i)*m_bnorm;
305 if( i>=q &&( (i<=10) || (i>10 && i%10 == 0) ))
307 residual = compute_universal_error( m_TH, i, q, f,
318 DG_RANK0 std::cout <<
"# ||r||_W = " << residual <<
"\tat i = " << i <<
"\n";
319 if (res_fac*residual< eps*(xnorm + nrmb_correction) )
327 set_iter( m_max_iter);
333 return TH.values(iter,2)*fabs(T1);
335 template<
class UnaryOp>
337 unsigned q, UnaryOp f,
HVec& yH)
339 unsigned new_iter = iter + 1 + q;
341 HDiaMatrix THtilde( new_iter, new_iter, 3*new_iter-2, 3);
342 THtilde.diagonal_offsets[0] = -1;
343 THtilde.diagonal_offsets[1] = 0;
344 THtilde.diagonal_offsets[2] = 1;
345 for(
unsigned u=0; u<iter+1; u++)
347 THtilde.values(u,0) = TH.values(u,0);
348 THtilde.values(u,1) = TH.values(u,1);
349 THtilde.values(u,2) = TH.values(u,2);
351 for(
unsigned u=1; u<=q; u++)
353 THtilde.values( iter+u, 0) = u==1 ? TH.values(iter,2) :
354 TH.values( iter+1-u, 1);
355 THtilde.values( iter+u, 1) = TH.values( iter-u, 1);
356 THtilde.values( iter+u, 2) = TH.values( iter-u, 0);
359 HVec yHtilde = f( THtilde);
360 for(
unsigned u=0; u<yH.size(); u++)
368 void set_iter(
unsigned new_iter) {
371 m_TH.resize(new_iter, new_iter, 3*new_iter-2, 3, m_max_iter);
372 m_TH.diagonal_offsets[0] = -1;
373 m_TH.diagonal_offsets[1] = 0;
374 m_TH.diagonal_offsets[2] = 1;
377 ContainerType m_v, m_vp, m_vm;
380 unsigned m_iter, m_max_iter;
381 bool m_verbose =
false;
Tridiagonalize and approximate via Lanczos algorithm. A is self-adjoint in the weights .
Definition: lanczos.h:68
cusp::dia_matrix< int, value_type, cusp::host_memory > HDiaMatrix
Definition: lanczos.h:72
value_type get_bnorm() const
Norm of b from last call to operator()
Definition: lanczos.h:115
get_value_type< ContainerType > value_type
value type of the ContainerType class
Definition: lanczos.h:70
cusp::coo_matrix< int, value_type, cusp::host_memory > HCooMatrix
Definition: lanczos.h:71
UniversalLanczos(const ContainerType ©able, unsigned max_iterations)
Allocate memory for the method.
Definition: lanczos.h:82
void set_verbose(bool verbose)
Set or unset debugging output during iterations.
Definition: lanczos.h:111
dg::HVec HVec
Definition: lanczos.h:73
unsigned get_max() const
Get the current maximum number of iterations.
Definition: lanczos.h:107
unsigned get_iter() const
Get the number of iterations in the last call to tridiag or solve (same as T.num_rows)
Definition: lanczos.h:200
UniversalLanczos()
Allocate nothing, Call construct method before usage.
Definition: lanczos.h:75
unsigned solve(ContainerType0 &x, FuncTe1 f, MatrixType &&A, const ContainerType1 &b, const ContainerType2 &weights, value_type eps, value_type nrmb_correction=1., std::string error_norm="universal", value_type res_fac=1., unsigned q=1)
via Lanczos and matrix function computation. A is self-adjoint in the weights .
Definition: lanczos.h:146
void set_max(unsigned new_max)
Set the maximum number of iterations.
Definition: lanczos.h:100
const HDiaMatrix & tridiag(MatrixType &&A, const ContainerType0 &b, const ContainerType1 &weights, value_type eps=1e-12, value_type nrmb_correction=1., std::string error_norm="universal", value_type res_fac=1., unsigned q=1)
Tridiagonalization of A using Lanczos method.
Definition: lanczos.h:184
void construct(Params &&...ps)
Perfect forward parameters to one of the constructors.
Definition: lanczos.h:92
void copy(const ContainerTypeIn &source, ContainerTypeOut &target)
void axpbypgz(get_value_type< ContainerType > alpha, const ContainerType1 &x, get_value_type< ContainerType > beta, const ContainerType2 &y, get_value_type< ContainerType > gamma, ContainerType &z)
void axpby(get_value_type< ContainerType > alpha, const ContainerType1 &x, get_value_type< ContainerType > beta, ContainerType &y)
void scal(ContainerType &x, get_value_type< ContainerType > alpha)
void symv(MatrixType &&M, const ContainerType1 &x, ContainerType2 &y)
get_value_type< MatrixType > dot(const ContainerType1 &x, const MatrixType &m, const ContainerType2 &y)
MPI_Vector< thrust::host_vector< real_type > > weights(const aRealMPITopology2d< real_type > &g)
value_type compute_Tinv_m1(const cusp::dia_matrix< int, value_type, cusp::host_memory > &T, unsigned size)
Computes the value of via a Thomas algorithm.
Definition: tridiaginv.h:32
typename TensorTraits< std::decay_t< Vector > >::value_type get_value_type
static auto make_Linear_Te1(int exp)
Create a functor that computes directly.
Definition: matrixfunction.h:190
thrust::host_vector< double > HVec
Classes for Krylov space approximations of a Matrix-Vector product.