suanPan
SparseMatCUDA.hpp
Go to the documentation of this file.
1/*******************************************************************************
2 * Copyright (C) 2017-2024 Theodore Chang
3 *
4 * This program is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program. If not, see <http://www.gnu.org/licenses/>.
16 ******************************************************************************/
31// ReSharper disable CppCStyleCast
32#ifndef SPARSEMATCUDA_HPP
33#define SPARSEMATCUDA_HPP
34
35#ifdef SUANPAN_CUDA
36
37#include <cusolverSp.h>
38#include <cusparse.h>
39#include "SparseMat.hpp"
40#include "csr_form.hpp"
41
42template<sp_d T> class SparseMatCUDA final : public SparseMat<T> {
43 cusolverSpHandle_t handle = nullptr;
44 cudaStream_t stream = nullptr;
45 cusparseMatDescr_t descr = nullptr;
46
47 void* d_val_idx = nullptr;
48 void* d_col_idx = nullptr;
49 void* d_row_ptr = nullptr;
50
51 void acquire() {
52 cusolverSpCreate(&handle);
53 cudaStreamCreate(&stream);
54 cusolverSpSetStream(handle, stream);
55 cusparseCreateMatDescr(&descr);
56 cusparseSetMatType(descr, CUSPARSE_MATRIX_TYPE_GENERAL);
57 cusparseSetMatIndexBase(descr, CUSPARSE_INDEX_BASE_ZERO);
58 cudaMalloc(&d_row_ptr, sizeof(int) * (this->n_rows + 1));
59 }
60
61 void release() const {
62 if(handle) cusolverSpDestroy(handle);
63 if(stream) cudaStreamDestroy(stream);
64 if(descr) cusparseDestroyMatDescr(descr);
65 if(d_row_ptr) cudaFree(d_row_ptr);
66 }
67
68 template<sp_d ET> void device_alloc(csr_form<ET, int>&& csr_mat) {
69 const size_t n_val = sizeof(ET) * csr_mat.n_elem;
70 const size_t n_col = sizeof(int) * csr_mat.n_elem;
71
72 cudaMalloc(&d_val_idx, n_val);
73 cudaMalloc(&d_col_idx, n_col);
74
75 cudaMemcpyAsync(d_val_idx, csr_mat.val_mem(), n_val, cudaMemcpyHostToDevice, stream);
76 cudaMemcpyAsync(d_col_idx, csr_mat.col_mem(), n_col, cudaMemcpyHostToDevice, stream);
77 cudaMemcpyAsync(d_row_ptr, csr_mat.row_mem(), sizeof(int) * (csr_mat.n_rows + 1llu), cudaMemcpyHostToDevice, stream);
78 }
79
80 void device_dealloc() const {
81 if(d_val_idx) cudaFree(d_val_idx);
82 if(d_col_idx) cudaFree(d_col_idx);
83 }
84
85protected:
86 using SparseMat<T>::direct_solve;
87
88 int direct_solve(Mat<T>&, const Mat<T>&) override;
89
90public:
91 SparseMatCUDA(const uword in_row, const uword in_col, const uword in_elem = 0)
92 : SparseMat<T>(in_row, in_col, in_elem) { acquire(); }
93
94 SparseMatCUDA(const SparseMatCUDA& other)
95 : SparseMat<T>(other) { acquire(); }
96
97 SparseMatCUDA(SparseMatCUDA&&) noexcept = delete;
98 SparseMatCUDA& operator=(const SparseMatCUDA&) = delete;
99 SparseMatCUDA& operator=(SparseMatCUDA&&) noexcept = delete;
100
101 ~SparseMatCUDA() override {
102 release();
103 device_dealloc();
104 }
105
106 unique_ptr<MetaMat<T>> make_copy() override { return std::make_unique<SparseMatCUDA>(*this); }
107};
108
109template<sp_d T> int SparseMatCUDA<T>::direct_solve(Mat<T>& X, const Mat<T>& B) {
110 if(!this->factored) {
111 // deallocate memory previously allocated for csr matrix
112 device_dealloc();
113
114 std::is_same_v<T, float> || Precision::MIXED == this->setting.precision ? device_alloc(csr_form<float, int>(this->triplet_mat)) : device_alloc(csr_form<double, int>(this->triplet_mat));
115
116 this->factored = true;
117 }
118
119 const size_t n_rhs = (std::is_same_v<T, float> || Precision::MIXED == this->setting.precision ? sizeof(float) : sizeof(double)) * B.n_elem;
120
121 void* d_b = nullptr;
122 void* d_x = nullptr;
123
124 cudaMalloc(&d_b, n_rhs);
125 cudaMalloc(&d_x, n_rhs);
126
127 int singularity;
128 auto code = 0;
129
130 if constexpr(std::is_same_v<T, float>) {
131 cudaMemcpyAsync(d_b, B.memptr(), n_rhs, cudaMemcpyHostToDevice, stream);
132
133 for(auto I = 0llu; I < B.n_elem; I += B.n_rows) code += cusolverSpScsrlsvqr(handle, int(this->n_rows), int(this->triplet_mat.n_elem), descr, (float*)d_val_idx, (int*)d_row_ptr, (int*)d_col_idx, (float*)d_b + I, float(this->setting.tolerance), 3, (float*)d_x + I, &singularity);
134
135 X.set_size(arma::size(B));
136
137 cudaMemcpyAsync(X.memptr(), d_x, n_rhs, cudaMemcpyDeviceToHost, stream);
138
139 cudaDeviceSynchronize();
140 }
141 else if(Precision::FULL == this->setting.precision) {
142 cudaMemcpyAsync(d_b, B.memptr(), n_rhs, cudaMemcpyHostToDevice, stream);
143
144 for(auto I = 0llu; I < B.n_elem; I += B.n_rows) code += cusolverSpDcsrlsvqr(handle, int(this->n_rows), int(this->triplet_mat.n_elem), descr, (double*)d_val_idx, (int*)d_row_ptr, (int*)d_col_idx, (double*)d_b + I, this->setting.tolerance, 3, (double*)d_x + I, &singularity);
145
146 X.set_size(arma::size(B));
147
148 cudaMemcpyAsync(X.memptr(), d_x, n_rhs, cudaMemcpyDeviceToHost, stream);
149
150 cudaDeviceSynchronize();
151 }
152 else {
153 X = arma::zeros(arma::size(B));
154
155 mat full_residual = B;
156
157 auto multiplier = norm(full_residual);
158
159 auto counter = 0u;
160 while(counter++ < this->setting.iterative_refinement) {
161 if(multiplier < this->setting.tolerance) break;
162
163 auto residual = conv_to<fmat>::from(full_residual / multiplier);
164
165 cudaMemcpyAsync(d_b, residual.memptr(), n_rhs, cudaMemcpyHostToDevice, stream);
166
167 code = 0;
168 for(auto I = 0llu; I < B.n_elem; I += B.n_rows) code += cusolverSpScsrlsvqr(handle, int(this->n_rows), int(this->triplet_mat.n_elem), descr, (float*)d_val_idx, (int*)d_row_ptr, (int*)d_col_idx, (float*)d_b + I, float(this->setting.tolerance), 3, (float*)d_x + I, &singularity);
169 if(0 != code) break;
170
171 cudaMemcpyAsync(residual.memptr(), d_x, n_rhs, cudaMemcpyDeviceToHost, stream);
172
173 cudaDeviceSynchronize();
174
175 const mat incre = multiplier * conv_to<mat>::from(residual);
176
177 X += incre;
178
179 suanpan_debug("Mixed precision algorithm multiplier: {:.5E}.\n", multiplier = arma::norm(full_residual -= this->operator*(incre)));
180 }
181 }
182
183 if(d_b) cudaFree(d_b);
184 if(d_x) cudaFree(d_x);
185
186 return 0 == code ? SUANPAN_SUCCESS : SUANPAN_FAIL;
187}
188
189#endif
190
191#endif
192
A SparseMatCUDA class that holds matrices.
A SparseMat class that holds matrices.
Definition: SparseMat.hpp:34
Definition: csr_form.hpp:25
unique_ptr< Material > make_copy(const shared_ptr< Material > &)
Definition: Material.cpp:370
double norm(const vec &)
Definition: tensor.cpp:370
#define suanpan_debug(...)
Definition: suanPan.h:307
constexpr auto SUANPAN_SUCCESS
Definition: suanPan.h:172
constexpr auto SUANPAN_FAIL
Definition: suanPan.h:173