StatMech
Loading...
Searching...
No Matches
magma_wrapper.hpp
Go to the documentation of this file.
1#pragma once
2
3#if defined(GPU) && __has_include(<magma_v2.h>)
4 #include "Headers/mytypes.hpp"
5 #include <magma_v2.h>
6 #include <magma_operators.h>
7
8// MAGMA_BLAS
9template<typename ScalarType>
10static inline void magma_axpy_wrapper(magma_int_t n, ScalarType alpha, ScalarType const* dx,
11 magma_int_t incx, ScalarType* dy, magma_int_t incy,
12 magma_queue_t queue) {
13 if constexpr(std::is_same_v<ScalarType, Complex_t<float>>) {
14 magma_caxpy(n, alpha, reinterpret_cast<magmaFloatComplex_const_ptr>(dx), incx,
15 reinterpret_cast<magmaFloatComplex_const_ptr>(dy), incy, queue);
16 }
17 else if constexpr(std::is_same_v<ScalarType, Complex_t<double>>) {
18 magma_zaxpy(n, alpha, reinterpret_cast<magmaDoubleComplex_const_ptr>(dx), incx,
19 reinterpret_cast<magmaDoubleComplex_const_ptr>(dy), incy, queue);
20 }
21 else {
22 static_assert([] { return false; }(),
23 "Error: magma_axpy_wrapper for the input scalar type is NOT implemented.");
24 }
25}
26
27// static inline magma_int_t magmablas_hemv_work_wrapper
28// (magma_uplo_t Uplo, magma_int_t n, Complex_t<float> alpha,
29// matrix_gpu<Complex_t<float>> const& dA,
30// Complex_t<float> const* dx, magma_int_t incx, Complex_t<float> beta,
31// Complex_t<float>* dy, magma_int_t incy,
32// Complex_t<float>* dwork, magma_int_t lwork, magma_queue_t queue)
33// {
34// return magmablas_chemv_work(Uplo, n, alpha, dA.ptr(), dA.LD(), dx, incx, beta, dy, incy, dwork, lwork, queue);
35// }
36//
37// static inline magma_int_t magmablas_hemv_work_wrapper
38// (magma_uplo_t Uplo, magma_int_t n, Complex_t<double> alpha,
39// matrix_gpu<Complex_t<double>> const& dA,
40// Complex_t<double> const* dx, magma_int_t incx, Complex_t<double> beta,
41// Complex_t<double>* dy, magma_int_t incy,
42// Complex_t<double>* dwork, magma_int_t lwork, magma_queue_t queue)
43// {
44// return magmablas_zhemv_work(Uplo, n, alpha, dA.ptr(), dA.LD(), dx, incx, beta, dy, incy, dwork, lwork, queue);
45// }
46
47template<typename ScalarType>
48static inline ScalarType magma_dotc_wrapper(magma_int_t n, ScalarType const* dx, magma_int_t incx,
49 ScalarType const* dy, magma_int_t incy,
50 magma_queue_t queue) {
51 if constexpr(std::is_same_v<ScalarType, float>) {
52 return magma_sdot(n, dx, incx, dy, incy, queue);
53 }
54 else if constexpr(std::is_same_v<ScalarType, double>) {
55 return magma_ddot(n, dx, incx, dy, incy, queue);
56 }
57 else if constexpr(std::is_same_v<ScalarType, Complex_t<float>>) {
58 return magma_cdotc(n, reinterpret_cast<magmaFloatComplex_const_ptr>(dx), incx,
59 reinterpret_cast<magmaFloatComplex_const_ptr>(dy), incy, queue);
60 }
61 else if constexpr(std::is_same_v<ScalarType, Complex_t<double>>) {
62 return magma_zdotc(n, reinterpret_cast<magmaDoubleComplex_const_ptr>(dx), incx,
63 reinterpret_cast<magmaDoubleComplex_const_ptr>(dy), incy, queue);
64 }
65 else {
66 static_assert([] { return false; }(),
67 "Error: magma_dotc_wrapper for the input scalar type is NOT implemented.");
68 }
69 // To suppress wrong compiler warning "warning: missing return statement at end of non-void function"
70 // Never be executed
71 return 0;
72}
73
74template<typename ScalarType>
75static inline void magma_hemm_wrapper(magma_side_t side, magma_uplo_t uplo, magma_int_t m,
76 magma_int_t n, ScalarType alpha, ScalarType const* dA,
77 magma_int_t ldda, ScalarType const* dB, magma_int_t lddb,
78 ScalarType beta, ScalarType* dC, magma_int_t lddc,
79 magma_queue_t queue) {
80 if constexpr(std::is_same_v<ScalarType, float>) {
81 return magmablas_ssymm(side, uplo, m, n, alpha, dA, ldda, dB, lddb, beta, dC, lddc, queue);
82 }
83 else if constexpr(std::is_same_v<ScalarType, double>) {
84 return magmablas_dsymm(side, uplo, m, n, alpha, dA, ldda, dB, lddb, beta, dC, lddc, queue);
85 }
86 else if constexpr(std::is_same_v<ScalarType, Complex_t<float>>) {
87 return magma_chemm(side, uplo, m, n, alpha,
88 reinterpret_cast<magmaFloatComplex_const_ptr>(dA), ldda,
89 reinterpret_cast<magmaFloatComplex_const_ptr>(dB), lddb, beta,
90 reinterpret_cast<magmaFloatComplex_ptr>(dC), lddc, queue);
91 }
92 else if constexpr(std::is_same_v<ScalarType, Complex_t<double>>) {
93 return magma_zhemm(side, uplo, m, n, alpha,
94 reinterpret_cast<magmaDoubleComplex_const_ptr>(dA), ldda,
95 reinterpret_cast<magmaDoubleComplex_const_ptr>(dB), lddb, beta,
96 reinterpret_cast<magmaDoubleComplex_ptr>(dC), lddc, queue);
97 }
98 else {
99 static_assert([] { return false; }(),
100 "Error: magma_hemm_wrapper for the input scalar type is NOT implemented.");
101 }
102 return;
103}
104
105// MAGMA_LAPACK
106template<typename ScalarType>
107static inline magma_int_t magma_heevd_gpu_wrapper(
108 magma_vec_t jobz, magma_uplo_t uplo, magma_int_t n, ScalarType* dA, magma_int_t ldda,
109 typename Eigen::NumTraits<ScalarType>::Real* w, ScalarType* wA, magma_int_t ldwa,
110 ScalarType* work, magma_int_t lwork, typename Eigen::NumTraits<ScalarType>::Real* rwork,
111 magma_int_t lrwork, magma_int_t* iwork, magma_int_t liwork, magma_int_t* info) {
112 if constexpr(std::is_same_v<ScalarType, float>) {
113 return magma_sheevd_gpu(jobz, uplo, n, dA, ldda, w, wA, ldwa, work, lwork, rwork, lrwork,
114 iwork, liwork, info);
115 }
116 else if constexpr(std::is_same_v<ScalarType, double>) {
117 return magma_dheevd_gpu(jobz, uplo, n, dA, ldda, w, wA, ldwa, work, lwork, rwork, lrwork,
118 iwork, liwork, info);
119 }
120 else if constexpr(std::is_same_v<ScalarType, Complex_t<float>>) {
121 return magma_cheevd_gpu(jobz, uplo, n, reinterpret_cast<magmaFloatComplex_ptr>(dA), ldda, w,
122 reinterpret_cast<magmaFloatComplex_ptr>(wA), ldwa,
123 reinterpret_cast<magmaFloatComplex_ptr>(work), lwork, rwork, lrwork,
124 iwork, liwork, info);
125 }
126 else if constexpr(std::is_same_v<ScalarType, Complex_t<double>>) {
127 return magma_zheevd_gpu(jobz, uplo, n, reinterpret_cast<magmaDoubleComplex_ptr>(dA), ldda,
128 w, reinterpret_cast<magmaDoubleComplex_ptr>(wA), ldwa,
129 reinterpret_cast<magmaDoubleComplex_ptr>(work), lwork, rwork,
130 lrwork, iwork, liwork, info);
131 }
132 else {
133 static_assert(
134 [] { return false; }(),
135 "Error: magma_heevd_gpu_wrapper for the input scalar type is NOT implemented.");
136 }
137 // To suppress wrong compiler warning "warning: missing return statement at end of non-void function"
138 // Never be executed
139 return EXIT_SUCCESS;
140}
141
142// template<class MatrixGPU_t, class VectorCPU_t,
143// typename std::enable_if_t<on_GPU_v<MatrixGPU_t>>* = nullptr >
144// static inline magma_int_t magma_heevd_wrapper(magma_vec_t Job, MatrixGPU_t& dMatrix,
145// VectorCPU_t& hEigVal,
146// std::vector<typename MatrixGPU_t::Scalar> wA,
147// std::vector<typename MatrixGPU_t::Scalar> work,
148// std::vector<typename MatrixGPU_t::RealScalar> rwork,
149// std::vector<magma_int_t> iwork) {
150// assert(dMatrix.rows() == dMatrix.cols());
151// magma_int_t info, lwork = -1, lrwork = -1, liwork = -1;
152// magma_int_t dim = dMatrix.rows();
153// hEigVal.resize(dim);
154// wA.resize(dim * dim);
155// work.resize(1);
156// rwork.resize(1);
157// iwork.resize(1);
158// //********** Query for work space sizes (lwork, lrwork, liwork) **********//
159// if constexpr(std::is_same_v<typename MatrixGPU_t::Scalar, Complex_t<float>>) {
160// static_assert(std::is_same_v<typename MatrixGPU_tRealScalaralar, float>,
161// "MatrixGPURealScalarScalar is not float.");
162// debug_print("# " << __func__ << "(c): Job=" << Job << ", dim=" << dim);
163// magma_cheevd_gpu(Job, MagmaLower, dim,
164// reinterpret_cast<magmaFloatComplex_ptr>(dMatrix.data()), dMatrix.LD(),
165// (float*)&*hEigVal.begin(),
166// reinterpret_cast<magmaFloatComplex_ptr>(&*wA.begin()), dim,
167// reinterpret_cast<magmaFloatComplex_ptr>(&*work.begin()), lwork,
168// &*rwork.begin(), lrwork, &liwork, liwork, &info);
169// }
170// else if constexpr(std::is_same_v<typename MatrixGPU_t::Scalar, Complex_t<double>>) {
171// static_assert(std::is_same_v<typename MatrixGRealScalaralScalar, double>,
172// "MatriRealScalarRealScalar is not double.");
173// debug_print("# " << __func__ << "(z): Job=" << Job << ", dim=" << dim);
174// magma_zheevd_gpu(Job, MagmaLower, dim,
175// reinterpret_cast<magmaDoubleComplex_ptr>(dMatrix.data()), dMatrix.LD(),
176// &*hEigVal.begin(), reinterpret_cast<magmaDoubleComplex_ptr>(&*wA.begin()),
177// dim, reinterpret_cast<magmaDoubleComplex_ptr>(&*work.begin()), lwork,
178// &*rwork.begin(), lrwork, &liwork, liwork, &info);
179// }
180// else {
181// static_assert([] { return false; }(),
182// "Error: magma_heevd_wrapper for the input scalar type is NOT implemented.");
183// }
184// lwork = (magma_int_t)real(work[0]);
185// lrwork = (magma_int_t)rwork[0];
186// work.resize(lwork);
187// rwork.resize(lrwork);
188// iwork.resize(liwork);
189// //********** (END)Query for work space sizes (lwork, lrwork, liwork) **********//
190
191// if constexpr(std::is_same_v< typename MatrixGPU_t::Scalar, Complex_t<float>>) {
192// Eigen::VectorX<float> hEigVal_f(dim);
193// magma_cheevd_gpu(Job, MagmaLower, dim,
194// reinterpret_cast<magmaFloatComplex_ptr>(dMatrix.data()), dMatrix.LD(),
195// &*hEigVal_f.begin(), reinterpret_cast<magmaFloatComplex_ptr>(&*wA.begin()),
196// dim, reinterpret_cast<magmaFloatComplex_ptr>(&*work.begin()), lwork,
197// &*rwork.begin(), lrwork, &*iwork.begin(), liwork, &info);
198// std::copy(hEigVal_f.begin(), hEigVal_f.end(), hEigVal.begin());
199// }
200// else if constexpr(std::is_same_v< typename MatrixGPU_t::Scalar, Complex_t<double>>) {
201// magma_zheevd_gpu(Job, MagmaLower, dim,
202// reinterpret_cast<magmaDoubleComplex_ptr>(dMatrix.data()), dMatrix.LD(),
203// &*hEigVal.begin(), reinterpret_cast<magmaDoubleComplex_ptr>(&*wA.begin()),
204// dim, reinterpret_cast<magmaDoubleComplex_ptr>(&*work.begin()), lwork,
205// &*rwork.begin(), lrwork, &*iwork.begin(), liwork, &info);
206// }
207
208// if(info != 0) {
209// if(dim < 20) std::cout << dMatrix << std::endl;
210// std::cerr << "# Error: magma_cheevd_gpu failed. (info=" << info << ")" << std::endl;
211// std::cout << "# Inputs: 3.8.dim=" << dim << ", 5.LDH=" << dMatrix.LD()
212// << ", 10.lwork=" << lwork << ", 12.lrwork=" << lrwork << ", 13.liwork=" << liwork
213// << std::endl;
214// std::cerr << "# ";
215// magma_xerbla("magma_cheevd_gpu", -info);
216// std::cerr << "# " << magma_strerror(info) << "\n" << std::endl;
217// std::exit(EXIT_FAILURE);
218// }
219// return info;
220// }
221
222// template<class MatrixGPU_t, class VectorCPU_t>
223// static inline magma_int_t magma_heevd_wrapper(magma_vec_t Job, MatrixGPU_t& dMatrix,
224// VectorCPU_t& hEigVal) {
225// using Scalar = typename MatrixGPU_t::Scalar;
226// using RealScalar = Eigen::NumTraits<Scalar>::Real;
227// std::vector<Scalar> wA;
228// std::vector<Scalar> work;
229// std::vector<RealScalar> rwork;
230// std::vector<magma_int_t> iwork;
231// return magma_heevd_wrapper(Job, dMatrix, hEigVal, wA, work, rwork, iwork);
232// }
233
234#endif