Actual source code: gcreate.c

petsc-master 2020-11-24
Report Typos and Errors

  2: #include <petsc/private/matimpl.h>

  4: PETSC_INTERN PetscErrorCode MatSetBlockSizes_Default(Mat mat,PetscInt rbs, PetscInt cbs)
  5: {
  7:   if (!mat->preallocated) return(0);
  8:   if (mat->rmap->bs > 0 && mat->rmap->bs != rbs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot change row block size %D to %D\n",mat->rmap->bs,rbs);
  9:   if (mat->cmap->bs > 0 && mat->cmap->bs != cbs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot change column block size %D to %D\n",mat->cmap->bs,cbs);
 10:   return(0);
 11: }

 13: PETSC_INTERN PetscErrorCode MatShift_Basic(Mat Y,PetscScalar a)
 14: {
 16:   PetscInt       i,start,end;
 17:   PetscScalar    alpha = a;
 18:   PetscBool      prevoption;

 21:   MatGetOption(Y,MAT_NO_OFF_PROC_ENTRIES,&prevoption);
 22:   MatSetOption(Y,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);
 23:   MatGetOwnershipRange(Y,&start,&end);
 24:   for (i=start; i<end; i++) {
 25:     if (i < Y->cmap->N) {
 26:       MatSetValues(Y,1,&i,1,&i,&alpha,ADD_VALUES);
 27:     }
 28:   }
 29:   MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
 30:   MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
 31:   MatSetOption(Y,MAT_NO_OFF_PROC_ENTRIES,prevoption);
 32:   return(0);
 33: }

 35: /*@
 36:    MatCreate - Creates a matrix where the type is determined
 37:    from either a call to MatSetType() or from the options database
 38:    with a call to MatSetFromOptions(). The default matrix type is
 39:    AIJ, using the routines MatCreateSeqAIJ() or MatCreateAIJ()
 40:    if you do not set a type in the options database. If you never
 41:    call MatSetType() or MatSetFromOptions() it will generate an
 42:    error when you try to use the matrix.

 44:    Collective

 46:    Input Parameter:
 47: .  comm - MPI communicator

 49:    Output Parameter:
 50: .  A - the matrix

 52:    Options Database Keys:
 53: +    -mat_type seqaij   - AIJ type, uses MatCreateSeqAIJ()
 54: .    -mat_type mpiaij   - AIJ type, uses MatCreateAIJ()
 55: .    -mat_type seqdense - dense type, uses MatCreateSeqDense()
 56: .    -mat_type mpidense - dense type, uses MatCreateDense()
 57: .    -mat_type seqbaij  - block AIJ type, uses MatCreateSeqBAIJ()
 58: -    -mat_type mpibaij  - block AIJ type, uses MatCreateBAIJ()

 60:    Even More Options Database Keys:
 61:    See the manpages for particular formats (e.g., MatCreateSeqAIJ())
 62:    for additional format-specific options.

 64:    Level: beginner

 66: .seealso: MatCreateSeqAIJ(), MatCreateAIJ(),
 67:           MatCreateSeqDense(), MatCreateDense(),
 68:           MatCreateSeqBAIJ(), MatCreateBAIJ(),
 69:           MatCreateSeqSBAIJ(), MatCreateSBAIJ(),
 70:           MatConvert()
 71: @*/
 72: PetscErrorCode  MatCreate(MPI_Comm comm,Mat *A)
 73: {
 74:   Mat            B;


 80:   *A = NULL;
 81:   MatInitializePackage();

 83:   PetscHeaderCreate(B,MAT_CLASSID,"Mat","Matrix","Mat",comm,MatDestroy,MatView);
 84:   PetscLayoutCreate(comm,&B->rmap);
 85:   PetscLayoutCreate(comm,&B->cmap);
 86:   PetscStrallocpy(VECSTANDARD,&B->defaultvectype);

 88:   B->congruentlayouts = PETSC_DECIDE;
 89:   B->preallocated     = PETSC_FALSE;
 90:   *A                  = B;
 91:   return(0);
 92: }

 94: /*@
 95:    MatSetErrorIfFailure - Causes Mat to generate an error, for example a zero pivot, is detected.

 97:    Logically Collective on Mat

 99:    Input Parameters:
100: +  mat -  matrix obtained from MatCreate()
101: -  flg - PETSC_TRUE indicates you want the error generated

103:    Level: advanced

105: .seealso: PCSetErrorIfFailure()
106: @*/
107: PetscErrorCode  MatSetErrorIfFailure(Mat mat,PetscBool flg)
108: {
112:   mat->erroriffailure = flg;
113:   return(0);
114: }

116: /*@
117:   MatSetSizes - Sets the local and global sizes, and checks to determine compatibility

119:   Collective on Mat

121:   Input Parameters:
122: +  A - the matrix
123: .  m - number of local rows (or PETSC_DECIDE)
124: .  n - number of local columns (or PETSC_DECIDE)
125: .  M - number of global rows (or PETSC_DETERMINE)
126: -  N - number of global columns (or PETSC_DETERMINE)

128:    Notes:
129:    m (n) and M (N) cannot be both PETSC_DECIDE
130:    If one processor calls this with M (N) of PETSC_DECIDE then all processors must, otherwise the program will hang.

132:    If PETSC_DECIDE is not used for the arguments 'm' and 'n', then the
133:    user must ensure that they are chosen to be compatible with the
134:    vectors. To do this, one first considers the matrix-vector product
135:    'y = A x'. The 'm' that is used in the above routine must match the
136:    local size used in the vector creation routine VecCreateMPI() for 'y'.
137:    Likewise, the 'n' used must match that used as the local size in
138:    VecCreateMPI() for 'x'.

140:    You cannot change the sizes once they have been set.

142:    The sizes must be set before MatSetUp() or MatXXXSetPreallocation() is called.

144:   Level: beginner

146: .seealso: MatGetSize(), PetscSplitOwnership()
147: @*/
148: PetscErrorCode  MatSetSizes(Mat A, PetscInt m, PetscInt n, PetscInt M, PetscInt N)
149: {
154:   if (M > 0 && m > M) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local row size %D cannot be larger than global row size %D",m,M);
155:   if (N > 0 && n > N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local column size %D cannot be larger than global column size %D",n,N);
156:   if ((A->rmap->n >= 0 && A->rmap->N >= 0) && (A->rmap->n != m || (M > 0 && A->rmap->N != M))) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot change/reset row sizes to %D local %D global after previously setting them to %D local %D global",m,M,A->rmap->n,A->rmap->N);
157:   if ((A->cmap->n >= 0 && A->cmap->N >= 0) && (A->cmap->n != n || (N > 0 && A->cmap->N != N))) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot change/reset column sizes to %D local %D global after previously setting them to %D local %D global",n,N,A->cmap->n,A->cmap->N);
158:   A->rmap->n = m;
159:   A->cmap->n = n;
160:   A->rmap->N = M > -1 ? M : A->rmap->N;
161:   A->cmap->N = N > -1 ? N : A->cmap->N;
162:   return(0);
163: }

165: /*@
166:    MatSetFromOptions - Creates a matrix where the type is determined
167:    from the options database. Generates a parallel MPI matrix if the
168:    communicator has more than one processor.  The default matrix type is
169:    AIJ, using the routines MatCreateSeqAIJ() and MatCreateAIJ() if
170:    you do not select a type in the options database.

172:    Collective on Mat

174:    Input Parameter:
175: .  A - the matrix

177:    Options Database Keys:
178: +    -mat_type seqaij   - AIJ type, uses MatCreateSeqAIJ()
179: .    -mat_type mpiaij   - AIJ type, uses MatCreateAIJ()
180: .    -mat_type seqdense - dense type, uses MatCreateSeqDense()
181: .    -mat_type mpidense - dense type, uses MatCreateDense()
182: .    -mat_type seqbaij  - block AIJ type, uses MatCreateSeqBAIJ()
183: -    -mat_type mpibaij  - block AIJ type, uses MatCreateBAIJ()

185:    Even More Options Database Keys:
186:    See the manpages for particular formats (e.g., MatCreateSeqAIJ())
187:    for additional format-specific options.

189:    Level: beginner

191: .seealso: MatCreateSeqAIJ((), MatCreateAIJ(),
192:           MatCreateSeqDense(), MatCreateDense(),
193:           MatCreateSeqBAIJ(), MatCreateBAIJ(),
194:           MatCreateSeqSBAIJ(), MatCreateSBAIJ(),
195:           MatConvert()
196: @*/
197: PetscErrorCode  MatSetFromOptions(Mat B)
198: {
200:   const char     *deft = MATAIJ;
201:   char           type[256];
202:   PetscBool      flg,set;


207:   PetscObjectOptionsBegin((PetscObject)B);

209:   if (B->rmap->bs < 0) {
210:     PetscInt newbs = -1;
211:     PetscOptionsInt("-mat_block_size","Set the blocksize used to store the matrix","MatSetBlockSize",newbs,&newbs,&flg);
212:     if (flg) {
213:       PetscLayoutSetBlockSize(B->rmap,newbs);
214:       PetscLayoutSetBlockSize(B->cmap,newbs);
215:     }
216:   }

218:   PetscOptionsFList("-mat_type","Matrix type","MatSetType",MatList,deft,type,256,&flg);
219:   if (flg) {
220:     MatSetType(B,type);
221:   } else if (!((PetscObject)B)->type_name) {
222:     MatSetType(B,deft);
223:   }

225:   PetscOptionsName("-mat_is_symmetric","Checks if mat is symmetric on MatAssemblyEnd()","MatIsSymmetric",&B->checksymmetryonassembly);
226:   PetscOptionsReal("-mat_is_symmetric","Checks if mat is symmetric on MatAssemblyEnd()","MatIsSymmetric",B->checksymmetrytol,&B->checksymmetrytol,NULL);
227:   PetscOptionsBool("-mat_null_space_test","Checks if provided null space is correct in MatAssemblyEnd()","MatSetNullSpaceTest",B->checknullspaceonassembly,&B->checknullspaceonassembly,NULL);
228:   PetscOptionsBool("-mat_error_if_failure","Generate an error if an error occurs when factoring the matrix","MatSetErrorIfFailure",B->erroriffailure,&B->erroriffailure,NULL);

230:   if (B->ops->setfromoptions) {
231:     (*B->ops->setfromoptions)(PetscOptionsObject,B);
232:   }

234:   flg  = PETSC_FALSE;
235:   PetscOptionsBool("-mat_new_nonzero_location_err","Generate an error if new nonzeros are created in the matrix structure (useful to test preallocation)","MatSetOption",flg,&flg,&set);
236:   if (set) {MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,flg);}
237:   flg  = PETSC_FALSE;
238:   PetscOptionsBool("-mat_new_nonzero_allocation_err","Generate an error if new nonzeros are allocated in the matrix structure (useful to test preallocation)","MatSetOption",flg,&flg,&set);
239:   if (set) {MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,flg);}
240:   flg  = PETSC_FALSE;
241:   PetscOptionsBool("-mat_ignore_zero_entries","For AIJ/IS matrices this will stop zero values from creating a zero location in the matrix","MatSetOption",flg,&flg,&set);
242:   if (set) {MatSetOption(B,MAT_IGNORE_ZERO_ENTRIES,flg);}

244:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
245:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)B);
246:   PetscOptionsEnd();
247:   return(0);
248: }

