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 55 PetscErrorCode MatDestroy_MPIAIJViennaCL(Mat A) 56 { 57 PetscErrorCode ierr; 58 59 PetscFunctionBegin; 60 ierr = MatDestroy_MPIAIJ(A);CHKERRQ(ierr); 61 PetscFunctionReturn(0); 62 } 63 64 PETSC_EXTERN PetscErrorCode MatCreate_MPIAIJViennaCL(Mat A) 65 { 66 PetscErrorCode ierr; 67 68 PetscFunctionBegin; 69 ierr = MatCreate_MPIAIJ(A);CHKERRQ(ierr); 70 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",MatMPIAIJSetPreallocation_MPIAIJViennaCL);CHKERRQ(ierr); 71 A->ops->getvecs = MatCreateVecs_MPIAIJViennaCL; 72 73 ierr = PetscObjectChangeTypeName((PetscObject)A,MATMPIAIJVIENNACL);CHKERRQ(ierr); 74 PetscFunctionReturn(0); 75 } 76 77 78 /*@ 79 MatCreateAIJViennaCL - Creates a sparse matrix in AIJ (compressed row) format 80 (the default parallel PETSc format). This matrix will ultimately be pushed down 81 to GPUs and use the ViennaCL library for calculations. For good matrix 82 assembly performance the user should preallocate the matrix storage by setting 83 the parameter nz (or the array nnz). By setting these parameters accurately, 84 performance during matrix assembly can be increased substantially. 85 86 87 Collective on MPI_Comm 88 89 Input Parameters: 90 + comm - MPI communicator, set to PETSC_COMM_SELF 91 . m - number of rows 92 . n - number of columns 93 . nz - number of nonzeros per row (same for all rows) 94 - nnz - array containing the number of nonzeros in the various rows 95 (possibly different for each row) or NULL 96 97 Output Parameter: 98 . A - the matrix 99 100 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 101 MatXXXXSetPreallocation() paradigm instead of this routine directly. 102 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 103 104 Notes: 105 If nnz is given then nz is ignored 106 107 The AIJ format (also called the Yale sparse matrix format or 108 compressed row storage), is fully compatible with standard Fortran 77 109 storage. That is, the stored row and column indices can begin at 110 either one (as in Fortran) or zero. See the users' manual for details. 111 112 Specify the preallocated storage with either nz or nnz (not both). 113 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 114 allocation. For large problems you MUST preallocate memory or you 115 will get TERRIBLE performance, see the users' manual chapter on matrices. 116 117 Level: intermediate 118 119 .seealso: MatCreate(), MatCreateAIJ(), MatCreateAIJCUSP(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatCreateAIJ(), MATMPIAIJVIENNACL, MATAIJVIENNACL 120 @*/ 121 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) 122 { 123 PetscErrorCode ierr; 124 PetscMPIInt size; 125 126 PetscFunctionBegin; 127 ierr = MatCreate(comm,A);CHKERRQ(ierr); 128 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 129 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 130 if (size > 1) { 131 ierr = MatSetType(*A,MATMPIAIJVIENNACL);CHKERRQ(ierr); 132 ierr = MatMPIAIJSetPreallocation(*A,d_nz,d_nnz,o_nz,o_nnz);CHKERRQ(ierr); 133 } else { 134 ierr = MatSetType(*A,MATSEQAIJVIENNACL);CHKERRQ(ierr); 135 ierr = MatSeqAIJSetPreallocation(*A,d_nz,d_nnz);CHKERRQ(ierr); 136 } 137 PetscFunctionReturn(0); 138 } 139 140 /*M 141 MATAIJVIENNACL - MATMPIAIJVIENNACL= "aijviennacl" = "mpiaijviennacl" - A matrix type to be used for sparse matrices. 142 143 A matrix type (CSR format) whose data resides on GPUs. 144 All matrix calculations are performed using the ViennaCL library. 145 146 This matrix type is identical to MATSEQAIJVIENNACL when constructed with a single process communicator, 147 and MATMPIAIJVIENNACL otherwise. As a result, for single process communicators, 148 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported 149 for communicators controlling multiple processes. It is recommended that you call both of 150 the above preallocation routines for simplicity. 151 152 Options Database Keys: 153 + -mat_type mpiaijviennacl - sets the matrix type to "mpiaijviennacl" during a call to MatSetFromOptions() 154 155 Level: beginner 156 157 .seealso: MatCreateAIJViennaCL(), MATSEQAIJVIENNACL, MatCreateSeqAIJVIENNACL() 158 M*/ 159 160