monolish  0.14.2
MONOlithic LIner equation Solvers for Highly-parallel architecture
convert_linearoperator.cpp
Go to the documentation of this file.
1 #include "../../../include/common/monolish_dense.hpp"
2 #include "../../../include/common/monolish_logger.hpp"
3 #include "../../../include/common/monolish_matrix.hpp"
4 #include "../../../include/common/monolish_vector.hpp"
5 #include "../../../include/monolish_blas.hpp"
6 #include "../../internal/monolish_internal.hpp"
7 
8 #include <fstream>
9 #include <iomanip>
10 #include <limits>
11 #include <sstream>
12 
13 namespace monolish {
14 namespace matrix {
15 
16 template <typename T> void LinearOperator<T>::convert(COO<T> &coo) {
17  Logger &logger = Logger::get_instance();
18  logger.util_in(monolish_func);
19 
20  // todo coo err check (only square)
21 
22  rowN = coo.get_row();
23  colN = coo.get_col();
24 
25  gpu_status = coo.get_device_mem_stat();
26 
27  set_matvec([&](const vector<T> &VEC) {
28  CRS<T> crs(coo);
29  vector<T> vec(crs.get_row(), 0);
30  blas::matvec(crs, VEC, vec);
31  return vec;
32  });
33  rmatvec = nullptr;
34  set_matmul_dense([&](const Dense<T> &MAT) {
35  CRS<T> crs(coo);
36  Dense<T> mat(crs.get_row(), MAT.get_col(), 0.0);
37  blas::matmul(crs, MAT, mat);
38  return mat;
39  });
40  rmatmul_dense = nullptr;
41 
42  logger.util_out();
43 }
44 
45 template void LinearOperator<double>::convert(COO<double> &coo);
46 template void LinearOperator<float>::convert(COO<float> &coo);
47 
48 template <typename T> void LinearOperator<T>::convert(CRS<T> &crs) {
49  Logger &logger = Logger::get_instance();
50  logger.util_in(monolish_func);
51 
52  // todo crs err check (only square)
53 
54  rowN = crs.get_row();
55  colN = crs.get_col();
56 
57  gpu_status = crs.get_device_mem_stat();
58 
59  set_matvec([&](const vector<T> &VEC) {
60  vector<T> vec(crs.get_row(), 0);
61  if (gpu_status) {
62  util::send(vec);
63  }
64  blas::matvec(crs, VEC, vec);
65  return vec;
66  });
67  rmatvec = nullptr;
68  set_matmul_dense([&](const Dense<T> &MAT) {
69  Dense<T> mat(crs.get_row(), MAT.get_col(), 0.0);
70  if (gpu_status) {
71  util::send(mat);
72  }
73  blas::matmul(crs, MAT, mat);
74  return mat;
75  });
76  rmatmul_dense = nullptr;
77 
78  logger.util_out();
79 }
80 
81 template void LinearOperator<double>::convert(CRS<double> &crs);
82 template void LinearOperator<float>::convert(CRS<float> &crs);
83 
84 template <typename T> void LinearOperator<T>::convert(Dense<T> &dense) {
85  Logger &logger = Logger::get_instance();
86  logger.util_in(monolish_func);
87 
88  // todo dense err check (only square)
89 
90  rowN = dense.get_row();
91  colN = dense.get_col();
92 
93  gpu_status = dense.get_device_mem_stat();
94 
95  set_matvec([&](const vector<T> &VEC) {
96  vector<T> vec(dense.get_row(), 0);
97  if (gpu_status) {
98  util::send(vec);
99  }
100  blas::matvec(dense, VEC, vec);
101  return vec;
102  });
103  rmatvec = nullptr;
104  set_matmul_dense([&](const Dense<T> &MAT) {
105  Dense<T> mat(dense.get_row(), MAT.get_col(), 0.0);
106  if (gpu_status) {
107  util::send(mat);
108  }
109  blas::matmul(dense, MAT, mat);
110  return mat;
111  });
112  rmatmul_dense = nullptr;
113 
114  logger.util_out();
115 }
116 
117 template void LinearOperator<double>::convert(Dense<double> &crs);
118 template void LinearOperator<float>::convert(Dense<float> &crs);
119 
120 template <typename T>
122  if (!get_matvec_init_flag()) {
123  Dense<T> A(rowN, colN);
124  dense = A;
125  return;
126  }
127 
128  std::vector<T> values(rowN * colN);
129  for (size_t i = 0; i < colN; ++i) {
130  std::vector<T> vec_tmp(colN, 0);
131  vec_tmp[i] = 1;
132  vector<T> vec(vec_tmp);
133  vector<T> ans(rowN);
134  if (gpu_status) {
135  util::send(ans, vec);
136  }
137  ans = matvec(vec);
138  if (gpu_status) {
139  util::recv(ans);
140  }
141  for (size_t j = 0; j < rowN; ++j) {
142  values[j * colN + i] = ans[j];
143  }
144  }
145 
146  dense = Dense<T>(rowN, colN, values);
147 }
148 
151 
152 } // namespace matrix
153 } // namespace monolish
monolish::blas::rmatvec
void rmatvec(const matrix::LinearOperator< double > &A, const vector< double > &x, vector< double > &y)
matrix (LinearOperator) and vector multiplication: y = Ax
Definition: matvec_blas.cpp:250
monolish_func
#define monolish_func
Definition: monolish_logger.hpp:9
monolish::matrix::LinearOperator::convert_to_Dense
void convert_to_Dense(Dense< Float > &dense) const
Definition: convert_linearoperator.cpp:121
monolish::matrix::Dense
Dense format Matrix.
Definition: monolish_coo.hpp:35
monolish
Definition: monolish_matrix_blas.hpp:10
monolish::matrix::LinearOperator::convert
void convert(COO< Float > &coo)
Convert LinearOperator from COO.
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::util::recv
auto recv(T &x)
recv. and free data from GPU
Definition: monolish_common.hpp:656
monolish::vector
vector class
Definition: monolish_coo.hpp:32
monolish::blas::matmul
void matmul(const matrix::Dense< double > &A, const matrix::Dense< double > &B, matrix::Dense< double > &C)
Dense matrix multiplication: C = AB.
Definition: dense-dense_matmul.cpp:7
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