250: /*@C
251:    MatXAIJSetPreallocation - set preallocation for serial and parallel AIJ, BAIJ, and SBAIJ matrices and their unassembled versions.

253:    Collective on Mat

255:    Input Arguments:
256: +  A - matrix being preallocated
257: .  bs - block size
258: .  dnnz - number of nonzero column blocks per block row of diagonal part of parallel matrix
259: .  onnz - number of nonzero column blocks per block row of off-diagonal part of parallel matrix
260: .  dnnzu - number of nonzero column blocks per block row of upper-triangular part of diagonal part of parallel matrix
261: -  onnzu - number of nonzero column blocks per block row of upper-triangular part of off-diagonal part of parallel matrix

263:    Level: beginner

265: .seealso: MatSeqAIJSetPreallocation(), MatMPIAIJSetPreallocation(), MatSeqBAIJSetPreallocation(), MatMPIBAIJSetPreallocation(), MatSeqSBAIJSetPreallocation(), MatMPISBAIJSetPreallocation(),
266:           PetscSplitOwnership()
267: @*/
268: PetscErrorCode MatXAIJSetPreallocation(Mat A,PetscInt bs,const PetscInt dnnz[],const PetscInt onnz[],const PetscInt dnnzu[],const PetscInt onnzu[])
269: {
271:   PetscInt       cbs;
272:   void           (*aij)(void);
273:   void           (*is)(void);
274:   void           (*hyp)(void) = NULL;

277:   if (bs != PETSC_DECIDE) { /* don't mess with an already set block size */
278:     MatSetBlockSize(A,bs);
279:   }
280:   PetscLayoutSetUp(A->rmap);
281:   PetscLayoutSetUp(A->cmap);
282:   MatGetBlockSizes(A,&bs,&cbs);
283:   /* these routines assumes bs == cbs, this should be checked somehow */
284:   MatSeqBAIJSetPreallocation(A,bs,0,dnnz);
285:   MatMPIBAIJSetPreallocation(A,bs,0,dnnz,0,onnz);
286:   MatSeqSBAIJSetPreallocation(A,bs,0,dnnzu);
287:   MatMPISBAIJSetPreallocation(A,bs,0,dnnzu,0,onnzu);
288:   /*
289:     In general, we have to do extra work to preallocate for scalar (AIJ) or unassembled (IS) matrices so we check whether it will do any
290:     good before going on with it.
291:   */
292:   PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);
293:   PetscObjectQueryFunction((PetscObject)A,"MatISSetPreallocation_C",&is);
294: #if defined(PETSC_HAVE_HYPRE)
295:   PetscObjectQueryFunction((PetscObject)A,"MatHYPRESetPreallocation_C",&hyp);
296: #endif
297:   if (!aij && !is && !hyp) {
298:     PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);
299:   }
300:   if (aij || is || hyp) {
301:     if (bs == cbs && bs == 1) {
302:       MatSeqAIJSetPreallocation(A,0,dnnz);
303:       MatMPIAIJSetPreallocation(A,0,dnnz,0,onnz);
304:       MatISSetPreallocation(A,0,dnnz,0,onnz);
305: #if defined(PETSC_HAVE_HYPRE)
306:       MatHYPRESetPreallocation(A,0,dnnz,0,onnz);
307: #endif
308:     } else { /* Convert block-row precallocation to scalar-row */
309:       PetscInt i,m,*sdnnz,*sonnz;
310:       MatGetLocalSize(A,&m,NULL);
311:       PetscMalloc2((!!dnnz)*m,&sdnnz,(!!onnz)*m,&sonnz);
312:       for (i=0; i<m; i++) {
313:         if (dnnz) sdnnz[i] = dnnz[i/bs] * cbs;
314:         if (onnz) sonnz[i] = onnz[i/bs] * cbs;
315:       }
316:       MatSeqAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL);
317:       MatMPIAIJSetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);
318:       MatISSetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);
319: #if defined(PETSC_HAVE_HYPRE)
320:       MatHYPRESetPreallocation(A,0,dnnz ? sdnnz : NULL,0,onnz ? sonnz : NULL);
321: #endif
322:       PetscFree2(sdnnz,sonnz);
323:     }
324:   }
325:   return(0);
326: }

