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