Reference documentation for deal.II version GIT e8f57130d0 2023-09-28 07:00:03+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
cuda_precondition.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2018 - 2023 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
18 
20 
21 namespace
22 {
29  template <typename Number>
30  cusparseStatus_t
31  cusparseXcsrilu02(cusparseHandle_t /*handle*/,
32  int /*m*/,
33  int /*nnz*/,
34  const cusparseMatDescr_t /*descrA*/,
35  Number * /*csrValA_valM*/,
36  const int * /*csrRowPtrA*/,
37  const int * /*csrColIndA*/,
38  csrilu02Info_t /*info*/,
39  cusparseSolvePolicy_t /*policy*/,
40  void * /*pBuffer*/)
41  {
43  return CUSPARSE_STATUS_INVALID_VALUE;
44  }
45 
46  template <>
47  cusparseStatus_t
48  cusparseXcsrilu02<float>(cusparseHandle_t handle,
49  int m,
50  int nnz,
51  const cusparseMatDescr_t descrA,
52  float *csrValA_valM,
53  const int *csrRowPtrA,
54  const int *csrColIndA,
55  csrilu02Info_t info,
56  cusparseSolvePolicy_t policy,
57  void *pBuffer)
58  {
59  return cusparseScsrilu02(handle,
60  m,
61  nnz,
62  descrA,
63  csrValA_valM,
64  csrRowPtrA,
65  csrColIndA,
66  info,
67  policy,
68  pBuffer);
69  }
70 
71  template <>
72  cusparseStatus_t
73  cusparseXcsrilu02<double>(cusparseHandle_t handle,
74  int m,
75  int nnz,
76  const cusparseMatDescr_t descrA,
77  double *csrValA_valM,
78  const int *csrRowPtrA,
79  const int *csrColIndA,
80  csrilu02Info_t info,
81  cusparseSolvePolicy_t policy,
82  void *pBuffer)
83  {
84  return cusparseDcsrilu02(handle,
85  m,
86  nnz,
87  descrA,
88  csrValA_valM,
89  csrRowPtrA,
90  csrColIndA,
91  info,
92  policy,
93  pBuffer);
94  }
95 
96  /*
97  template <>
98  cusparseStatus_t
99  cusparseXcsrilu02<cuComplex>(cusparseHandle_t handle,
100  int m,
101  int nnz,
102  const cusparseMatDescr_t descrA,
103  cuComplex * csrValA_valM,
104  const int * csrRowPtrA,
105  const int * csrColIndA,
106  csrilu02Info_t info,
107  cusparseSolvePolicy_t policy,
108  void * pBuffer)
109  {
110  return cusparseCcsrilu02(handle,
111  m,
112  nnz,
113  descrA,
114  csrValA_valM,
115  csrRowPtrA,
116  csrColIndA,
117  info,
118  policy,
119  pBuffer);
120  }
121 
122  template <>
123  cusparseStatus_t
124  cusparseXcsrilu02<cuDoubleComplex>(cusparseHandle_t handle,
125  int m,
126  int nnz,
127  const cusparseMatDescr_t descrA,
128  cuDoubleComplex * csrValA_valM,
129  const int * csrRowPtrA,
130  const int * csrColIndA,
131  csrilu02Info_t info,
132  cusparseSolvePolicy_t policy,
133  void * pBuffer)
134  {
135  return cusparseZcsrilu02(handle,
136  m,
137  nnz,
138  descrA,
139  csrValA_valM,
140  csrRowPtrA,
141  csrColIndA,
142  info,
143  policy,
144  pBuffer);
145  }
146  */
147 
148 
149 
156  template <typename Number>
157  cusparseStatus_t
158  cusparseXcsrilu02_analysis(cusparseHandle_t /*handle*/,
159  int /*m*/,
160  int /*nnz*/,
161  const cusparseMatDescr_t /*descrA*/,
162  const Number * /*csrValA*/,
163  const int * /*csrRowPtrA*/,
164  const int * /*csrColIndA*/,
165  csrilu02Info_t /*info*/,
166  cusparseSolvePolicy_t /*policy*/,
167  void * /*pBuffer*/)
168  {
169  AssertThrow(false, ExcNotImplemented());
170  return CUSPARSE_STATUS_INVALID_VALUE;
171  }
172 
173  template <>
174  cusparseStatus_t
175  cusparseXcsrilu02_analysis<float>(cusparseHandle_t handle,
176  int m,
177  int nnz,
178  const cusparseMatDescr_t descrA,
179  const float *csrValA,
180  const int *csrRowPtrA,
181  const int *csrColIndA,
182  csrilu02Info_t info,
183  cusparseSolvePolicy_t policy,
184  void *pBuffer)
185  {
186  return cusparseScsrilu02_analysis(handle,
187  m,
188  nnz,
189  descrA,
190  csrValA,
191  csrRowPtrA,
192  csrColIndA,
193  info,
194  policy,
195  pBuffer);
196  }
197 
198  template <>
199  cusparseStatus_t
200  cusparseXcsrilu02_analysis<double>(cusparseHandle_t handle,
201  int m,
202  int nnz,
203  const cusparseMatDescr_t descrA,
204  const double *csrValA,
205  const int *csrRowPtrA,
206  const int *csrColIndA,
207  csrilu02Info_t info,
208  cusparseSolvePolicy_t policy,
209  void *pBuffer)
210  {
211  return cusparseDcsrilu02_analysis(handle,
212  m,
213  nnz,
214  descrA,
215  csrValA,
216  csrRowPtrA,
217  csrColIndA,
218  info,
219  policy,
220  pBuffer);
221  }
222 
223  /*
224  template <>
225  cusparseStatus_t
226  cusparseXcsrilu02_analysis<cuComplex>(cusparseHandle_t handle,
227  int m,
228  int nnz,
229  const cusparseMatDescr_t descrA,
230  const cuComplex * csrValA,
231  const int * csrRowPtrA,
232  const int * csrColIndA,
233  csrilu02Info_t info,
234  cusparseSolvePolicy_t policy,
235  void * pBuffer)
236  {
237  return cusparseCcsrilu02_analysis(handle,
238  m,
239  nnz,
240  descrA,
241  csrValA,
242  csrRowPtrA,
243  csrColIndA,
244  info,
245  policy,
246  pBuffer);
247  }
248 
249  template <>
250  cusparseStatus_t
251  cusparseXcsrilu02_analysis<cuDoubleComplex>(cusparseHandle_t handle,
252  int m,
253  int nnz,
254  const cusparseMatDescr_t descrA,
255  const cuDoubleComplex * csrValA,
256  const int * csrRowPtrA,
257  const int * csrColIndA,
258  csrilu02Info_t info,
259  cusparseSolvePolicy_t policy,
260  void * pBuffer)
261  {
262  return cusparseZcsrilu02_analysis(handle,
263  m,
264  nnz,
265  descrA,
266  csrValA,
267  csrRowPtrA,
268  csrColIndA,
269  info,
270  policy,
271  pBuffer);
272  }
273 */
274 
275 
276 
283  template <typename Number>
284  cusparseStatus_t
285  cusparseXcsrilu02_bufferSize(cusparseHandle_t /*handle*/,
286  int /*m*/,
287  int /*nnz*/,
288  const cusparseMatDescr_t /*descrA*/,
289  Number * /*csrValA*/,
290  const int * /*csrRowPtrA*/,
291  const int * /*csrColIndA*/,
292  csrilu02Info_t /*info*/,
293  int * /*pBufferSizeInBytes*/)
294  {
295  AssertThrow(false, ExcNotImplemented());
296  return CUSPARSE_STATUS_INVALID_VALUE;
297  }
298 
299  template <>
300  cusparseStatus_t
301  cusparseXcsrilu02_bufferSize<float>(cusparseHandle_t handle,
302  int m,
303  int nnz,
304  const cusparseMatDescr_t descrA,
305  float *csrValA,
306  const int *csrRowPtrA,
307  const int *csrColIndA,
308  csrilu02Info_t info,
309  int *pBufferSizeInBytes)
310  {
311  return cusparseScsrilu02_bufferSize(handle,
312  m,
313  nnz,
314  descrA,
315  csrValA,
316  csrRowPtrA,
317  csrColIndA,
318  info,
319  pBufferSizeInBytes);
320  }
321 
322  template <>
323  cusparseStatus_t
324  cusparseXcsrilu02_bufferSize<double>(cusparseHandle_t handle,
325  int m,
326  int nnz,
327  const cusparseMatDescr_t descrA,
328  double *csrValA,
329  const int *csrRowPtrA,
330  const int *csrColIndA,
331  csrilu02Info_t info,
332  int *pBufferSizeInBytes)
333  {
334  return cusparseDcsrilu02_bufferSize(handle,
335  m,
336  nnz,
337  descrA,
338  csrValA,
339  csrRowPtrA,
340  csrColIndA,
341  info,
342  pBufferSizeInBytes);
343  }
344 
345  /*
346  template <>
347  cusparseStatus_t
348  cusparseXcsrilu02_bufferSize<cuComplex>(cusparseHandle_t handle,
349  int m,
350  int nnz,
351  const cusparseMatDescr_t descrA,
352  cuComplex * csrValA,
353  const int * csrRowPtrA,
354  const int * csrColIndA,
355  csrilu02Info_t info,
356  int *pBufferSizeInBytes)
357  {
358  return cusparseCcsrilu02_bufferSize(handle,
359  m,
360  nnz,
361  descrA,
362  csrValA,
363  csrRowPtrA,
364  csrColIndA,
365  info,
366  pBufferSizeInBytes);
367  }
368 
369  template <>
370  cusparseStatus_t
371  cusparseXcsrilu02_bufferSize<cuDoubleComplex>(cusparseHandle_t handle, int
372  m, int nnz, const cusparseMatDescr_t descrA,
373  cuDoubleComplex *csrValA,
374  const int * csrRowPtrA,
375  const int * csrColIndA,
376  csrilu02Info_t info,
377  int *pBufferSizeInBytes)
378  {
379  return cusparseZcsrilu02_bufferSize(handle,
380  m,
381  nnz,
382  descrA,
383  csrValA,
384  csrRowPtrA,
385  csrColIndA,
386  info,
387  pBufferSizeInBytes);
388  }
389 */
390 
391 
392 
399  template <typename Number>
400  cusparseStatus_t
401  cusparseXcsric02(cusparseHandle_t /*handle*/,
402  int /*m*/,
403  int /*nnz*/,
404  const cusparseMatDescr_t /*descrA*/,
405  Number * /*csrValA_valM*/,
406  const int * /*csrRowPtrA*/,
407  const int * /*csrColIndA*/,
408  csric02Info_t /*info*/,
409  cusparseSolvePolicy_t /*policy*/,
410  void * /*pBuffer*/)
411  {
412  AssertThrow(false, ExcNotImplemented());
413  return CUSPARSE_STATUS_INVALID_VALUE;
414  }
415 
416  template <>
417  cusparseStatus_t
418  cusparseXcsric02<float>(cusparseHandle_t handle,
419  int m,
420  int nnz,
421  const cusparseMatDescr_t descrA,
422  float *csrValA_valM,
423  const int *csrRowPtrA,
424  const int *csrColIndA,
425  csric02Info_t info,
426  cusparseSolvePolicy_t policy,
427  void *pBuffer)
428  {
429  return cusparseScsric02(handle,
430  m,
431  nnz,
432  descrA,
433  csrValA_valM,
434  csrRowPtrA,
435  csrColIndA,
436  info,
437  policy,
438  pBuffer);
439  }
440 
441  template <>
442  cusparseStatus_t
443  cusparseXcsric02<double>(cusparseHandle_t handle,
444  int m,
445  int nnz,
446  const cusparseMatDescr_t descrA,
447  double *csrValA_valM,
448  const int *csrRowPtrA,
449  const int *csrColIndA,
450  csric02Info_t info,
451  cusparseSolvePolicy_t policy,
452  void *pBuffer)
453  {
454  return cusparseDcsric02(handle,
455  m,
456  nnz,
457  descrA,
458  csrValA_valM,
459  csrRowPtrA,
460  csrColIndA,
461  info,
462  policy,
463  pBuffer);
464  }
465 
466  /*
467  template <>
468  cusparseStatus_t
469  cusparseXcsric02<cuComplex>(cusparseHandle_t handle,
470  int m,
471  int nnz,
472  const cusparseMatDescr_t descrA,
473  cuComplex * csrValA_valM,
474  const int * csrRowPtrA,
475  const int * csrColIndA,
476  csric02Info_t info,
477  cusparseSolvePolicy_t policy,
478  void * pBuffer)
479  {
480  return cusparseCcsric02(handle,
481  m,
482  nnz,
483  descrA,
484  csrValA_valM,
485  csrRowPtrA,
486  csrColIndA,
487  info,
488  policy,
489  pBuffer);
490  }
491 
492  template <>
493  cusparseStatus_t
494  cusparseXcsric02<cuDoubleComplex>(cusparseHandle_t handle,
495  int m,
496  int nnz,
497  const cusparseMatDescr_t descrA,
498  cuDoubleComplex * csrValA_valM,
499  const int * csrRowPtrA,
500  const int * csrColIndA,
501  csric02Info_t info,
502  cusparseSolvePolicy_t policy,
503  void * pBuffer)
504  {
505  return cusparseZcsric02(handle,
506  m,
507  nnz,
508  descrA,
509  csrValA_valM,
510  csrRowPtrA,
511  csrColIndA,
512  info,
513  policy,
514  pBuffer);
515  }
516  */
517 
518 
519 
526  template <typename Number>
527  cusparseStatus_t
528  cusparseXcsrsv2_solve(cusparseHandle_t /*handle*/,
529  cusparseOperation_t /*transA*/,
530  int /*m*/,
531  int /*nnz*/,
532  const Number * /*alpha*/,
533  const cusparseMatDescr_t /*descra*/,
534  const Number * /*csrValA*/,
535  const int * /*csrRowPtrA*/,
536  const int * /*csrColIndA*/,
537  csrsv2Info_t /*info*/,
538  const Number * /*x*/,
539  Number * /*y*/,
540  cusparseSolvePolicy_t /*policy*/,
541  void * /*pBuffer*/)
542  {
543  AssertThrow(false, ExcNotImplemented());
544  return CUSPARSE_STATUS_INVALID_VALUE;
545  }
546 
547  template <>
548  cusparseStatus_t
549  cusparseXcsrsv2_solve<float>(cusparseHandle_t handle,
550  cusparseOperation_t transA,
551  int m,
552  int nnz,
553  const float *alpha,
554  const cusparseMatDescr_t descra,
555  const float *csrValA,
556  const int *csrRowPtrA,
557  const int *csrColIndA,
558  csrsv2Info_t info,
559  const float *x,
560  float *y,
561  cusparseSolvePolicy_t policy,
562  void *pBuffer)
563  {
564  return cusparseScsrsv2_solve(handle,
565  transA,
566  m,
567  nnz,
568  alpha,
569  descra,
570  csrValA,
571  csrRowPtrA,
572  csrColIndA,
573  info,
574  x,
575  y,
576  policy,
577  pBuffer);
578  }
579 
580  template <>
581  cusparseStatus_t
582  cusparseXcsrsv2_solve<double>(cusparseHandle_t handle,
583  cusparseOperation_t transA,
584  int m,
585  int nnz,
586  const double *alpha,
587  const cusparseMatDescr_t descra,
588  const double *csrValA,
589  const int *csrRowPtrA,
590  const int *csrColIndA,
591  csrsv2Info_t info,
592  const double *x,
593  double *y,
594  cusparseSolvePolicy_t policy,
595  void *pBuffer)
596  {
597  return cusparseDcsrsv2_solve(handle,
598  transA,
599  m,
600  nnz,
601  alpha,
602  descra,
603  csrValA,
604  csrRowPtrA,
605  csrColIndA,
606  info,
607  x,
608  y,
609  policy,
610  pBuffer);
611  }
612 
613  /*
614  template <>
615  cusparseStatus_t
616  cusparseXcsrsv2_solve<cuComplex>(cusparseHandle_t handle,
617  cusparseOperation_t transA,
618  int m,
619  int nnz,
620  const cuComplex * alpha,
621  const cusparseMatDescr_t descra,
622  const cuComplex * csrValA,
623  const int * csrRowPtrA,
624  const int * csrColIndA,
625  csrsv2Info_t info,
626  const cuComplex * x,
627  cuComplex * y,
628  cusparseSolvePolicy_t policy,
629  void * pBuffer)
630  {
631  return cusparseCcsrsv2_solve(handle,
632  transA,
633  m,
634  nnz,
635  alpha,
636  descra,
637  csrValA,
638  csrRowPtrA,
639  csrColIndA,
640  info,
641  x,
642  y,
643  policy,
644  pBuffer);
645  }
646 
647  template <>
648  cusparseStatus_t
649  cusparseXcsrsv2_solve<cuDoubleComplex>(cusparseHandle_t handle,
650  cusparseOperation_t transA,
651  int m,
652  int nnz,
653  const cuDoubleComplex * alpha,
654  const cusparseMatDescr_t descra,
655  const cuDoubleComplex * csrValA,
656  const int * csrRowPtrA,
657  const int * csrColIndA,
658  csrsv2Info_t info,
659  const cuDoubleComplex * x,
660  cuDoubleComplex * y,
661  cusparseSolvePolicy_t policy,
662  void * pBuffer)
663  {
664  return cusparseZcsrsv2_solve(handle,
665  transA,
666  m,
667  nnz,
668  alpha,
669  descra,
670  csrValA,
671  csrRowPtrA,
672  csrColIndA,
673  info,
674  x,
675  y,
676  policy,
677  pBuffer);
678  }
679 */
680 
681 
682 
689  template <typename Number>
690  cusparseStatus_t
691  cusparseXcsrsv2_analysis(cusparseHandle_t /*handle*/,
692  cusparseOperation_t /*transA*/,
693  int /*m*/,
694  int /*nnz*/,
695  const cusparseMatDescr_t /*descrA*/,
696  const Number * /*csrValA*/,
697  const int * /*csrRowPtrA*/,
698  const int * /*csrColIndA*/,
699  csrsv2Info_t /*info*/,
700  cusparseSolvePolicy_t /*policy*/,
701  void * /*pBuffer*/)
702  {
703  AssertThrow(false, ExcNotImplemented());
704  return CUSPARSE_STATUS_INVALID_VALUE;
705  }
706 
707  template <>
708  cusparseStatus_t
709  cusparseXcsrsv2_analysis<float>(cusparseHandle_t handle,
710  cusparseOperation_t transA,
711  int m,
712  int nnz,
713  const cusparseMatDescr_t descrA,
714  const float *csrValA,
715  const int *csrRowPtrA,
716  const int *csrColIndA,
717  csrsv2Info_t info,
718  cusparseSolvePolicy_t policy,
719  void *pBuffer)
720  {
721  return cusparseScsrsv2_analysis(handle,
722  transA,
723  m,
724  nnz,
725  descrA,
726  csrValA,
727  csrRowPtrA,
728  csrColIndA,
729  info,
730  policy,
731  pBuffer);
732  }
733 
734  template <>
735  cusparseStatus_t
736  cusparseXcsrsv2_analysis<double>(cusparseHandle_t handle,
737  cusparseOperation_t transA,
738  int m,
739  int nnz,
740  const cusparseMatDescr_t descrA,
741  const double *csrValA,
742  const int *csrRowPtrA,
743  const int *csrColIndA,
744  csrsv2Info_t info,
745  cusparseSolvePolicy_t policy,
746  void *pBuffer)
747  {
748  return cusparseDcsrsv2_analysis(handle,
749  transA,
750  m,
751  nnz,
752  descrA,
753  csrValA,
754  csrRowPtrA,
755  csrColIndA,
756  info,
757  policy,
758  pBuffer);
759  }
760 
761  /*
762  template <>
763  cusparseStatus_t
764  cusparseXcsrsv2_analysis<cuComplex>(cusparseHandle_t handle,
765  cusparseOperation_t transA,
766  int m,
767  int nnz,
768  const cusparseMatDescr_t descrA,
769  const cuComplex * csrValA,
770  const int * csrRowPtrA,
771  const int * csrColIndA,
772  csrsv2Info_t info,
773  cusparseSolvePolicy_t policy,
774  void * pBuffer)
775  {
776  return cusparseCcsrsv2_analysis(handle,
777  transA,
778  m,
779  nnz,
780  descrA,
781  csrValA,
782  csrRowPtrA,
783  csrColIndA,
784  info,
785  policy,
786  pBuffer);
787  }
788 
789  template <>
790  cusparseStatus_t
791  cusparseXcsrsv2_analysis<cuDoubleComplex>(cusparseHandle_t handle,
792  cusparseOperation_t transA,
793  int m,
794  int nnz,
795  const cusparseMatDescr_t descrA,
796  const cuDoubleComplex * csrValA,
797  const int * csrRowPtrA,
798  const int * csrColIndA,
799  csrsv2Info_t info,
800  cusparseSolvePolicy_t policy,
801  void * pBuffer)
802  {
803  return cusparseZcsrsv2_analysis(handle,
804  transA,
805  m,
806  nnz,
807  descrA,
808  csrValA,
809  csrRowPtrA,
810  csrColIndA,
811  info,
812  policy,
813  pBuffer);
814  }
815 */
816 
817 
824  template <typename Number>
825  cusparseStatus_t
826  cusparseXcsric02_analysis(cusparseHandle_t /*handle*/,
827  int /*m*/,
828  int /*nnz*/,
829  const cusparseMatDescr_t /*descrA*/,
830  const Number * /*csrValA*/,
831  const int * /*csrRowPtrA*/,
832  const int * /*csrColIndA*/,
833  csric02Info_t /*info*/,
834  cusparseSolvePolicy_t /*policy*/,
835  void * /*pBuffer*/)
836  {
837  AssertThrow(false, ExcNotImplemented());
838  return CUSPARSE_STATUS_INVALID_VALUE;
839  }
840 
841  template <>
842  cusparseStatus_t
843  cusparseXcsric02_analysis<float>(cusparseHandle_t handle,
844  int m,
845  int nnz,
846  const cusparseMatDescr_t descrA,
847  const float *csrValA,
848  const int *csrRowPtrA,
849  const int *csrColIndA,
850  csric02Info_t info,
851  cusparseSolvePolicy_t policy,
852  void *pBuffer)
853  {
854  return cusparseScsric02_analysis(handle,
855  m,
856  nnz,
857  descrA,
858  csrValA,
859  csrRowPtrA,
860  csrColIndA,
861  info,
862  policy,
863  pBuffer);
864  }
865 
866  template <>
867  cusparseStatus_t
868  cusparseXcsric02_analysis<double>(cusparseHandle_t handle,
869  int m,
870  int nnz,
871  const cusparseMatDescr_t descrA,
872  const double *csrValA,
873  const int *csrRowPtrA,
874  const int *csrColIndA,
875  csric02Info_t info,
876  cusparseSolvePolicy_t policy,
877  void *pBuffer)
878  {
879  return cusparseDcsric02_analysis(handle,
880  m,
881  nnz,
882  descrA,
883  csrValA,
884  csrRowPtrA,
885  csrColIndA,
886  info,
887  policy,
888  pBuffer);
889  }
890 
891  /*
892  template <>
893  cusparseStatus_t
894  cusparseXcsric02_analysis<cuComplex>(cusparseHandle_t handle,
895  int m,
896  int nnz,
897  const cusparseMatDescr_t descrA,
898  const cuComplex * csrValA,
899  const int * csrRowPtrA,
900  const int * csrColIndA,
901  csric02Info_t info,
902  cusparseSolvePolicy_t policy,
903  void * pBuffer)
904  {
905  return cusparseCcsric02_analysis(handle,
906  m,
907  nnz,
908  descrA,
909  csrValA,
910  csrRowPtrA,
911  csrColIndA,
912  info,
913  policy,
914  pBuffer);
915  }
916 
917  template <>
918  cusparseStatus_t
919  cusparseXcsric02_analysis<cuDoubleComplex>(cusparseHandle_t handle,
920  int m,
921  int nnz,
922  const cusparseMatDescr_t descrA,
923  const cuDoubleComplex * csrValA,
924  const int * csrRowPtrA,
925  const int * csrColIndA,
926  csric02Info_t info,
927  cusparseSolvePolicy_t policy,
928  void * pBuffer)
929  {
930  return cusparseZcsric02_analysis(handle,
931  m,
932  nnz,
933  descrA,
934  csrValA,
935  csrRowPtrA,
936  csrColIndA,
937  info,
938  policy,
939  pBuffer);
940  }
941 */
942 
943 
944 
951  template <typename Number>
952  cusparseStatus_t
953  cusparseXcsrsv2_bufferSize(cusparseHandle_t /*handle*/,
954  cusparseOperation_t /*transA*/,
955  int /*m*/,
956  int /*nnz*/,
957  const cusparseMatDescr_t /*descrA*/,
958  Number * /*csrValA*/,
959  const int * /*csrRowPtrA*/,
960  const int * /*csrColIndA*/,
961  csrsv2Info_t /*info*/,
962  int * /*pBufferSizeInBytes*/)
963  {
964  AssertThrow(false, ExcNotImplemented());
965  return CUSPARSE_STATUS_INVALID_VALUE;
966  }
967 
968  template <>
969  cusparseStatus_t
970  cusparseXcsrsv2_bufferSize<float>(cusparseHandle_t handle,
971  cusparseOperation_t transA,
972  int m,
973  int nnz,
974  const cusparseMatDescr_t descrA,
975  float *csrValA,
976  const int *csrRowPtrA,
977  const int *csrColIndA,
978  csrsv2Info_t info,
979  int *pBufferSizeInBytes)
980  {
981  return cusparseScsrsv2_bufferSize(handle,
982  transA,
983  m,
984  nnz,
985  descrA,
986  csrValA,
987  csrRowPtrA,
988  csrColIndA,
989  info,
990  pBufferSizeInBytes);
991  }
992 
993  template <>
994  cusparseStatus_t
995  cusparseXcsrsv2_bufferSize<double>(cusparseHandle_t handle,
996  cusparseOperation_t transA,
997  int m,
998  int nnz,
999  const cusparseMatDescr_t descrA,
1000  double *csrValA,
1001  const int *csrRowPtrA,
1002  const int *csrColIndA,
1003  csrsv2Info_t info,
1004  int *pBufferSizeInBytes)
1005  {
1006  return cusparseDcsrsv2_bufferSize(handle,
1007  transA,
1008  m,
1009  nnz,
1010  descrA,
1011  csrValA,
1012  csrRowPtrA,
1013  csrColIndA,
1014  info,
1015  pBufferSizeInBytes);
1016  }
1017 
1018  /*
1019  template <>
1020  cusparseStatus_t
1021  cusparseXcsrsv2_bufferSize<cuComplex>(cusparseHandle_t handle,
1022  cusparseOperation_t transA,
1023  int m,
1024  int nnz,
1025  const cusparseMatDescr_t descrA,
1026  cuComplex * csrValA,
1027  const int * csrRowPtrA,
1028  const int * csrColIndA,
1029  csrsv2Info_t info,
1030  int *pBufferSizeInBytes)
1031  {
1032  return cusparseCcsrsv2_bufferSize(handle,
1033  transA,
1034  m,
1035  nnz,
1036  descrA,
1037  csrValA,
1038  csrRowPtrA,
1039  csrColIndA,
1040  info,
1041  pBufferSizeInBytes);
1042  }
1043 
1044  template <>
1045  cusparseStatus_t
1046  cusparseXcsrsv2_bufferSize<cuDoubleComplex>(cusparseHandle_t handle,
1047  cusparseOperation_t transA,
1048  int m,
1049  int nnz,
1050  const cusparseMatDescr_t descrA,
1051  cuDoubleComplex * csrValA,
1052  const int * csrRowPtrA,
1053  const int * csrColIndA,
1054  csrsv2Info_t info,
1055  int * pBufferSizeInBytes)
1056  {
1057  return cusparseZcsrsv2_bufferSize(handle,
1058  transA,
1059  m,
1060  nnz,
1061  descrA,
1062  csrValA,
1063  csrRowPtrA,
1064  csrColIndA,
1065  info,
1066  pBufferSizeInBytes);
1067  }
1068 */
1069 
1070 
1071 
1078  template <typename Number>
1079  cusparseStatus_t
1080  cusparseXcsric02_bufferSize(cusparseHandle_t /*handle*/,
1081  int /*m*/,
1082  int /*nnz*/,
1083  const cusparseMatDescr_t /*descrA*/,
1084  Number * /*csrValA*/,
1085  const int * /*csrRowPtrA*/,
1086  const int * /*csrColIndA*/,
1087  csric02Info_t /*info*/,
1088  int * /*pBufferSizeInBytes*/)
1089  {
1090  AssertThrow(false, ExcNotImplemented());
1091  return CUSPARSE_STATUS_INVALID_VALUE;
1092  }
1093 
1094  template <>
1095  cusparseStatus_t
1096  cusparseXcsric02_bufferSize<float>(cusparseHandle_t handle,
1097  int m,
1098  int nnz,
1099  const cusparseMatDescr_t descrA,
1100  float *csrValA,
1101  const int *csrRowPtrA,
1102  const int *csrColIndA,
1103  csric02Info_t info,
1104  int *pBufferSizeInBytes)
1105  {
1106  return cusparseScsric02_bufferSize(handle,
1107  m,
1108  nnz,
1109  descrA,
1110  csrValA,
1111  csrRowPtrA,
1112  csrColIndA,
1113  info,
1114  pBufferSizeInBytes);
1115  }
1116 
1117  template <>
1118  cusparseStatus_t
1119  cusparseXcsric02_bufferSize<double>(cusparseHandle_t handle,
1120  int m,
1121  int nnz,
1122  const cusparseMatDescr_t descrA,
1123  double *csrValA,
1124  const int *csrRowPtrA,
1125  const int *csrColIndA,
1126  csric02Info_t info,
1127  int *pBufferSizeInBytes)
1128  {
1129  return cusparseDcsric02_bufferSize(handle,
1130  m,
1131  nnz,
1132  descrA,
1133  csrValA,
1134  csrRowPtrA,
1135  csrColIndA,
1136  info,
1137  pBufferSizeInBytes);
1138  }
1139 
1140  /*
1141  template <>
1142  cusparseStatus_t
1143  cusparseXcsric02_bufferSize<cuComplex>(cusparseHandle_t handle,
1144  int m,
1145  int nnz,
1146  const cusparseMatDescr_t descrA,
1147  cuComplex * csrValA,
1148  const int * csrRowPtrA,
1149  const int * csrColIndA,
1150  csric02Info_t info,
1151  int *pBufferSizeInBytes)
1152  {
1153  return cusparseCcsric02_bufferSize(handle,
1154  m,
1155  nnz,
1156  descrA,
1157  csrValA,
1158  csrRowPtrA,
1159  csrColIndA,
1160  info,
1161  pBufferSizeInBytes);
1162  }
1163 
1164  template <>
1165  cusparseStatus_t
1166  cusparseXcsric02_bufferSize<cuDoubleComplex>(cusparseHandle_t handle,
1167  int m,
1168  int nnz,
1169  const cusparseMatDescr_t descrA,
1170  cuDoubleComplex * csrValA,
1171  const int * csrRowPtrA,
1172  const int * csrColIndA,
1173  csric02Info_t info,
1174  int * pBufferSizeInBytes)
1175  {
1176  return cusparseZcsric02_bufferSize(handle,
1177  m,
1178  nnz,
1179  descrA,
1180  csrValA,
1181  csrRowPtrA,
1182  csrColIndA,
1183  info,
1184  pBufferSizeInBytes);
1185  }
1186  */
1187 } // namespace
1188 
1189 namespace CUDAWrappers
1190 {
1191  template <typename Number>
1193  bool use_level_analysis_)
1194  : use_level_analysis(use_level_analysis_)
1195  {}
1196 
1197 
1198 
1199  template <typename Number>
1202  , P_val_dev(nullptr, Utilities::CUDA::delete_device_data<Number>)
1203  , P_row_ptr_dev(nullptr)
1204  , P_column_index_dev(nullptr)
1205  , tmp_dev(nullptr, Utilities::CUDA::delete_device_data<Number>)
1206  , buffer_dev(nullptr, Utilities::CUDA::delete_device_data<void>)
1207  , policy_L(CUSPARSE_SOLVE_POLICY_USE_LEVEL)
1208  , policy_Lt(CUSPARSE_SOLVE_POLICY_USE_LEVEL)
1209  , policy_M(CUSPARSE_SOLVE_POLICY_USE_LEVEL)
1210  , n_rows(0)
1211  , n_nonzero_elements(0)
1212  {
1213  // step 1: create a descriptor which contains
1214  // - matrix M is base-0
1215  // - matrix L is base-0
1216  // - matrix L is lower triangular
1217  // - matrix L has non-unit diagonal
1218  cusparseStatus_t status = cusparseCreateMatDescr(&descr_M);
1219  AssertCusparse(status);
1220  status = cusparseSetMatIndexBase(descr_M, CUSPARSE_INDEX_BASE_ZERO);
1221  AssertCusparse(status);
1222  status = cusparseSetMatType(descr_M, CUSPARSE_MATRIX_TYPE_GENERAL);
1223  AssertCusparse(status);
1224 
1225  status = cusparseCreateMatDescr(&descr_L);
1226  AssertCusparse(status);
1227  status = cusparseSetMatIndexBase(descr_L, CUSPARSE_INDEX_BASE_ZERO);
1228  AssertCusparse(status);
1229  status = cusparseSetMatType(descr_L, CUSPARSE_MATRIX_TYPE_GENERAL);
1230  AssertCusparse(status);
1231  status = cusparseSetMatFillMode(descr_L, CUSPARSE_FILL_MODE_LOWER);
1232  AssertCusparse(status);
1233  status = cusparseSetMatDiagType(descr_L, CUSPARSE_DIAG_TYPE_NON_UNIT);
1234  AssertCusparse(status);
1235 
1236  // step 2: create a empty info structure
1237  // we need one info for csric02 and two info's for csrsv2
1238  status = cusparseCreateCsric02Info(&info_M);
1239  AssertCusparse(status);
1240  status = cusparseCreateCsrsv2Info(&info_L);
1241  AssertCusparse(status);
1242  status = cusparseCreateCsrsv2Info(&info_Lt);
1243  AssertCusparse(status);
1244  }
1245 
1246 
1247 
1248  template <typename Number>
1250  {
1251  // step 8: free resources
1252  cusparseStatus_t status = cusparseDestroyMatDescr(descr_M);
1253  AssertNothrowCusparse(status);
1254 
1255  status = cusparseDestroyMatDescr(descr_L);
1256  AssertNothrowCusparse(status);
1257 
1258  status = cusparseDestroyCsric02Info(info_M);
1259  AssertNothrowCusparse(status);
1260 
1261  status = cusparseDestroyCsrsv2Info(info_L);
1262  AssertNothrowCusparse(status);
1263 
1264  status = cusparseDestroyCsrsv2Info(info_Lt);
1265  AssertNothrowCusparse(status);
1266  }
1267 
1268 
1269 
1270  template <typename Number>
1271  void
1273  const AdditionalData &additional_data)
1274  {
1275  if (additional_data.use_level_analysis)
1276  {
1277  policy_L = CUSPARSE_SOLVE_POLICY_USE_LEVEL;
1278  policy_Lt = CUSPARSE_SOLVE_POLICY_USE_LEVEL;
1279  policy_M = CUSPARSE_SOLVE_POLICY_USE_LEVEL;
1280  }
1281  else
1282  {
1283  policy_L = CUSPARSE_SOLVE_POLICY_NO_LEVEL;
1284  policy_Lt = CUSPARSE_SOLVE_POLICY_NO_LEVEL;
1285  policy_M = CUSPARSE_SOLVE_POLICY_NO_LEVEL;
1286  }
1287 
1288  n_rows = A.m();
1289  n_nonzero_elements = A.n_nonzero_elements();
1290  AssertDimension(A.m(), A.n());
1291 
1292  matrix_pointer = &A;
1293  const Number *A_val_dev;
1294  std::tie(A_val_dev,
1295  P_column_index_dev,
1296  P_row_ptr_dev,
1297  std::ignore,
1298  std::ignore) = A.get_cusparse_matrix();
1299 
1300  // create a copy of the matrix entries since the algorithm works in-place.
1301  P_val_dev.reset(
1302  Utilities::CUDA::allocate_device_data<Number>(n_nonzero_elements));
1303  cudaError_t cuda_status = cudaMemcpy(P_val_dev.get(),
1304  A_val_dev,
1305  n_nonzero_elements * sizeof(Number),
1306  cudaMemcpyDeviceToDevice);
1307  AssertCuda(cuda_status);
1308 
1309  // initialize an internal buffer we need later on
1310  tmp_dev.reset(Utilities::CUDA::allocate_device_data<Number>(n_rows));
1311 
1312  // step 3: query how much memory used in csric02 and csrsv2, and allocate
1313  // the buffer
1314  int BufferSize_M;
1315  cusparseStatus_t status = cusparseXcsric02_bufferSize(cusparse_handle,
1316  n_rows,
1317  n_nonzero_elements,
1318  descr_M,
1319  P_val_dev.get(),
1320  P_row_ptr_dev,
1321  P_column_index_dev,
1322  info_M,
1323  &BufferSize_M);
1324  AssertCusparse(status);
1325 
1326  int BufferSize_L;
1327  status = cusparseXcsrsv2_bufferSize(cusparse_handle,
1328  CUSPARSE_OPERATION_NON_TRANSPOSE,
1329  n_rows,
1330  n_nonzero_elements,
1331  descr_L,
1332  P_val_dev.get(),
1333  P_row_ptr_dev,
1334  P_column_index_dev,
1335  info_L,
1336  &BufferSize_L);
1337  AssertCusparse(status);
1338 
1339  int BufferSize_Lt;
1340  status = cusparseXcsrsv2_bufferSize(cusparse_handle,
1341  CUSPARSE_OPERATION_TRANSPOSE,
1342  n_rows,
1343  n_nonzero_elements,
1344  descr_L,
1345  P_val_dev.get(),
1346  P_row_ptr_dev,
1347  P_column_index_dev,
1348  info_Lt,
1349  &BufferSize_Lt);
1350  AssertCusparse(status);
1351 
1352  const int BufferSize =
1353  std::max(BufferSize_M, std::max(BufferSize_L, BufferSize_Lt));
1354  // workaround: since allocate_device_data needs a type, we pass char
1355  // which is required to have size 1.
1356  buffer_dev.reset(static_cast<void *>(
1357  Utilities::CUDA::allocate_device_data<char>(BufferSize / sizeof(char))));
1358 
1359  // step 4: perform analysis of incomplete Cholesky on M
1360  // perform analysis of triangular solve on L
1361  // perform analysis of triangular solve on L'
1362  // The lower triangular part of M has the same sparsity pattern as L, so
1363  // we can do analysis of csric02 and csrsv2 simultaneously.
1364 
1365  status = cusparseXcsric02_analysis(cusparse_handle,
1366  n_rows,
1367  n_nonzero_elements,
1368  descr_M,
1369  P_val_dev.get(),
1370  P_row_ptr_dev,
1371  P_column_index_dev,
1372  info_M,
1373  policy_M,
1374  buffer_dev.get());
1375  AssertCusparse(status);
1376 
1377  int structural_zero;
1378  status =
1379  cusparseXcsric02_zeroPivot(cusparse_handle, info_M, &structural_zero);
1380  AssertCusparse(status);
1381 
1382  status = cusparseXcsrsv2_analysis(cusparse_handle,
1383  CUSPARSE_OPERATION_TRANSPOSE,
1384  n_rows,
1385  n_nonzero_elements,
1386  descr_L,
1387  P_val_dev.get(),
1388  P_row_ptr_dev,
1389  P_column_index_dev,
1390  info_Lt,
1391  policy_Lt,
1392  buffer_dev.get());
1393  AssertCusparse(status);
1394 
1395  status = cusparseXcsrsv2_analysis(cusparse_handle,
1396  CUSPARSE_OPERATION_NON_TRANSPOSE,
1397  n_rows,
1398  n_nonzero_elements,
1399  descr_L,
1400  P_val_dev.get(),
1401  P_row_ptr_dev,
1402  P_column_index_dev,
1403  info_L,
1404  policy_L,
1405  buffer_dev.get());
1406  AssertCusparse(status);
1407 
1408  // step 5: M = L * L'
1409  status = cusparseXcsric02(cusparse_handle,
1410  n_rows,
1411  n_nonzero_elements,
1412  descr_M,
1413  P_val_dev.get(),
1414  P_row_ptr_dev,
1415  P_column_index_dev,
1416  info_M,
1417  policy_M,
1418  buffer_dev.get());
1419  AssertCusparse(status);
1420 
1421  int numerical_zero;
1422  status =
1423  cusparseXcsric02_zeroPivot(cusparse_handle, info_M, &numerical_zero);
1424  AssertCusparse(status);
1425  }
1426 
1427 
1428 
1429  template <typename Number>
1430  void
1434  {
1435  Assert(P_val_dev != nullptr, ExcNotInitialized());
1436  Assert(P_row_ptr_dev != nullptr, ExcNotInitialized());
1437  Assert(P_column_index_dev != nullptr, ExcNotInitialized());
1438  AssertDimension(dst.size(), static_cast<unsigned int>(n_rows));
1439  AssertDimension(src.size(), static_cast<unsigned int>(n_rows));
1440  Assert(tmp_dev != nullptr, ExcInternalError());
1441 
1442  const Number *const src_dev = src.get_values();
1443  Number *const dst_dev = dst.get_values();
1444  // step 6: solve L*z = alpha*x
1445  const Number alpha = internal::NumberType<Number>::value(1.);
1446  cusparseStatus_t status =
1447  cusparseXcsrsv2_solve(cusparse_handle,
1448  CUSPARSE_OPERATION_NON_TRANSPOSE,
1449  n_rows,
1450  n_nonzero_elements,
1451  &alpha,
1452  descr_L,
1453  P_val_dev.get(),
1454  P_row_ptr_dev,
1455  P_column_index_dev,
1456  info_L,
1457  src_dev,
1458  tmp_dev.get(),
1459  policy_L,
1460  buffer_dev.get());
1461  AssertCusparse(status);
1462 
1463  // step 7: solve L'*y = alpha*z
1464  status = cusparseXcsrsv2_solve(cusparse_handle,
1465  CUSPARSE_OPERATION_TRANSPOSE,
1466  n_rows,
1467  n_nonzero_elements,
1468  &alpha,
1469  descr_L,
1470  P_val_dev.get(),
1471  P_row_ptr_dev,
1472  P_column_index_dev,
1473  info_Lt,
1474  tmp_dev.get(),
1475  dst_dev,
1476  policy_Lt,
1477  buffer_dev.get());
1478  AssertCusparse(status);
1479  }
1480 
1481 
1482 
1483  template <typename Number>
1484  void
1488  {
1489  // the constructed preconditioner is symmetric
1490  vmult(dst, src);
1491  }
1492 
1493 
1494 
1495  template <typename Number>
1497  bool use_level_analysis_)
1498  : use_level_analysis(use_level_analysis_)
1499  {}
1500 
1501 
1502 
1503  template <typename Number>
1505  const Utilities::CUDA::Handle &handle)
1507  , P_val_dev(nullptr, Utilities::CUDA::delete_device_data<Number>)
1508  , P_row_ptr_dev(nullptr)
1509  , P_column_index_dev(nullptr)
1510  , tmp_dev(nullptr, Utilities::CUDA::delete_device_data<Number>)
1511  , buffer_dev(nullptr, Utilities::CUDA::delete_device_data<void>)
1512  , policy_L(CUSPARSE_SOLVE_POLICY_USE_LEVEL)
1513  , policy_U(CUSPARSE_SOLVE_POLICY_USE_LEVEL)
1514  , policy_M(CUSPARSE_SOLVE_POLICY_USE_LEVEL)
1515  , n_rows(0)
1516  , n_nonzero_elements(0)
1517  {
1518  // step 1: create a descriptor which contains
1519  // - matrix M is base-0
1520  // - matrix L is base-0
1521  // - matrix L is lower triangular
1522  // - matrix L has unit diagonal
1523  // - matrix U is base-0
1524  // - matrix U is upper triangular
1525  // - matrix U has non-unit diagonal
1526  cusparseStatus_t status = cusparseCreateMatDescr(&descr_M);
1527  AssertCusparse(status);
1528  status = cusparseSetMatIndexBase(descr_M, CUSPARSE_INDEX_BASE_ZERO);
1529  AssertCusparse(status);
1530  status = cusparseSetMatType(descr_M, CUSPARSE_MATRIX_TYPE_GENERAL);
1531  AssertCusparse(status);
1532 
1533  status = cusparseCreateMatDescr(&descr_L);
1534  AssertCusparse(status);
1535  status = cusparseSetMatIndexBase(descr_L, CUSPARSE_INDEX_BASE_ZERO);
1536  AssertCusparse(status);
1537  status = cusparseSetMatType(descr_L, CUSPARSE_MATRIX_TYPE_GENERAL);
1538  AssertCusparse(status);
1539  status = cusparseSetMatFillMode(descr_L, CUSPARSE_FILL_MODE_LOWER);
1540  AssertCusparse(status);
1541  status = cusparseSetMatDiagType(descr_L, CUSPARSE_DIAG_TYPE_UNIT);
1542  AssertCusparse(status);
1543 
1544  status = cusparseCreateMatDescr(&descr_U);
1545  AssertCusparse(status);
1546  status = cusparseSetMatIndexBase(descr_U, CUSPARSE_INDEX_BASE_ZERO);
1547  AssertCusparse(status);
1548  status = cusparseSetMatType(descr_U, CUSPARSE_MATRIX_TYPE_GENERAL);
1549  AssertCusparse(status);
1550  status = cusparseSetMatFillMode(descr_U, CUSPARSE_FILL_MODE_UPPER);
1551  AssertCusparse(status);
1552  status = cusparseSetMatDiagType(descr_U, CUSPARSE_DIAG_TYPE_NON_UNIT);
1553  AssertCusparse(status);
1554 
1555  // step 2: create a empty info structure
1556  // we need one info for csrilu02 and two info's for csrsv2
1557  status = cusparseCreateCsrilu02Info(&info_M);
1558  AssertCusparse(status);
1559  status = cusparseCreateCsrsv2Info(&info_L);
1560  AssertCusparse(status);
1561  status = cusparseCreateCsrsv2Info(&info_U);
1562  AssertCusparse(status);
1563  }
1564 
1565 
1566 
1567  template <typename Number>
1569  {
1570  // step 8: free resources
1571  cusparseStatus_t status = cusparseDestroyMatDescr(descr_M);
1572  AssertNothrowCusparse(status);
1573 
1574  status = cusparseDestroyMatDescr(descr_L);
1575  AssertNothrowCusparse(status);
1576 
1577  status = cusparseDestroyMatDescr(descr_U);
1578  AssertNothrowCusparse(status);
1579 
1580  status = cusparseDestroyCsrilu02Info(info_M);
1581  AssertNothrowCusparse(status);
1582 
1583  status = cusparseDestroyCsrsv2Info(info_L);
1584  AssertNothrowCusparse(status);
1585 
1586  status = cusparseDestroyCsrsv2Info(info_U);
1587  AssertNothrowCusparse(status);
1588  }
1589 
1590 
1591 
1592  template <typename Number>
1593  void
1595  const AdditionalData &additional_data)
1596  {
1597  if (additional_data.use_level_analysis)
1598  {
1599  policy_L = CUSPARSE_SOLVE_POLICY_USE_LEVEL;
1600  policy_U = CUSPARSE_SOLVE_POLICY_USE_LEVEL;
1601  policy_M = CUSPARSE_SOLVE_POLICY_USE_LEVEL;
1602  }
1603  else
1604  {
1605  policy_L = CUSPARSE_SOLVE_POLICY_NO_LEVEL;
1606  policy_U = CUSPARSE_SOLVE_POLICY_NO_LEVEL;
1607  policy_M = CUSPARSE_SOLVE_POLICY_NO_LEVEL;
1608  }
1609 
1610  matrix_pointer = &A;
1611  n_rows = A.m();
1612  n_nonzero_elements = A.n_nonzero_elements();
1613  AssertDimension(A.m(), A.n());
1614 
1615  const Number *A_val_dev;
1616  std::tie(A_val_dev,
1617  P_column_index_dev,
1618  P_row_ptr_dev,
1619  std::ignore,
1620  std::ignore) = A.get_cusparse_matrix();
1621 
1622  // create a copy of the matrix entries since the algorithm works in-place.
1623  P_val_dev.reset(
1624  Utilities::CUDA::allocate_device_data<Number>(n_nonzero_elements));
1625  cudaError_t cuda_status = cudaMemcpy(P_val_dev.get(),
1626  A_val_dev,
1627  n_nonzero_elements * sizeof(Number),
1628  cudaMemcpyDeviceToDevice);
1629  AssertCuda(cuda_status);
1630 
1631  // initialize an internal buffer we need later on
1632  tmp_dev.reset(Utilities::CUDA::allocate_device_data<Number>(n_rows));
1633 
1634  // step 3: query how much memory used in csrilu02 and csrsv2, and allocate
1635  // the buffer
1636  int BufferSize_M;
1637  cusparseStatus_t status = cusparseXcsrilu02_bufferSize(cusparse_handle,
1638  n_rows,
1639  n_nonzero_elements,
1640  descr_M,
1641  P_val_dev.get(),
1642  P_row_ptr_dev,
1643  P_column_index_dev,
1644  info_M,
1645  &BufferSize_M);
1646  AssertCusparse(status);
1647 
1648  int BufferSize_L;
1649  status = cusparseXcsrsv2_bufferSize(cusparse_handle,
1650  CUSPARSE_OPERATION_NON_TRANSPOSE,
1651  n_rows,
1652  n_nonzero_elements,
1653  descr_L,
1654  P_val_dev.get(),
1655  P_row_ptr_dev,
1656  P_column_index_dev,
1657  info_L,
1658  &BufferSize_L);
1659  AssertCusparse(status);
1660 
1661  int BufferSize_U;
1662  status = cusparseXcsrsv2_bufferSize(cusparse_handle,
1663  CUSPARSE_OPERATION_NON_TRANSPOSE,
1664  n_rows,
1665  n_nonzero_elements,
1666  descr_U,
1667  P_val_dev.get(),
1668  P_row_ptr_dev,
1669  P_column_index_dev,
1670  info_U,
1671  &BufferSize_U);
1672  AssertCusparse(status);
1673 
1674  const int BufferSize =
1675  std::max(BufferSize_M, std::max(BufferSize_L, BufferSize_U));
1676  // workaround: since allocate_device_data needs a type, we pass char
1677  // which is required to have size 1.
1678  buffer_dev.reset(static_cast<void *>(
1679  Utilities::CUDA::allocate_device_data<char>(BufferSize / sizeof(char))));
1680 
1681  // step 4: perform analysis of incomplete Cholesky on M
1682  // perform analysis of triangular solve on L
1683  // perform analysis of triangular solve on U
1684  // The lower(upper) triangular part of M has the same sparsity pattern as
1685  // L(U), we can do analysis of csrilu0 and csrsv2 simultaneously.
1686 
1687  status = cusparseXcsrilu02_analysis(cusparse_handle,
1688  n_rows,
1689  n_nonzero_elements,
1690  descr_M,
1691  P_val_dev.get(),
1692  P_row_ptr_dev,
1693  P_column_index_dev,
1694  info_M,
1695  policy_M,
1696  buffer_dev.get());
1697  AssertCusparse(status);
1698 
1699  int structural_zero;
1700  status =
1701  cusparseXcsrilu02_zeroPivot(cusparse_handle, info_M, &structural_zero);
1702  AssertCusparse(status);
1703 
1704  status = cusparseXcsrsv2_analysis(cusparse_handle,
1705  CUSPARSE_OPERATION_NON_TRANSPOSE,
1706  n_rows,
1707  n_nonzero_elements,
1708  descr_L,
1709  P_val_dev.get(),
1710  P_row_ptr_dev,
1711  P_column_index_dev,
1712  info_L,
1713  policy_L,
1714  buffer_dev.get());
1715  AssertCusparse(status);
1716 
1717  status = cusparseXcsrsv2_analysis(cusparse_handle,
1718  CUSPARSE_OPERATION_NON_TRANSPOSE,
1719  n_rows,
1720  n_nonzero_elements,
1721  descr_U,
1722  P_val_dev.get(),
1723  P_row_ptr_dev,
1724  P_column_index_dev,
1725  info_U,
1726  policy_U,
1727  buffer_dev.get());
1728 
1729  // step 5: M = L * U
1730  status = cusparseXcsrilu02(cusparse_handle,
1731  n_rows,
1732  n_nonzero_elements,
1733  descr_M,
1734  P_val_dev.get(),
1735  P_row_ptr_dev,
1736  P_column_index_dev,
1737  info_M,
1738  policy_M,
1739  buffer_dev.get());
1740  AssertCusparse(status);
1741 
1742  int numerical_zero;
1743  status =
1744  cusparseXcsrilu02_zeroPivot(cusparse_handle, info_M, &numerical_zero);
1745  AssertCusparse(status);
1746  }
1747 
1748 
1749 
1750  template <typename Number>
1751  void
1755  {
1756  Assert(P_val_dev != nullptr, ExcNotInitialized());
1757  Assert(P_row_ptr_dev != nullptr, ExcNotInitialized());
1758  Assert(P_column_index_dev != nullptr, ExcNotInitialized());
1759  AssertDimension(dst.size(), static_cast<unsigned int>(n_rows));
1760  AssertDimension(src.size(), static_cast<unsigned int>(n_rows));
1761  Assert(tmp_dev != nullptr, ExcInternalError());
1762 
1763  const Number *const src_dev = src.get_values();
1764  Number *const dst_dev = dst.get_values();
1765 
1766  // step 6: solve L*z = alpha*x
1767  const Number alpha = internal::NumberType<Number>::value(1.);
1768  cusparseStatus_t status =
1769  cusparseXcsrsv2_solve(cusparse_handle,
1770  CUSPARSE_OPERATION_NON_TRANSPOSE,
1771  n_rows,
1772  n_nonzero_elements,
1773  &alpha,
1774  descr_L,
1775  P_val_dev.get(),
1776  P_row_ptr_dev,
1777  P_column_index_dev,
1778  info_L,
1779  src_dev,
1780  tmp_dev.get(),
1781  policy_L,
1782  buffer_dev.get());
1783  AssertCusparse(status);
1784 
1785  // step 7: solve U*y = alpha*z
1786  status = cusparseXcsrsv2_solve(cusparse_handle,
1787  CUSPARSE_OPERATION_NON_TRANSPOSE,
1788  n_rows,
1789  n_nonzero_elements,
1790  &alpha,
1791  descr_U,
1792  P_val_dev.get(),
1793  P_row_ptr_dev,
1794  P_column_index_dev,
1795  info_U,
1796  tmp_dev.get(),
1797  dst_dev,
1798  policy_U,
1799  buffer_dev.get());
1800  AssertCusparse(status);
1801  }
1802 
1803 
1804 
1805  template <typename Number>
1806  void
1809  const LinearAlgebra::CUDAWrappers::Vector<Number> & /*src*/) const
1810  {
1811  Assert(false, ExcNotImplemented());
1812  }
1813 
1814 
1815 
1816  // explicit instantiations
1817  template class PreconditionIC<float>;
1818  template class PreconditionIC<double>;
1819  // template class PreconditionIC<cuComplex>;
1820  // template class PreconditionIC<cuDoubleComplex>;
1821  template class PreconditionILU<float>;
1822  template class PreconditionILU<double>;
1823  // template class PreconditionILU<cuComplex>;
1824  // template class PreconditionILU<cuDoubleComplex>;
1825 } // namespace CUDAWrappers
1826 
cusparseSolvePolicy_t policy_M
std::unique_ptr< Number[], void(*)(Number *)> P_val_dev
void Tvmult(LinearAlgebra::CUDAWrappers::Vector< Number > &dst, const LinearAlgebra::CUDAWrappers::Vector< Number > &src) const
std::unique_ptr< Number[], void(*)(Number *)> tmp_dev
void initialize(const SparseMatrix< Number > &matrix, const AdditionalData &additional_data=AdditionalData())
cusparseSolvePolicy_t policy_Lt
cusparseSolvePolicy_t policy_L
void vmult(LinearAlgebra::CUDAWrappers::Vector< Number > &dst, const LinearAlgebra::CUDAWrappers::Vector< Number > &src) const
PreconditionIC(const Utilities::CUDA::Handle &handle)
std::unique_ptr< void, void(*)(void *)> buffer_dev
cusparseSolvePolicy_t policy_M
void vmult(LinearAlgebra::CUDAWrappers::Vector< Number > &dst, const LinearAlgebra::CUDAWrappers::Vector< Number > &src) const
std::unique_ptr< void, void(*)(void *)> buffer_dev
void initialize(const SparseMatrix< Number > &matrix, const AdditionalData &additional_data=AdditionalData())
cusparseSolvePolicy_t policy_U
PreconditionILU(const Utilities::CUDA::Handle &handle)
std::unique_ptr< Number[], void(*)(Number *)> P_val_dev
std::unique_ptr< Number[], void(*)(Number *)> tmp_dev
void Tvmult(LinearAlgebra::CUDAWrappers::Vector< Number > &dst, const LinearAlgebra::CUDAWrappers::Vector< Number > &src) const
cusparseSolvePolicy_t policy_L
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
#define AssertCusparse(error_code)
Definition: exceptions.h:2027
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
#define AssertNothrowCusparse(error_code)
Definition: exceptions.h:2057
#define AssertCuda(error_code)
Definition: exceptions.h:1942
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
static const char A
Default CUDA
Definition: memory_space.h:53
void delete_device_data(Number *device_ptr) noexcept
Definition: cuda.h:121
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T & value(const T &t)
Definition: numbers.h:703