328: /*
329:         Merges some information from Cs header to A; the C object is then destroyed

331:         This is somewhat different from MatHeaderReplace() it would be nice to merge the code
332: */
333: PetscErrorCode MatHeaderMerge(Mat A,Mat *C)
334: {
336:   PetscInt       refct;
337:   PetscOps       Abops;
338:   struct _MatOps Aops;
339:   char           *mtype,*mname,*mprefix;
340:   Mat_Product    *product;

345:   if (A == *C) return(0);
347:   /* save the parts of A we need */
348:   Abops = ((PetscObject)A)->bops[0];
349:   Aops  = A->ops[0];
350:   refct = ((PetscObject)A)->refct;
351:   mtype = ((PetscObject)A)->type_name;
352:   mname = ((PetscObject)A)->name;
353:   mprefix = ((PetscObject)A)->prefix;
354:   product = A->product;

356:   /* zero these so the destroy below does not free them */
357:   ((PetscObject)A)->type_name = NULL;
358:   ((PetscObject)A)->name      = NULL;

360:   /* free all the interior data structures from mat */
361:   (*A->ops->destroy)(A);

363:   PetscFree(A->defaultvectype);
364:   PetscLayoutDestroy(&A->rmap);
365:   PetscLayoutDestroy(&A->cmap);
366:   PetscFunctionListDestroy(&((PetscObject)A)->qlist);
367:   PetscObjectListDestroy(&((PetscObject)A)->olist);

369:   /* copy C over to A */
370:   PetscMemcpy(A,*C,sizeof(struct _p_Mat));

372:   /* return the parts of A we saved */
373:   ((PetscObject)A)->bops[0]   = Abops;
374:   A->ops[0]                   = Aops;
375:   ((PetscObject)A)->refct     = refct;
376:   ((PetscObject)A)->type_name = mtype;
377:   ((PetscObject)A)->name      = mname;
378:   ((PetscObject)A)->prefix    = mprefix;
379:   A->product                  = product;

381:   /* since these two are copied into A we do not want them destroyed in C */
382:   ((PetscObject)*C)->qlist = NULL;
383:   ((PetscObject)*C)->olist = NULL;

385:   PetscHeaderDestroy(C);
386:   return(0);
387: }
388: /*
389:         Replace A's header with that of C; the C object is then destroyed

391:         This is essentially code moved from MatDestroy()

393:         This is somewhat different from MatHeaderMerge() it would be nice to merge the code

395:         Used in DM hence is declared PETSC_EXTERN
396: */
397: PETSC_EXTERN PetscErrorCode MatHeaderReplace(Mat A,Mat *C)
398: {
399:   PetscErrorCode   ierr;
400:   PetscInt         refct;
401:   PetscObjectState state;
402:   struct _p_Mat    buffer;
403:   MatStencilInfo   stencil;

408:   if (A == *C) return(0);
410:   if (((PetscObject)*C)->refct != 1) SETERRQ1(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Object C has refct %D > 1, would leave hanging reference",((PetscObject)*C)->refct);

412:   /* swap C and A */
413:   refct   = ((PetscObject)A)->refct;
414:   state   = ((PetscObject)A)->state;
415:   stencil = A->stencil;
416:   PetscMemcpy(&buffer,A,sizeof(struct _p_Mat));
417:   PetscMemcpy(A,*C,sizeof(struct _p_Mat));
418:   PetscMemcpy(*C,&buffer,sizeof(struct _p_Mat));
419:   ((PetscObject)A)->refct   = refct;
420:   ((PetscObject)A)->state   = state + 1;
421:   A->stencil                = stencil;

423:   ((PetscObject)*C)->refct = 1;
424:   MatShellSetOperation(*C,MATOP_DESTROY,(void(*)(void))NULL);
425:   MatDestroy(C);
426:   return(0);
427: }

429: /*@
430:      MatBindToCPU - marks a matrix to temporarily stay on the CPU and perform computations on the CPU

432:    Input Parameters:
433: +   A - the matrix
434: -   flg - bind to the CPU if value of PETSC_TRUE

436:    Level: intermediate
437: @*/
438: PetscErrorCode MatBindToCPU(Mat A,PetscBool flg)
439: {
440: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)

446:   if (A->boundtocpu == flg) return(0);
447:   A->boundtocpu = flg;
448:   if (A->ops->bindtocpu) {
449:     (*A->ops->bindtocpu)(A,flg);
450:   }
451:   return(0);
452: #else
456:   return(0);
457: #endif
458: }

