1 #define PETSC_SKIP_IMMINTRIN_H_CUDAWORKAROUND 1 2 3 #include <petscconf.h> 4 #include <../src/mat/impls/aij/mpi/mpiaij.h> /*I "petscmat.h" I*/ 5 #include <../src/mat/impls/aij/seq/seqviennacl/viennaclmatimpl.h> 6 7 PetscErrorCode MatMPIAIJSetPreallocation_MPIAIJViennaCL(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[]) { 8 Mat_MPIAIJ *b = (Mat_MPIAIJ *)B->data; 9 10 PetscFunctionBegin; 11 PetscCall(PetscLayoutSetUp(B->rmap)); 12 PetscCall(PetscLayoutSetUp(B->cmap)); 13 if (!B->preallocated) { 14 /* Explicitly create the two MATSEQAIJVIENNACL matrices. */ 15 PetscCall(MatCreate(PETSC_COMM_SELF, &b->A)); 16 PetscCall(MatSetSizes(b->A, B->rmap->n, B->cmap->n, B->rmap->n, B->cmap->n)); 17 PetscCall(MatSetType(b->A, MATSEQAIJVIENNACL)); 18 PetscCall(PetscLogObjectParent((PetscObject)B, (PetscObject)b->A)); 19 PetscCall(MatCreate(PETSC_COMM_SELF, &b->B)); 20 PetscCall(MatSetSizes(b->B, B->rmap->n, B->cmap->N, B->rmap->n, B->cmap->N)); 21 PetscCall(MatSetType(b->B, MATSEQAIJVIENNACL)); 22 PetscCall(PetscLogObjectParent((PetscObject)B, (PetscObject)b->B)); 23 } 24 PetscCall(MatSeqAIJSetPreallocation(b->A, d_nz, d_nnz)); 25 PetscCall(MatSeqAIJSetPreallocation(b->B, o_nz, o_nnz)); 26 B->preallocated = PETSC_TRUE; 27 PetscFunctionReturn(0); 28 } 29 30 PetscErrorCode MatAssemblyEnd_MPIAIJViennaCL(Mat A, MatAssemblyType mode) { 31 Mat_MPIAIJ *b = (Mat_MPIAIJ *)A->data; 32 PetscBool v; 33 34 PetscFunctionBegin; 35 PetscCall(MatAssemblyEnd_MPIAIJ(A, mode)); 36 PetscCall(PetscObjectTypeCompare((PetscObject)b->lvec, VECSEQVIENNACL, &v)); 37 if (!v) { 38 PetscInt m; 39 PetscCall(VecGetSize(b->lvec, &m)); 40 PetscCall(VecDestroy(&b->lvec)); 41 PetscCall(VecCreateSeqViennaCL(PETSC_COMM_SELF, m, &b->lvec)); 42 } 43 PetscFunctionReturn(0); 44 } 45 46 PetscErrorCode MatDestroy_MPIAIJViennaCL(Mat A) { 47 PetscFunctionBegin; 48 PetscCall(MatDestroy_MPIAIJ(A)); 49 PetscFunctionReturn(0); 50 } 51 52 PETSC_EXTERN PetscErrorCode MatCreate_MPIAIJViennaCL(Mat A) { 53 PetscFunctionBegin; 54 PetscCall(MatCreate_MPIAIJ(A)); 55 A->boundtocpu = PETSC_FALSE; 56 PetscCall(PetscFree(A->defaultvectype)); 57 PetscCall(PetscStrallocpy(VECVIENNACL, &A->defaultvectype)); 58 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMPIAIJSetPreallocation_C", MatMPIAIJSetPreallocation_MPIAIJViennaCL)); 59 A->ops->assemblyend = MatAssemblyEnd_MPIAIJViennaCL; 60 PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATMPIAIJVIENNACL)); 61 PetscFunctionReturn(0); 62 } 63 64 /*@C 65 MatCreateAIJViennaCL - Creates a sparse matrix in AIJ (compressed row) format 66 (the default parallel PETSc format). This matrix will ultimately be pushed down 67 to GPUs and use the ViennaCL library for calculations. For good matrix 68 assembly performance the user should preallocate the matrix storage by setting 69 the parameter nz (or the array nnz). By setting these parameters accurately, 70 performance during matrix assembly can be increased substantially. 71 72 Collective 73 74 Input Parameters: 75 + comm - MPI communicator, set to PETSC_COMM_SELF 76 . m - number of rows 77 . n - number of columns 78 . nz - number of nonzeros per row (same for all rows) 79 - nnz - array containing the number of nonzeros in the various rows 80 (possibly different for each row) or NULL 81 82 Output Parameter: 83 . A - the matrix 84 85 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 86 MatXXXXSetPreallocation() paradigm instead of this routine directly. 87 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 88 89 Notes: 90 If nnz is given then nz is ignored 91 92 The AIJ format (also called the Yale sparse matrix format or 93 compressed row storage), is fully compatible with standard Fortran 77 94 storage. That is, the stored row and column indices can begin at 95 either one (as in Fortran) or zero. See the users' manual for details. 96 97 Specify the preallocated storage with either nz or nnz (not both). 98 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 99 allocation. For large problems you MUST preallocate memory or you 100 will get TERRIBLE performance, see the users' manual chapter on matrices. 101 102 Level: intermediate 103 104 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatCreateAIJCUSPARSE()`, `MatSetValues()`, `MatSeqAIJSetColumnIndices()`, `MatCreateSeqAIJWithArrays()`, `MatCreateAIJ()`, `MATMPIAIJVIENNACL`, `MATAIJVIENNACL` 105 @*/ 106 PetscErrorCode MatCreateAIJViennaCL(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[], Mat *A) { 107 PetscMPIInt size; 108 109 PetscFunctionBegin; 110 PetscCall(MatCreate(comm, A)); 111 PetscCall(MatSetSizes(*A, m, n, M, N)); 112 PetscCallMPI(MPI_Comm_size(comm, &size)); 113 if (size > 1) { 114 PetscCall(MatSetType(*A, MATMPIAIJVIENNACL)); 115 PetscCall(MatMPIAIJSetPreallocation(*A, d_nz, d_nnz, o_nz, o_nnz)); 116 } else { 117 PetscCall(MatSetType(*A, MATSEQAIJVIENNACL)); 118 PetscCall(MatSeqAIJSetPreallocation(*A, d_nz, d_nnz)); 119 } 120 PetscFunctionReturn(0); 121 } 122 123 /*MC 124 MATAIJVIENNACL - MATMPIAIJVIENNACL= "aijviennacl" = "mpiaijviennacl" - A matrix type to be used for sparse matrices. 125 126 A matrix type (CSR format) whose data resides on GPUs. 127 All matrix calculations are performed using the ViennaCL library. 128 129 This matrix type is identical to MATSEQAIJVIENNACL when constructed with a single process communicator, 130 and MATMPIAIJVIENNACL otherwise. As a result, for single process communicators, 131 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported 132 for communicators controlling multiple processes. It is recommended that you call both of 133 the above preallocation routines for simplicity. 134 135 Options Database Keys: 136 . -mat_type mpiaijviennacl - sets the matrix type to "mpiaijviennacl" during a call to MatSetFromOptions() 137 138 Level: beginner 139 140 .seealso: `MatCreateAIJViennaCL()`, `MATSEQAIJVIENNACL`, `MatCreateSeqAIJVIENNACL()` 141 M*/ 142