Actual source code: mpiaijcusparse.cu

petsc-master 2019-07-20
Report Typos and Errors
  1: #define PETSC_SKIP_SPINLOCK

  3: #include <petscconf.h>
  4:  #include <../src/mat/impls/aij/mpi/mpiaij.h>
  5:  #include <../src/mat/impls/aij/mpi/mpicusparse/mpicusparsematimpl.h>

  7: PetscErrorCode  MatMPIAIJSetPreallocation_MPIAIJCUSPARSE(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:   Mat_MPIAIJCUSPARSE * cusparseStruct = (Mat_MPIAIJCUSPARSE*)b->spptr;
 11:   PetscErrorCode     ierr;
 12:   PetscInt           i;

 15:   PetscLayoutSetUp(B->rmap);
 16:   PetscLayoutSetUp(B->cmap);
 17:   if (d_nnz) {
 18:     for (i=0; i<B->rmap->n; i++) {
 19:       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]);
 20:     }
 21:   }
 22:   if (o_nnz) {
 23:     for (i=0; i<B->rmap->n; i++) {
 24:       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]);
 25:     }
 26:   }
 27:   if (!B->preallocated) {
 28:     /* Explicitly create 2 MATSEQAIJCUSPARSE matrices. */
 29:     MatCreate(PETSC_COMM_SELF,&b->A);
 30:     MatPinToCPU(b->A,B->pinnedtocpu);
 31:     MatSetSizes(b->A,B->rmap->n,B->cmap->n,B->rmap->n,B->cmap->n);
 32:     MatSetType(b->A,MATSEQAIJCUSPARSE);
 33:     PetscLogObjectParent((PetscObject)B,(PetscObject)b->A);
 34:     MatCreate(PETSC_COMM_SELF,&b->B);
 35:     MatPinToCPU(b->B,B->pinnedtocpu);
 36:     MatSetSizes(b->B,B->rmap->n,B->cmap->N,B->rmap->n,B->cmap->N);
 37:     MatSetType(b->B,MATSEQAIJCUSPARSE);
 38:     PetscLogObjectParent((PetscObject)B,(PetscObject)b->B);
 39:   }
 40:   MatSeqAIJSetPreallocation(b->A,d_nz,d_nnz);
 41:   MatSeqAIJSetPreallocation(b->B,o_nz,o_nnz);
 42:   MatCUSPARSESetFormat(b->A,MAT_CUSPARSE_MULT,cusparseStruct->diagGPUMatFormat);
 43:   MatCUSPARSESetFormat(b->B,MAT_CUSPARSE_MULT,cusparseStruct->offdiagGPUMatFormat);
 44:   MatCUSPARSESetHandle(b->A,cusparseStruct->handle);
 45:   MatCUSPARSESetHandle(b->B,cusparseStruct->handle);
 46:   MatCUSPARSESetStream(b->A,cusparseStruct->stream);
 47:   MatCUSPARSESetStream(b->B,cusparseStruct->stream);

 49:   B->preallocated = PETSC_TRUE;
 50:   return(0);
 51: }

 53: PetscErrorCode MatMult_MPIAIJCUSPARSE(Mat A,Vec xx,Vec yy)
 54: {
 55:   /*
 56:      This multiplication sequence is different sequence
 57:      than the CPU version. In particular, the diagonal block
 58:      multiplication kernel is launched in one stream. Then,
 59:      in a separate stream, the data transfers from DeviceToHost
 60:      (with MPI messaging in between), then HostToDevice are
 61:      launched. Once the data transfer stream is synchronized,
 62:      to ensure messaging is complete, the MatMultAdd kernel
 63:      is launched in the original (MatMult) stream to protect
 64:      against race conditions.
 65:   */
 66:   Mat_MPIAIJ     *a = (Mat_MPIAIJ*)A->data;
 68:   PetscInt       nt;

 71:   VecGetLocalSize(xx,&nt);
 72:   if (nt != A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible partition of A (%D) and xx (%D)",A->cmap->n,nt);
 73:   VecScatterInitializeForGPU(a->Mvctx,xx);
 74:   (*a->A->ops->mult)(a->A,xx,yy);
 75:   VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
 76:   VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
 77:   (*a->B->ops->multadd)(a->B,a->lvec,yy,yy);
 78:   VecScatterFinalizeForGPU(a->Mvctx);
 79:   return(0);
 80: }

 82: PetscErrorCode MatMultAdd_MPIAIJCUSPARSE(Mat A,Vec xx,Vec yy,Vec zz)
 83: {
 84:   /*
 85:      This multiplication sequence is different sequence
 86:      than the CPU version. In particular, the diagonal block
 87:      multiplication kernel is launched in one stream. Then,
 88:      in a separate stream, the data transfers from DeviceToHost
 89:      (with MPI messaging in between), then HostToDevice are
 90:      launched. Once the data transfer stream is synchronized,
 91:      to ensure messaging is complete, the MatMultAdd kernel
 92:      is launched in the original (MatMult) stream to protect
 93:      against race conditions.
 94:   */
 95:   Mat_MPIAIJ     *a = (Mat_MPIAIJ*)A->data;
 97:   PetscInt       nt;

100:   VecGetLocalSize(xx,&nt);
101:   if (nt != A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible partition of A (%D) and xx (%D)",A->cmap->n,nt);
102:   VecScatterInitializeForGPU(a->Mvctx,xx);
103:   (*a->A->ops->multadd)(a->A,xx,yy,zz);
104:   VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
105:   VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);
106:   (*a->B->ops->multadd)(a->B,a->lvec,zz,zz);
107:   VecScatterFinalizeForGPU(a->Mvctx);
108:   return(0);
109: }

111: PetscErrorCode MatMultTranspose_MPIAIJCUSPARSE(Mat A,Vec xx,Vec yy)
112: {
113:   /* This multiplication sequence is different sequence
114:      than the CPU version. In particular, the diagonal block
115:      multiplication kernel is launched in one stream. Then,
116:      in a separate stream, the data transfers from DeviceToHost
117:      (with MPI messaging in between), then HostToDevice are
118:      launched. Once the data transfer stream is synchronized,
119:      to ensure messaging is complete, the MatMultAdd kernel
120:      is launched in the original (MatMult) stream to protect
121:      against race conditions.

123:      This sequence should only be called for GPU computation. */
124:   Mat_MPIAIJ     *a = (Mat_MPIAIJ*)A->data;
126:   PetscInt       nt;

129:   VecGetLocalSize(xx,&nt);
130:   if (nt != A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible partition of A (%D) and xx (%D)",A->rmap->n,nt);
131:   VecScatterInitializeForGPU(a->Mvctx,xx);
132:   (*a->B->ops->multtranspose)(a->B,xx,a->lvec);
133:   (*a->A->ops->multtranspose)(a->A,xx,yy);
134:   VecScatterBegin(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE);
135:   VecScatterEnd(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE);
136:   VecScatterFinalizeForGPU(a->Mvctx);
137:   return(0);
138: }

140: PetscErrorCode MatCUSPARSESetFormat_MPIAIJCUSPARSE(Mat A,MatCUSPARSEFormatOperation op,MatCUSPARSEStorageFormat format)
141: {
142:   Mat_MPIAIJ         *a               = (Mat_MPIAIJ*)A->data;
143:   Mat_MPIAIJCUSPARSE * cusparseStruct = (Mat_MPIAIJCUSPARSE*)a->spptr;

146:   switch (op) {
147:   case MAT_CUSPARSE_MULT_DIAG:
148:     cusparseStruct->diagGPUMatFormat = format;
149:     break;
150:   case MAT_CUSPARSE_MULT_OFFDIAG:
151:     cusparseStruct->offdiagGPUMatFormat = format;
152:     break;
153:   case MAT_CUSPARSE_ALL:
154:     cusparseStruct->diagGPUMatFormat    = format;
155:     cusparseStruct->offdiagGPUMatFormat = format;
156:     break;
157:   default:
158:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unsupported operation %d for MatCUSPARSEFormatOperation. Only MAT_CUSPARSE_MULT_DIAG, MAT_CUSPARSE_MULT_DIAG, and MAT_CUSPARSE_MULT_ALL are currently supported.",op);
159:   }
160:   return(0);
161: }

163: PetscErrorCode MatSetFromOptions_MPIAIJCUSPARSE(PetscOptionItems *PetscOptionsObject,Mat A)
164: {
165:   MatCUSPARSEStorageFormat format;
166:   PetscErrorCode           ierr;
167:   PetscBool                flg;
168:   Mat_MPIAIJ               *a = (Mat_MPIAIJ*)A->data;
169:   Mat_MPIAIJCUSPARSE       *cusparseStruct = (Mat_MPIAIJCUSPARSE*)a->spptr;

172:   PetscOptionsHead(PetscOptionsObject,"MPIAIJCUSPARSE options");
173:   if (A->factortype==MAT_FACTOR_NONE) {
174:     PetscOptionsEnum("-mat_cusparse_mult_diag_storage_format","sets storage format of the diagonal blocks of (mpi)aijcusparse gpu matrices for SpMV",
175:                             "MatCUSPARSESetFormat",MatCUSPARSEStorageFormats,(PetscEnum)cusparseStruct->diagGPUMatFormat,(PetscEnum*)&format,&flg);
176:     if (flg) {
177:       MatCUSPARSESetFormat(A,MAT_CUSPARSE_MULT_DIAG,format);
178:     }
179:     PetscOptionsEnum("-mat_cusparse_mult_offdiag_storage_format","sets storage format of the off-diagonal blocks (mpi)aijcusparse gpu matrices for SpMV",
180:                             "MatCUSPARSESetFormat",MatCUSPARSEStorageFormats,(PetscEnum)cusparseStruct->offdiagGPUMatFormat,(PetscEnum*)&format,&flg);
181:     if (flg) {
182:       MatCUSPARSESetFormat(A,MAT_CUSPARSE_MULT_OFFDIAG,format);
183:     }
184:     PetscOptionsEnum("-mat_cusparse_storage_format","sets storage format of the diagonal and off-diagonal blocks (mpi)aijcusparse gpu matrices for SpMV",
185:                             "MatCUSPARSESetFormat",MatCUSPARSEStorageFormats,(PetscEnum)cusparseStruct->diagGPUMatFormat,(PetscEnum*)&format,&flg);
186:     if (flg) {
187:       MatCUSPARSESetFormat(A,MAT_CUSPARSE_ALL,format);
188:     }
189:   }
190:   PetscOptionsTail();
191:   return(0);
192: }

194: PetscErrorCode MatAssemblyEnd_MPIAIJCUSPARSE(Mat A,MatAssemblyType mode)
195: {
197:   Mat_MPIAIJ     *mpiaij;

200:   mpiaij = (Mat_MPIAIJ*)A->data;
201:   MatAssemblyEnd_MPIAIJ(A,mode);
202:   if (!A->was_assembled && mode == MAT_FINAL_ASSEMBLY) {
203:     VecSetType(mpiaij->lvec,VECSEQCUDA);
204:   }
205:   return(0);
206: }

208: PetscErrorCode MatDestroy_MPIAIJCUSPARSE(Mat A)
209: {
210:   PetscErrorCode     ierr;
211:   Mat_MPIAIJ         *a              = (Mat_MPIAIJ*)A->data;
212:   Mat_MPIAIJCUSPARSE *cusparseStruct = (Mat_MPIAIJCUSPARSE*)a->spptr;
213:   cudaError_t        err;
214:   cusparseStatus_t   stat;

217:   try {
218:     MatCUSPARSEClearHandle(a->A);
219:     MatCUSPARSEClearHandle(a->B);
220:     stat = cusparseDestroy(cusparseStruct->handle);CHKERRCUDA(stat);
221:     err = cudaStreamDestroy(cusparseStruct->stream);CHKERRCUDA(err);
222:     delete cusparseStruct;
223:   } catch(char *ex) {
224:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Mat_MPIAIJCUSPARSE error: %s", ex);
225:   }
226:   MatDestroy_MPIAIJ(A);
227:   return(0);
228: }

230: PETSC_EXTERN PetscErrorCode MatCreate_MPIAIJCUSPARSE(Mat A)
231: {
232:   PetscErrorCode     ierr;
233:   Mat_MPIAIJ         *a;
234:   Mat_MPIAIJCUSPARSE * cusparseStruct;
235:   cudaError_t        err;
236:   cusparseStatus_t   stat;

239:   MatCreate_MPIAIJ(A);
240:   PetscObjectComposeFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",MatMPIAIJSetPreallocation_MPIAIJCUSPARSE);
241:   PetscFree(A->defaultvectype);
242:   PetscStrallocpy(VECCUDA,&A->defaultvectype);

244:   a        = (Mat_MPIAIJ*)A->data;
245:   a->spptr = new Mat_MPIAIJCUSPARSE;

247:   cusparseStruct                      = (Mat_MPIAIJCUSPARSE*)a->spptr;
248:   cusparseStruct->diagGPUMatFormat    = MAT_CUSPARSE_CSR;
249:   cusparseStruct->offdiagGPUMatFormat = MAT_CUSPARSE_CSR;
250:   stat = cusparseCreate(&(cusparseStruct->handle));CHKERRCUDA(stat);
251:   err = cudaStreamCreate(&(cusparseStruct->stream));CHKERRCUDA(err);

253:   A->ops->assemblyend    = MatAssemblyEnd_MPIAIJCUSPARSE;
254:   A->ops->mult           = MatMult_MPIAIJCUSPARSE;
255:   A->ops->multadd        = MatMultAdd_MPIAIJCUSPARSE;
256:   A->ops->multtranspose  = MatMultTranspose_MPIAIJCUSPARSE;
257:   A->ops->setfromoptions = MatSetFromOptions_MPIAIJCUSPARSE;
258:   A->ops->destroy        = MatDestroy_MPIAIJCUSPARSE;

260:   PetscObjectChangeTypeName((PetscObject)A,MATMPIAIJCUSPARSE);
261:   PetscObjectComposeFunction((PetscObject)A,"MatCUSPARSESetFormat_C",  MatCUSPARSESetFormat_MPIAIJCUSPARSE);
262:   return(0);
263: }

265: /*@
266:    MatCreateAIJCUSPARSE - Creates a sparse matrix in AIJ (compressed row) format
267:    (the default parallel PETSc format).  This matrix will ultimately pushed down
268:    to NVidia GPUs and use the CUSPARSE library for calculations. For good matrix
269:    assembly performance the user should preallocate the matrix storage by setting
270:    the parameter nz (or the array nnz).  By setting these parameters accurately,
271:    performance during matrix assembly can be increased by more than a factor of 50.

273:    Collective

275:    Input Parameters:
276: +  comm - MPI communicator, set to PETSC_COMM_SELF
277: .  m - number of rows
278: .  n - number of columns
279: .  nz - number of nonzeros per row (same for all rows)
280: -  nnz - array containing the number of nonzeros in the various rows
281:          (possibly different for each row) or NULL

283:    Output Parameter:
284: .  A - the matrix

286:    It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(),
287:    MatXXXXSetPreallocation() paradigm instead of this routine directly.
288:    [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation]

290:    Notes:
291:    If nnz is given then nz is ignored

293:    The AIJ format (also called the Yale sparse matrix format or
294:    compressed row storage), is fully compatible with standard Fortran 77
295:    storage.  That is, the stored row and column indices can begin at
296:    either one (as in Fortran) or zero.  See the users' manual for details.

298:    Specify the preallocated storage with either nz or nnz (not both).
299:    Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
300:    allocation.  For large problems you MUST preallocate memory or you
301:    will get TERRIBLE performance, see the users' manual chapter on matrices.

303:    By default, this format uses inodes (identical nodes) when possible, to
304:    improve numerical efficiency of matrix-vector products and solves. We
305:    search for consecutive rows with the same nonzero structure, thereby
306:    reusing matrix information to achieve increased efficiency.

308:    Level: intermediate

310: .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatCreateAIJ(), MATMPIAIJCUSPARSE, MATAIJCUSPARSE
311: @*/
312: PetscErrorCode  MatCreateAIJCUSPARSE(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)
313: {
315:   PetscMPIInt    size;

318:   MatCreate(comm,A);
319:   MatSetSizes(*A,m,n,M,N);
320:   MPI_Comm_size(comm,&size);
321:   if (size > 1) {
322:     MatSetType(*A,MATMPIAIJCUSPARSE);
323:     MatMPIAIJSetPreallocation(*A,d_nz,d_nnz,o_nz,o_nnz);
324:   } else {
325:     MatSetType(*A,MATSEQAIJCUSPARSE);
326:     MatSeqAIJSetPreallocation(*A,d_nz,d_nnz);
327:   }
328:   return(0);
329: }

331: /*MC
332:    MATAIJCUSPARSE - MATMPIAIJCUSPARSE = "aijcusparse" = "mpiaijcusparse" - A matrix type to be used for sparse matrices.

334:    A matrix type type whose data resides on Nvidia GPUs. These matrices can be in either
335:    CSR, ELL, or Hybrid format. The ELL and HYB formats require CUDA 4.2 or later.
336:    All matrix calculations are performed on Nvidia GPUs using the CUSPARSE library.

338:    This matrix type is identical to MATSEQAIJCUSPARSE when constructed with a single process communicator,
339:    and MATMPIAIJCUSPARSE otherwise.  As a result, for single process communicators,
340:    MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported
341:    for communicators controlling multiple processes.  It is recommended that you call both of
342:    the above preallocation routines for simplicity.

344:    Options Database Keys:
345: +  -mat_type mpiaijcusparse - sets the matrix type to "mpiaijcusparse" during a call to MatSetFromOptions()
346: .  -mat_cusparse_storage_format csr - sets the storage format of diagonal and off-diagonal matrices during a call to MatSetFromOptions(). Other options include ell (ellpack) or hyb (hybrid).
347: .  -mat_cusparse_mult_diag_storage_format csr - sets the storage format of diagonal matrix during a call to MatSetFromOptions(). Other options include ell (ellpack) or hyb (hybrid).
348: -  -mat_cusparse_mult_offdiag_storage_format csr - sets the storage format of off-diagonal matrix during a call to MatSetFromOptions(). Other options include ell (ellpack) or hyb (hybrid).

350:   Level: beginner

352:  .seealso: MatCreateAIJCUSPARSE(), MATSEQAIJCUSPARSE, MatCreateSeqAIJCUSPARSE(), MatCUSPARSESetFormat(), MatCUSPARSEStorageFormat, MatCUSPARSEFormatOperation
353: M
354: M*/