460: PetscErrorCode MatSetValuesCOO_Basic(Mat A,const PetscScalar coo_v[],InsertMode imode)
461: {
462:   IS             is_coo_i,is_coo_j;
463:   const PetscInt *coo_i,*coo_j;
464:   PetscInt       n,n_i,n_j;
465:   PetscScalar    zero = 0.;

469:   PetscObjectQuery((PetscObject)A,"__PETSc_coo_i",(PetscObject*)&is_coo_i);
470:   PetscObjectQuery((PetscObject)A,"__PETSc_coo_j",(PetscObject*)&is_coo_j);
471:   if (!is_coo_i) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_COR,"Missing coo_i IS");
472:   if (!is_coo_j) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_COR,"Missing coo_j IS");
473:   ISGetLocalSize(is_coo_i,&n_i);
474:   ISGetLocalSize(is_coo_j,&n_j);
475:   if (n_i != n_j)  SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_COR,"Wrong local size %D != %D",n_i,n_j);
476:   ISGetIndices(is_coo_i,&coo_i);
477:   ISGetIndices(is_coo_j,&coo_j);
478:   for (n = 0; n < n_i; n++) {
479:     MatSetValue(A,coo_i[n],coo_j[n],coo_v ? coo_v[n] : zero,imode);
480:   }
481:   ISRestoreIndices(is_coo_i,&coo_i);
482:   ISRestoreIndices(is_coo_j,&coo_j);
483:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
484:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
485:   return(0);
486: }

