monolish  0.14.0
MONOlithic LIner equation Solvers for Highly-parallel architecture
cg.cpp
Go to the documentation of this file.
1 #include "../../include/monolish_blas.hpp"
2 #include "../../include/monolish_equation.hpp"
3 #include "../../include/monolish_vml.hpp"
4 #include "../internal/monolish_internal.hpp"
5 
6 #include <fstream>
7 #include <iomanip>
8 #include <string>
9 
10 namespace monolish {
11 
12 template <typename MATRIX, typename T>
14  vector<T> &b) {
15  Logger &logger = Logger::get_instance();
16  logger.solver_in(monolish_func);
17 
18  if (A.get_row() != A.get_col()) {
19  throw std::runtime_error("error, A.row != A.col");
20  }
21  if (A.get_device_mem_stat() != x.get_device_mem_stat() &&
22  A.get_device_mem_stat() != b.get_device_mem_stat()) {
23  throw std::runtime_error("error, A.get_device_mem_stat != "
24  "x.get_device_mem_stat != b.get_device_mem_stat");
25  }
26 
27  vector<T> r(A.get_row(), 0.0);
28  vector<T> p(A.get_row(), 0.0);
29  vector<T> q(A.get_row(), 0.0);
30  vector<T> z(A.get_row(), 0.0);
31 
32  if (A.get_device_mem_stat() == true) {
33  monolish::util::send(r, p, q, z);
34  }
35 
36  this->precond.create_precond(A);
37 
38  // r = b-Ax
39  blas::matvec(A, x, q);
40  vml::sub(b, q, r);
41 
42  // p0 = Mr0
43  this->precond.apply_precond(r, z);
44  blas::copy(z, p);
45 
46  for (size_t iter = 0; iter < this->get_maxiter(); iter++) {
47  blas::matvec(A, p, q);
48 
49  auto tmp = blas::dot(z, r);
50  auto alpha = tmp / blas::dot(p, q);
51 
52  blas::axpy(alpha, p, x);
53 
54  blas::axpy(-alpha, q, r);
55 
56  this->precond.apply_precond(r, z);
57  auto beta = blas::dot(z, r) / tmp;
58 
59  blas::xpay(beta, z, p); // p = z + beta*p
60 
61  T resid = this->get_residual(r);
62  if (this->get_print_rhistory() == true) {
63  *this->rhistory_stream << iter + 1 << "\t" << std::scientific << resid
64  << std::endl;
65  }
66 
67  if (resid < this->get_tol() && this->get_miniter() <= iter + 1) {
68  logger.solver_out();
70  }
71 
72  if (std::isnan(resid)) {
74  }
75  }
76 
77  logger.solver_out();
79 }
80 template int equation::CG<matrix::Dense<double>, double>::monolish_CG(
82 template int equation::CG<matrix::Dense<float>, float>::monolish_CG(
84 
85 template int equation::CG<matrix::CRS<double>, double>::monolish_CG(
87 template int equation::CG<matrix::CRS<float>, float>::monolish_CG(
89 
90 template int equation::CG<matrix::LinearOperator<double>, double>::monolish_CG(
92 template int equation::CG<matrix::LinearOperator<float>, float>::monolish_CG(
94 
96 
97 template <typename MATRIX, typename T>
99  Logger &logger = Logger::get_instance();
100  logger.solver_in(monolish_func);
101 
102  int ret = 0;
103  if (this->get_lib() == 0) {
104  ret = monolish_CG(A, x, b);
105  }
106 
107  logger.solver_out();
108  return ret; // err code
109 }
110 template int equation::CG<matrix::Dense<double>, double>::solve(
112 template int equation::CG<matrix::Dense<float>, float>::solve(
114 
115 template int equation::CG<matrix::CRS<double>, double>::solve(
117 template int equation::CG<matrix::CRS<float>, float>::solve(
119 
120 template int equation::CG<matrix::LinearOperator<double>, double>::solve(
122 template int equation::CG<matrix::LinearOperator<float>, float>::solve(
124 } // namespace monolish
monolish::matrix::LinearOperator
Linear Operator imitating Matrix.
Definition: monolish_coo.hpp:30
monolish_func
#define monolish_func
Definition: monolish_logger.hpp:9
monolish::equation::CG::solve
int solve(MATRIX &A, vector< Float > &x, vector< Float > &b)
solve Ax = b by BiCGSTAB method(lib=0: monolish)
Definition: cg.cpp:98
monolish::Logger
logger class (singleton, for developper class)
Definition: monolish_logger.hpp:19
monolish::blas::axpy
void axpy(const double alpha, const vector< double > &x, vector< double > &y)
axpy: y = ax + y
Definition: vector_blas.cpp:595
monolish::equation::CG::monolish_CG
int monolish_CG(MATRIX &A, vector< Float > &x, vector< Float > &b)
Definition: cg.cpp:13
monolish::Logger::solver_in
void solver_in(const std::string func_name)
Definition: logger_utils.cpp:7
monolish::vml::sub
void sub(const matrix::CRS< double > &A, const matrix::CRS< double > &B, matrix::CRS< double > &C)
element by element subtract CRS matrix A and CRS matrix B.
Definition: matrix_vml.cpp:228
monolish::matrix::Dense
Dense format Matrix.
Definition: monolish_coo.hpp:28
MONOLISH_SOLVER_SUCCESS
#define MONOLISH_SOLVER_SUCCESS
Definition: monolish_common.hpp:10
MONOLISH_SOLVER_MAXITER
#define MONOLISH_SOLVER_MAXITER
Definition: monolish_common.hpp:12
monolish::equation::CG
CG solver class.
Definition: monolish_equation.hpp:63
monolish::blas::copy
void copy(const matrix::Dense< double > &A, matrix::Dense< double > &C)
Dense matrix copy (C=A)
Definition: dense_copy.cpp:25
monolish
Definition: monolish_matrix_blas.hpp:9
monolish::vector::get_device_mem_stat
bool get_device_mem_stat() const
true: sended, false: not send
Definition: monolish_vector.hpp:204
monolish::blas::matvec
void matvec(const matrix::Dense< double > &A, const vector< double > &x, vector< double > &y)
Dense matrix and vector multiplication: y = Ax.
Definition: matvec_blas.cpp:11
monolish::blas::dot
void dot(const vector< double > &x, const vector< double > &y, double &ans)
inner product (dot)
Definition: vector_blas.cpp:974
monolish::blas::xpay
void xpay(const double alpha, const vector< double > &x, vector< double > &y)
xpay: y = x + ay
Definition: vector_blas.cpp:1108
monolish::Logger::solver_out
void solver_out()
Definition: logger_utils.cpp:35
monolish::vector
vector class
Definition: monolish_coo.hpp:25
MONOLISH_SOLVER_RESIDUAL_NAN
#define MONOLISH_SOLVER_RESIDUAL_NAN
Definition: monolish_common.hpp:14
monolish::util::send
auto send(T &x)
send data to GPU
Definition: monolish_common.hpp:598
monolish::Logger::get_instance
static Logger & get_instance()
Definition: monolish_logger.hpp:42
monolish::matrix::CRS
Compressed Row Storage (CRS) format Matrix.
Definition: monolish_coo.hpp:29