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