488: PetscErrorCode MatSetPreallocationCOO_Basic(Mat A,PetscInt ncoo,const PetscInt coo_i[],const PetscInt coo_j[])
489: {
490:   Mat            preallocator;
491:   IS             is_coo_i,is_coo_j;
492:   PetscScalar    zero = 0.0;
493:   PetscInt       n;

497:   PetscLayoutSetUp(A->rmap);
498:   PetscLayoutSetUp(A->cmap);
499:   MatCreate(PetscObjectComm((PetscObject)A),&preallocator);
500:   MatSetType(preallocator,MATPREALLOCATOR);
501:   MatSetSizes(preallocator,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
502:   MatSetLayouts(preallocator,A->rmap,A->cmap);
503:   MatSetUp(preallocator);
504:   for (n = 0; n < ncoo; n++) {
505:     MatSetValue(preallocator,coo_i[n],coo_j[n],zero,INSERT_VALUES);
506:   }
507:   MatAssemblyBegin(preallocator,MAT_FINAL_ASSEMBLY);
508:   MatAssemblyEnd(preallocator,MAT_FINAL_ASSEMBLY);
509:   MatPreallocatorPreallocate(preallocator,PETSC_TRUE,A);
510:   MatDestroy(&preallocator);
511:   ISCreateGeneral(PETSC_COMM_SELF,ncoo,coo_i,PETSC_COPY_VALUES,&is_coo_i);
512:   ISCreateGeneral(PETSC_COMM_SELF,ncoo,coo_j,PETSC_COPY_VALUES,&is_coo_j);
513:   PetscObjectCompose((PetscObject)A,"__PETSc_coo_i",(PetscObject)is_coo_i);
514:   PetscObjectCompose((PetscObject)A,"__PETSc_coo_j",(PetscObject)is_coo_j);
515:   ISDestroy(&is_coo_i);
516:   ISDestroy(&is_coo_j);
517:   return(0);
518: }

520: /*@C
521:    MatSetPreallocationCOO - set preallocation for matrices using a coordinate format of the entries

523:    Collective on Mat

525:    Input Arguments:
526: +  A - matrix being preallocated
527: .  ncoo - number of entries in the locally owned part of the parallel matrix
528: .  coo_i - row indices
529: -  coo_j - column indices

531:    Level: beginner

533:    Notes: Entries can be repeated. Currently optimized for cuSPARSE matrices only.

535: .seealso: MatSetValuesCOO(), MatSeqAIJSetPreallocation(), MatMPIAIJSetPreallocation(), MatSeqBAIJSetPreallocation(), MatMPIBAIJSetPreallocation(), MatSeqSBAIJSetPreallocation(), MatMPISBAIJSetPreallocation()
536: @*/
537: PetscErrorCode MatSetPreallocationCOO(Mat A,PetscInt ncoo,const PetscInt coo_i[],const PetscInt coo_j[])
538: {
539:   PetscErrorCode (*f)(Mat,PetscInt,const PetscInt[],const PetscInt[]) = NULL;

547:   PetscLayoutSetUp(A->rmap);
548:   PetscLayoutSetUp(A->cmap);
549:   if (PetscDefined(USE_DEBUG)) {
550:     PetscInt i;
551:     for (i = 0; i < ncoo; i++) {
552:       if (coo_i[i] < A->rmap->rstart || coo_i[i] >= A->rmap->rend) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_USER,"Invalid row index %D! Must be in [%D,%D)",coo_i[i],A->rmap->rstart,A->rmap->rend);
553:       if (coo_j[i] < 0 || coo_j[i] >= A->cmap->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Invalid col index %D! Must be in [0,%D)",coo_j[i],A->cmap->N);
554:     }
555:   }
556:   PetscObjectQueryFunction((PetscObject)A,"MatSetPreallocationCOO_C",&f);
557:   if (f) {
558:     (*f)(A,ncoo,coo_i,coo_j);
559:   } else { /* allow fallback, very slow */
560:     MatSetPreallocationCOO_Basic(A,ncoo,coo_i,coo_j);
561:   }
562:   return(0);
563: }

565: /*@C
566:    MatSetValuesCOO - set values at once in a matrix preallocated using MatSetPreallocationCOO

568:    Collective on Mat

570:    Input Arguments:
571: +  A - matrix being preallocated
572: .  coo_v - the matrix values
573: -  imode - the insert mode

575:    Level: beginner

577:    Notes: The values must follow the order of the indices prescribed with MatSetPreallocationCOO().
578:           Currently optimized for cuSPARSE matrices only.

580: .seealso: MatSetPreallocationCOO(), InsertMode, INSERT_VALUES, ADD_VALUES
581: @*/
582: PetscErrorCode MatSetValuesCOO(Mat A, const PetscScalar coo_v[], InsertMode imode)
583: {
584:   PetscErrorCode (*f)(Mat,const PetscScalar[],InsertMode) = NULL;

590:   MatCheckPreallocated(A,1);
592:   PetscObjectQueryFunction((PetscObject)A,"MatSetValuesCOO_C",&f);
593:   if (f) {
594:     (*f)(A,coo_v,imode);
595:   } else { /* allow fallback */
596:     MatSetValuesCOO_Basic(A,coo_v,imode);
597:   }
598:   return(0);
599: }