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