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