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