QUDA  v1.1.0
A library for QCD on GPUs
blas_magma.cu
Go to the documentation of this file.
1 #include <blas_magma.h>
2 #include <string.h>
3 
4 #include <util_quda.h>
5 #include <quda_internal.h>
6 
7 #ifdef MAGMA_LIB
8 
9 #include <magma.h>
10 
11 #define _cV MagmaVec
12 #define _cU MagmaUpper
13 #define _cR MagmaRight
14 #define _cL MagmaLeft
15 #define _cC MagmaConjTrans
16 #define _cN MagmaNoTrans
17 #define _cNV MagmaNoVec
18 
19 #endif
20 
21 #ifdef MAGMA_LIB
22 
23 
24  template<typename magmaFloat> void magma_gesv(void *sol, const int ldn, const int n, void *Mat, const int ldm)
25  {
26  cudaPointerAttributes ptr_attr;
27  if(cudaPointerGetAttributes(&ptr_attr, Mat) == cudaErrorInvalidValue) errorQuda("In magma_gesv, a pointer was not allocated in, mapped by or registered with current CUDA context.\n");
28 
29  magma_int_t *ipiv;
30  magma_int_t err, info;
31 
32  magma_imalloc_pinned(&ipiv, n);
33 
34  void *tmp;
35 
36  magma_malloc_pinned((void**)&tmp, ldm*n*sizeof(magmaFloat));
37  memcpy(tmp, Mat, ldm*n*sizeof(magmaFloat));
38 
39  if ( ptr_attr.memoryType == cudaMemoryTypeDevice ) {
40  if(sizeof(magmaFloat) == sizeof(magmaFloatComplex))
41  {
42  err = magma_cgesv_gpu(n, 1, static_cast<magmaFloatComplex* >(tmp), ldm, ipiv, static_cast<magmaFloatComplex* >(sol), ldn, &info);
43  if(err != 0) errorQuda("\nError in SolveGPUProjMatrix (magma_cgesv_gpu), exit ...\n");
44  }
45  else
46  {
47  err = magma_zgesv_gpu(n, 1, static_cast<magmaDoubleComplex*>(tmp), ldm, ipiv, static_cast<magmaDoubleComplex*>(sol), ldn, &info);
48  if(err != 0) errorQuda("\nError in SolveGPUProjMatrix (magma_zgesv_gpu), exit ...\n");
49  }
50  } else if ( ptr_attr.memoryType == cudaMemoryTypeHost ) {
51 
52  if(sizeof(magmaFloat) == sizeof(magmaFloatComplex))
53  {
54  err = magma_cgesv(n, 1, static_cast<magmaFloatComplex* >(tmp), ldm, ipiv, static_cast<magmaFloatComplex* >(sol), ldn, &info);
55  if(err != 0) errorQuda("\nError in SolveGPUProjMatrix (magma_cgesv), exit ...\n");
56  }
57  else
58  {
59  err = magma_zgesv(n, 1, static_cast<magmaDoubleComplex*>(tmp), ldm, ipiv, static_cast<magmaDoubleComplex*>(sol), ldn, &info);
60  if(err != 0) errorQuda("\nError in SolveGPUProjMatrix (magma_zgesv), exit ...\n");
61  }
62  }
63 
64  magma_free_pinned(ipiv);
65  magma_free_pinned(tmp);
66 
67  return;
68  }
69 
70 
71 /////
72 
73  template<typename magmaFloat> void magma_geev(void *Mat, const int m, const int ldm, void *vr, void *evalues, const int ldv)
74  {
75  cudaPointerAttributes ptr_attr;
76  if(cudaPointerGetAttributes(&ptr_attr, Mat) == cudaErrorInvalidValue) errorQuda("In magma_geev, a pointer was not allocated in, mapped by or registered with current CUDA context.\n");
77 
78  magma_int_t err, info;
79 
80  void *work_ = nullptr, *rwork_ = nullptr;
81 
82  if ( ptr_attr.memoryType == cudaMemoryTypeDevice ) {
83  errorQuda("\nGPU version is not supported.\n");
84  } else if ( ptr_attr.memoryType == cudaMemoryTypeHost ) {
85 
86  if(sizeof(magmaFloat) == sizeof(magmaFloatComplex))
87  {
88  magmaFloatComplex qwork;
89 
90  magmaFloatComplex *work = static_cast<magmaFloatComplex*>(work_);
91  float *rwork = static_cast<float*>(rwork_);
92 
93  err = magma_cgeev(_cNV, _cV, m, nullptr, ldm, nullptr, nullptr, ldv, nullptr, ldv, &qwork, -1, nullptr, &info);
94  if( err != 0 ) errorQuda( "Error: CGEEVX, info %d\n",info);
95 
96  magma_int_t lwork = static_cast<magma_int_t>( MAGMA_C_REAL(qwork));
97 
98  magma_smalloc_pinned(&rwork, 2*m);
99  magma_cmalloc_pinned(&work, lwork);
100 
101  err = magma_cgeev(_cNV, _cV, m, static_cast<magmaFloatComplex*>(Mat), ldm, static_cast<magmaFloatComplex*>(evalues), nullptr, ldv, static_cast<magmaFloatComplex*>(vr), ldv, work, lwork, rwork, &info);
102  if( err != 0 ) errorQuda( "Error: CGEEVX, info %d\n",info);
103 
104  }
105  else
106  {
107  magmaDoubleComplex qwork;
108 
109  magmaDoubleComplex *work = static_cast<magmaDoubleComplex*>(work_);
110  double *rwork = static_cast<double*>(rwork_);
111  err = magma_zgeev(_cNV, _cV, m, nullptr, ldm, nullptr, nullptr, ldv, nullptr, ldv, &qwork, -1, nullptr, &info);
112  if( err != 0 ) errorQuda( "Error: ZGEEVX, info %d\n",info);
113 
114  magma_int_t lwork = static_cast<magma_int_t>( MAGMA_Z_REAL(qwork));
115 
116  magma_dmalloc_pinned(&rwork, 2*m);
117  magma_zmalloc_pinned(&work, lwork);
118 
119  err = magma_zgeev(_cNV, _cV, m, static_cast<magmaDoubleComplex*>(Mat), ldm, static_cast<magmaDoubleComplex*>(evalues), nullptr, ldv, static_cast<magmaDoubleComplex*>(vr), ldv, work, lwork, rwork, &info);
120  if( err != 0 ) errorQuda( "Error: ZGEEVX, info %d\n",info);
121  }
122  }
123 
124  if(rwork_) magma_free_pinned(rwork_);
125  if(work_ ) magma_free_pinned(work_);
126 
127  return;
128  }
129 
130 
131 /////
132 
133  template<typename magmaFloat> void magma_gels(void *Mat, void *c, int rows, int cols, int ldm)
134  {
135  cudaPointerAttributes ptr_attr;
136  if(cudaPointerGetAttributes(&ptr_attr, Mat) == cudaErrorInvalidValue) errorQuda("In magma_gels, a pointer was not allocated in, mapped by or registered with current CUDA context.\n");
137 
138  magma_int_t err, info, lwork;
139  void *hwork_ = nullptr;
140 
141  if ( ptr_attr.memoryType == cudaMemoryTypeDevice )
142  {
143  if(sizeof(magmaFloat) == sizeof(magmaFloatComplex))
144  {
145  magma_int_t nb = magma_get_cgeqrf_nb( rows, cols );
146  lwork = std::max( cols*nb, 2*nb*nb );
147 
148  magmaFloatComplex *hwork = static_cast<magmaFloatComplex*>(hwork_);
149  magma_cmalloc_cpu( &hwork, lwork);
150 
151  err = magma_cgels_gpu( _cN, rows, cols, 1, static_cast<magmaFloatComplex*>(Mat), ldm, static_cast<magmaFloatComplex*>(c),
152  ldm, hwork, lwork, &info );
153  if (err != 0) errorQuda("\nError in magma_cgels_gpu, %d, exit ...\n", info);
154  } else {
155  magma_int_t nb = magma_get_zgeqrf_nb( rows, cols );
156 
157  lwork = std::max( cols*nb, 2*nb*nb );
158  magmaDoubleComplex *hwork = static_cast<magmaDoubleComplex*>(hwork_);
159  magma_zmalloc_cpu( &hwork, lwork);
160 
161  err = magma_zgels_gpu( _cN, rows, cols, 1, static_cast<magmaDoubleComplex*>(Mat), ldm, static_cast<magmaDoubleComplex*>(c),
162  ldm, hwork, lwork, &info );
163  if (err != 0) errorQuda("\nError in magma_zgels_gpu, %d, exit ...\n", info);
164  }
165  } else if ( ptr_attr.memoryType == cudaMemoryTypeHost ) {
166 
167 
168  if(sizeof(magmaFloat) == sizeof(magmaFloatComplex))
169  {
170  magma_int_t nb = magma_get_cgeqrf_nb( rows, cols );
171 
172  lwork = std::max( cols*nb, 2*nb*nb );
173  magmaFloatComplex *hwork = static_cast<magmaFloatComplex*>(hwork_);
174  magma_cmalloc_cpu( &hwork, lwork);
175 
176  err = magma_cgels( _cN, rows, cols, 1, static_cast<magmaFloatComplex*>(Mat), ldm, static_cast<magmaFloatComplex*>(c),
177  ldm, hwork, lwork, &info );
178  if (err != 0) errorQuda("\nError in magma_cgels_cpu, %d, exit ...\n", info);
179  } else {
180  magma_int_t nb = magma_get_zgeqrf_nb( rows, cols );
181 
182  lwork = std::max( cols*nb, 2*nb*nb );
183  magmaDoubleComplex *hwork = static_cast<magmaDoubleComplex*>(hwork_);
184  magma_zmalloc_cpu( &hwork, lwork);
185 
186  err = magma_zgels( _cN, rows, cols, 1, static_cast<magmaDoubleComplex*>(Mat), ldm, static_cast<magmaDoubleComplex*>(c),
187  ldm, hwork, lwork, &info );
188  if (err != 0) errorQuda("\nError in magma_zgels_cpu, %d, exit ...\n", info);
189  }
190  }
191 
192  if(hwork_) magma_free_cpu(hwork_);
193 
194  return;
195  }
196 
197 
198  template<typename magmaFloat> void magma_heev(void *Mat, const int m, const int ldm, void *evalues)
199  {
200  cudaPointerAttributes ptr_attr;
201  if(cudaPointerGetAttributes(&ptr_attr, Mat) == cudaErrorInvalidValue) errorQuda("In magma_heev, a pointer was not allocated in, mapped by or registered with current CUDA context.\n");
202 
203  magma_int_t err, info;
204 
205  void *work_ = nullptr, *rwork_ = nullptr;
206  int *iwork = nullptr;
207  int qiwork;
208 
209  if ( ptr_attr.memoryType == cudaMemoryTypeDevice ) {
210  errorQuda("\nGPU version is not supported.\n");
211  } else if ( ptr_attr.memoryType == cudaMemoryTypeHost ) {
212  if(sizeof(magmaFloat) == sizeof(magmaFloatComplex))
213  {
214  magmaFloatComplex qwork;
215  float qrwork;
216 
217  magmaFloatComplex *work = static_cast<magmaFloatComplex*>(work_);
218  float *rwork = static_cast<float*>(rwork_);
219 
220  err = magma_cheevd(_cV, _cU, m, nullptr, ldm, nullptr, &qwork, -1, &qrwork, -1, &qiwork, -1, &info);
221  if( err != 0 ) errorQuda( "Error: CHEEVD, info %d\n",info);
222 
223  magma_int_t lwork = static_cast<magma_int_t>( MAGMA_C_REAL(qwork));
224  magma_int_t lrwork = static_cast<magma_int_t>( qrwork );
225  magma_int_t liwork = static_cast<magma_int_t>( qiwork );
226 
227  magma_cmalloc_pinned(&work, lwork);
228  magma_smalloc_pinned(&rwork, lrwork);
229  magma_imalloc_pinned(&iwork, liwork);
230 
231  err = magma_cheevd(_cV, _cU, m, static_cast<magmaFloatComplex*>(Mat), ldm, static_cast<float*>(evalues), work, lwork, rwork, lrwork, iwork, liwork, &info);
232  if( err != 0 ) errorQuda( "Error: CHEEVD, info %d\n",info);
233  } else {
234  magmaDoubleComplex qwork;
235  double qrwork;
236 
237  magmaDoubleComplex *work = static_cast<magmaDoubleComplex*>(work_);
238  double *rwork = static_cast<double*>(rwork_);
239 
240  err = magma_zheevd(_cV, _cU, m, nullptr, ldm, nullptr, &qwork, -1, &qrwork, -1, &qiwork, -1, &info);
241  if( err != 0 ) errorQuda( "Error: ZHEEVD, info %d\n",info);
242 
243  magma_int_t lwork = static_cast<magma_int_t>( MAGMA_Z_REAL(qwork));
244  magma_int_t lrwork = static_cast<magma_int_t>( qrwork );
245  magma_int_t liwork = static_cast<magma_int_t>( qiwork );
246 
247  magma_zmalloc_pinned(&work, lwork);
248  magma_dmalloc_pinned(&rwork, lrwork);
249  magma_imalloc_pinned(&iwork, liwork);
250 
251  err = magma_zheevd(_cV, _cU, m, static_cast<magmaDoubleComplex*>(Mat), ldm, static_cast<double*>(evalues), work, lwork, rwork, lrwork, iwork, liwork, &info);
252  if( err != 0 ) errorQuda( "Error: ZHEEVD, info %d\n",info);
253  }
254  }
255 
256  if(rwork_) magma_free_pinned(rwork_);
257  if(work_ ) magma_free_pinned(work_);
258  if(iwork ) magma_free_pinned(iwork);
259 
260  return;
261  }
262 
263 #endif // MAGMA_LIB
264 
265  void magma_Xgesv(void* sol, const int ldn, const int n, void* Mat, const int ldm, const int prec)
266  {
267 #ifdef MAGMA_LIB
268  if (prec == sizeof(std::complex< double >)) magma_gesv<magmaDoubleComplex>(sol, ldn, n, Mat, ldm);
269  else if (prec == sizeof(std::complex< float >)) magma_gesv<magmaFloatComplex >(sol, ldn, n, Mat, ldm);
270  else errorQuda("\nPrecision is not supported.\n");
271 #endif
272  return;
273  }
274 
275  void magma_Xgeev(void *Mat, const int m, const int ldm, void *vr, void *evalues, const int ldv, const int prec)
276  {
277 #ifdef MAGMA_LIB
278  if (prec == sizeof(std::complex< double >)) magma_geev<magmaDoubleComplex>(Mat, m, ldm, vr, evalues, ldv);
279  else if (prec == sizeof(std::complex< float >)) magma_geev<magmaFloatComplex >(Mat, m, ldm, vr, evalues, ldv);
280  else errorQuda("\nPrecision is not supported.\n");
281 #endif
282  return;
283  }
284 
285 
286  void magma_Xgels(void *Mat, void *c, int rows, int cols, int ldm, const int prec)
287  {
288 #ifdef MAGMA_LIB
289  if (prec == sizeof(std::complex< double >)) magma_gels<magmaDoubleComplex>(Mat, c, rows, cols, ldm);
290  else if (prec == sizeof(std::complex< float >)) magma_gels<magmaFloatComplex >(Mat, c, rows, cols, ldm);
291  else errorQuda("\nPrecision is not supported.\n");
292 #endif
293  return;
294  }
295 
296  void magma_Xheev(void *Mat, const int m, const int ldm, void *evalues, const int prec)
297  {
298 #ifdef MAGMA_LIB
299  if (prec == sizeof(std::complex< double >)) magma_heev<magmaDoubleComplex>(Mat, m, ldm, evalues);
300  else if (prec == sizeof(std::complex< float >)) magma_heev<magmaFloatComplex >(Mat, m, ldm, evalues);
301  else errorQuda("\nPrecision is not supported.\n");
302 #endif
303  return;
304  }
305 
306 
307  void OpenMagma(){
308 #ifdef MAGMA_LIB
309  magma_int_t err = magma_init();
310 
311  if(err != MAGMA_SUCCESS) errorQuda("\nError: cannot initialize MAGMA library\n");
312 
313  int major, minor, micro;
314 
315  magma_version( &major, &minor, &micro);
316  if (getVerbosity() >= QUDA_VERBOSE) printfQuda("\nMAGMA library version: %d.%d\n\n", major, minor);
317 #else
318  errorQuda("\nError: MAGMA library was not compiled, check your compilation options...\n");
319 #endif
320  return;
321  }
322 
323  void CloseMagma(){
324 #ifdef MAGMA_LIB
325  if(magma_finalize() != MAGMA_SUCCESS) errorQuda("\nError: cannot close MAGMA library\n");
326 #else
327  errorQuda("\nError: MAGMA library was not compiled, check your compilation options...\n");
328 #endif
329  return;
330  }
331 
332 
333 #ifdef MAGMA_LIB
334 
335 #undef _cV
336 #undef _cU
337 #undef _cR
338 #undef _cL
339 #undef _cC
340 #undef _cN
341 #undef _cNV
342 
343 #endif
Generated on Thu Oct 28 2021 16:10:26 for QUDA by doxygen 1.9.1