monolish  0.14.2
MONOlithic LIner equation Solvers for Highly-parallel architecture
jacobi.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 namespace monolish {
7 
9 // precondition
11 template <typename MATRIX, typename T>
12 void equation::Jacobi<MATRIX, T>::create_precond(MATRIX &A) {
13  Logger &logger = Logger::get_instance();
14  logger.solver_in(monolish_func);
15 
16  if (A.get_row() != A.get_col()) {
17  throw std::runtime_error("error A.row != A.col");
18  }
19 
20  this->precond.M.resize(A.get_row());
21  // send M
22  if (A.get_device_mem_stat() == true) {
23  this->precond.M.send();
24  }
25 
26  A.diag(this->precond.M);
27  vml::reciprocal(this->precond.M, this->precond.M);
28 
29  logger.solver_out();
30 }
31 template void equation::Jacobi<matrix::Dense<float>, float>::create_precond(
32  matrix::Dense<float> &A);
33 template void equation::Jacobi<matrix::Dense<double>, double>::create_precond(
34  matrix::Dense<double> &A);
35 
36 template void equation::Jacobi<matrix::CRS<float>, float>::create_precond(
37  matrix::CRS<float> &A);
38 template void equation::Jacobi<matrix::CRS<double>, double>::create_precond(
39  matrix::CRS<double> &A);
40 
41 template void
42 equation::Jacobi<matrix::LinearOperator<float>, float>::create_precond(
43  matrix::LinearOperator<float> &A);
44 template void
45 equation::Jacobi<matrix::LinearOperator<double>, double>::create_precond(
46  matrix::LinearOperator<double> &A);
47 
48 template <typename MATRIX, typename T>
49 void equation::Jacobi<MATRIX, T>::apply_precond(const vector<T> &r,
50  vector<T> &z) {
51  Logger &logger = Logger::get_instance();
52  logger.solver_in(monolish_func);
53 
54  vml::mul(this->precond.M, r, z); // x = Ab
55 
56  logger.solver_out();
57 }
58 template void equation::Jacobi<matrix::Dense<float>, float>::apply_precond(
59  const vector<float> &r, vector<float> &z);
60 template void equation::Jacobi<matrix::Dense<double>, double>::apply_precond(
61  const vector<double> &r, vector<double> &z);
62 
63 template void equation::Jacobi<matrix::CRS<float>, float>::apply_precond(
64  const vector<float> &r, vector<float> &z);
65 template void equation::Jacobi<matrix::CRS<double>, double>::apply_precond(
66  const vector<double> &r, vector<double> &z);
67 
68 template void
69 equation::Jacobi<matrix::LinearOperator<float>, float>::apply_precond(
70  const vector<float> &r, vector<float> &z);
71 template void
72 equation::Jacobi<matrix::LinearOperator<double>, double>::apply_precond(
73  const vector<double> &r, vector<double> &z);
74 
76 // solver
78 
79 template <typename MATRIX, typename T>
80 int equation::Jacobi<MATRIX, T>::monolish_Jacobi(MATRIX &A, vector<T> &x,
81  vector<T> &b) {
82  Logger &logger = Logger::get_instance();
83  logger.solver_in(monolish_func);
84 
85  if (A.get_row() != A.get_col()) {
86  throw std::runtime_error("error, A.row != A.col");
87  }
88  if (A.get_device_mem_stat() != x.get_device_mem_stat() &&
89  A.get_device_mem_stat() != b.get_device_mem_stat()) {
90  throw std::runtime_error("error, A.get_device_mem_stat != "
91  "x.get_device_mem_stat != b.get_device_mem_stat");
92  }
93 
94  vector<T> r(A.get_row(), 0.0);
95  vector<T> t(A.get_row(), 0.0);
96  vector<T> s(A.get_row(), 0.0);
97  vector<T> d(A.get_row(), 0.0);
98  util::send(r, t, s, d);
99 
100  auto bnrm2 = blas::nrm2(b);
101  bnrm2 = 1.0 / bnrm2;
102  T nrm2 = 0.0;
103 
104  A.diag(d);
105  vml::reciprocal(d, d); // d[i] = 1/d[i]
106 
107  this->precond.create_precond(A);
108 
109  for (size_t iter = 0; iter < this->get_maxiter(); iter++) {
110 
111  /* x += D^{-1}(b - Ax) */
112  this->precond.apply_precond(x, s);
113  blas::copy(x, s);
114  blas::matvec(A, s, t);
115  blas::axpyz(-1, t, b, r);
116  nrm2 = blas::nrm2(r);
117  vml::mul(r, d, r);
118  vml::add(x, r, x);
119 
120  nrm2 = nrm2 * bnrm2;
121 
122  if (this->get_print_rhistory() == true) {
123  *this->rhistory_stream << iter + 1 << "\t" << std::scientific << nrm2
124  << std::endl;
125  }
126 
127  if (nrm2 < this->get_tol() && this->get_miniter() <= iter + 1) {
128  logger.solver_out();
129  return MONOLISH_SOLVER_SUCCESS;
130  }
131 
132  if (std::isnan(nrm2)) {
133  return MONOLISH_SOLVER_RESIDUAL_NAN;
134  }
135  }
136 
137  logger.solver_out();
138  return MONOLISH_SOLVER_NOT_IMPL;
139 }
140 template int equation::Jacobi<matrix::Dense<double>, double>::monolish_Jacobi(
141  matrix::Dense<double> &A, vector<double> &x, vector<double> &b);
142 template int equation::Jacobi<matrix::Dense<float>, float>::monolish_Jacobi(
143  matrix::Dense<float> &A, vector<float> &x, vector<float> &b);
144 
145 template int equation::Jacobi<matrix::CRS<double>, double>::monolish_Jacobi(
146  matrix::CRS<double> &A, vector<double> &x, vector<double> &b);
147 template int equation::Jacobi<matrix::CRS<float>, float>::monolish_Jacobi(
148  matrix::CRS<float> &A, vector<float> &x, vector<float> &b);
149 
150 template int
151 equation::Jacobi<matrix::LinearOperator<double>, double>::monolish_Jacobi(
152  matrix::LinearOperator<double> &A, vector<double> &x, vector<double> &b);
153 template int
154 equation::Jacobi<matrix::LinearOperator<float>, float>::monolish_Jacobi(
155  matrix::LinearOperator<float> &A, vector<float> &x, vector<float> &b);
156 
157 template <typename MATRIX, typename T>
158 int equation::Jacobi<MATRIX, T>::solve(MATRIX &A, vector<T> &x, vector<T> &b) {
159  Logger &logger = Logger::get_instance();
160  logger.solver_in(monolish_func);
161 
162  int ret = 0;
163  if (this->get_lib() == 0) {
164  ret = monolish_Jacobi(A, x, b);
165  }
166 
167  logger.solver_out();
168  return ret; // err code
169 }
170 
171 template int equation::Jacobi<matrix::Dense<float>, float>::solve(
172  matrix::Dense<float> &A, vector<float> &x, vector<float> &b);
173 template int equation::Jacobi<matrix::Dense<double>, double>::solve(
174  matrix::Dense<double> &A, vector<double> &x, vector<double> &b);
175 
176 template int equation::Jacobi<matrix::CRS<float>, float>::solve(
177  matrix::CRS<float> &A, vector<float> &x, vector<float> &b);
178 template int equation::Jacobi<matrix::CRS<double>, double>::solve(
179  matrix::CRS<double> &A, vector<double> &x, vector<double> &b);
180 
181 template int equation::Jacobi<matrix::LinearOperator<float>, float>::solve(
182  matrix::LinearOperator<float> &A, vector<float> &x, vector<float> &b);
183 template int equation::Jacobi<matrix::LinearOperator<double>, double>::solve(
184  matrix::LinearOperator<double> &A, vector<double> &x, vector<double> &b);
185 } // namespace monolish
monolish::matrix::LinearOperator
Linear Operator imitating Matrix.
Definition: monolish_coo.hpp:37
monolish_func
#define monolish_func
Definition: monolish_logger.hpp:9
monolish::equation::Jacobi::solve
int solve(MATRIX &A, vector< Float > &x, vector< Float > &b)
solve Ax = b by jacobi method(lib=0: monolish)
Definition: jacobi.cpp:158
monolish::blas::nrm2
void nrm2(const vector< double > &x, double &ans)
nrm2: ||x||_2
Definition: vector_blas.cpp:1079
MONOLISH_SOLVER_NOT_IMPL
#define MONOLISH_SOLVER_NOT_IMPL
Definition: monolish_common.hpp:18
monolish::Logger
logger class (singleton, for developper class)
Definition: monolish_logger.hpp:19
monolish::Logger::solver_in
void solver_in(const std::string func_name)
Definition: logger_utils.cpp:7
monolish::matrix::Dense
Dense format Matrix.
Definition: monolish_coo.hpp:35
MONOLISH_SOLVER_SUCCESS
#define MONOLISH_SOLVER_SUCCESS
Definition: monolish_common.hpp:13
monolish::vml::reciprocal
void reciprocal(const matrix::CRS< double > &a, matrix::CRS< double > &y)
reciprocal to CRS matrix elements (C[0:nnz] = 1 / A[0:nnz])
Definition: matrix_vml.cpp:375
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:10
monolish::equation::Jacobi::create_precond
void create_precond(MATRIX &A)
Definition: jacobi.cpp:12
monolish::vector::get_device_mem_stat
bool get_device_mem_stat() const
true: sended, false: not send
Definition: monolish_vector.hpp:210
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:10
monolish::equation::Jacobi::monolish_Jacobi
int monolish_Jacobi(MATRIX &A, vector< Float > &x, vector< Float > &b)
Definition: jacobi.cpp:80
monolish::Logger::solver_out
void solver_out()
Definition: logger_utils.cpp:35
monolish::equation::Jacobi::apply_precond
void apply_precond(const vector< Float > &r, vector< Float > &z)
Definition: jacobi.cpp:49
monolish::vml::add
void add(const matrix::CRS< double > &A, const matrix::CRS< double > &B, matrix::CRS< double > &C)
element by element addition CRS matrix A and CRS matrix B.
Definition: matrix_vml.cpp:220
monolish::blas::axpyz
void axpyz(const double alpha, const vector< double > &x, const vector< double > &y, vector< double > &z)
axpyz: z = ax + y
Definition: vector_blas.cpp:665
monolish::vector
vector class
Definition: monolish_coo.hpp:32
monolish::vml::mul
void mul(const matrix::CRS< double > &A, const matrix::CRS< double > &B, matrix::CRS< double > &C)
element by element multiplication CRS matrix A and CRS matrix B.
Definition: matrix_vml.cpp:236
monolish::equation::Jacobi
Jacobi solver class.
Definition: monolish_equation.hpp:159
MONOLISH_SOLVER_RESIDUAL_NAN
#define MONOLISH_SOLVER_RESIDUAL_NAN
Definition: monolish_common.hpp:17
monolish::util::send
auto send(T &x)
send data to GPU
Definition: monolish_common.hpp:642
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:36
Generated by   doxygen 1.8.17