Actual source code: gcreate.c

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

  3: PETSC_INTERN PetscErrorCode MatSetBlockSizes_Default(Mat mat, PetscInt rbs, PetscInt cbs)
  4: {
  5:   if (!mat->preallocated) return 0;
  8:   return 0;
  9: }

 11: PETSC_INTERN PetscErrorCode MatShift_Basic(Mat Y, PetscScalar a)
 12: {
 13:   PetscInt    i, start, end;
 14:   PetscScalar alpha = a;
 15:   PetscBool   prevoption;

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

 29: /*@
 30:    MatCreate - Creates a matrix where the type is determined
 31:    from either a call to `MatSetType()` or from the options database
 32:    with a call to `MatSetFromOptions()`. The default matrix type is
 33:    `MATAIJ`, using the routines `MatCreateSeqAIJ()` or `MatCreateAIJ()`
 34:    if you do not set a type in the options database. If you never
 35:    call `MatSetType()` or `MatSetFromOptions()` it will generate an
 36:    error when you try to use the matrix.

 38:    Collective

 40:    Input Parameter:
 41: .  comm - MPI communicator

 43:    Output Parameter:
 44: .  A - the matrix

 46:    Options Database Keys:
 47: +    -mat_type seqaij   - `MATSEQAIJ` type, uses `MatCreateSeqAIJ()`
 48: .    -mat_type mpiaij   - `MATMPIAIJ` type, uses `MatCreateAIJ()`
 49: .    -mat_type seqdense - `MATSEQDENSE`, uses `MatCreateSeqDense()`
 50: .    -mat_type mpidense - `MATMPIDENSE` type, uses `MatCreateDense()`
 51: .    -mat_type seqbaij  - `MATSEQBAIJ` type, uses `MatCreateSeqBAIJ()`
 52: -    -mat_type mpibaij  - `MATMPIBAIJ` type, uses `MatCreateBAIJ()`

 54:    Even More Options Database Keys:
 55:    See the manpages for particular formats (e.g., `MatCreateSeqAIJ()`)
 56:    for additional format-specific options.

 58:    Level: beginner

 60: .seealso: `MatCreateSeqAIJ()`, `MatCreateAIJ()`,
 61:           `MatCreateSeqDense()`, `MatCreateDense()`,
 62:           `MatCreateSeqBAIJ()`, `MatCreateBAIJ()`,
 63:           `MatCreateSeqSBAIJ()`, `MatCreateSBAIJ()`,
 64:           `MatConvert()`
 65: @*/
 66: PetscErrorCode MatCreate(MPI_Comm comm, Mat *A)
 67: {
 68:   Mat B;


 72:   *A = NULL;
 73:   MatInitializePackage();

 75:   PetscHeaderCreate(B, MAT_CLASSID, "Mat", "Matrix", "Mat", comm, MatDestroy, MatView);
 76:   PetscLayoutCreate(comm, &B->rmap);
 77:   PetscLayoutCreate(comm, &B->cmap);
 78:   PetscStrallocpy(VECSTANDARD, &B->defaultvectype);
 79:   PetscStrallocpy(PETSCRANDER48, &B->defaultrandtype);

 81:   B->symmetric                   = PETSC_BOOL3_UNKNOWN;
 82:   B->hermitian                   = PETSC_BOOL3_UNKNOWN;
 83:   B->structurally_symmetric      = PETSC_BOOL3_UNKNOWN;
 84:   B->spd                         = PETSC_BOOL3_UNKNOWN;
 85:   B->symmetry_eternal            = PETSC_FALSE;
 86:   B->structural_symmetry_eternal = PETSC_FALSE;

 88:   B->congruentlayouts = PETSC_DECIDE;
 89:   B->preallocated     = PETSC_FALSE;
 90: #if defined(PETSC_HAVE_DEVICE)
 91:   B->boundtocpu = PETSC_TRUE;
 92: #endif
 93:   *A = B;
 94:   return 0;
 95: }

 97: /*@
 98:    MatSetErrorIfFailure - Causes `Mat` to generate an immediate error, for example a zero pivot, is detected.

100:    Logically Collective

102:    Input Parameters:
103: +  mat -  matrix obtained from `MatCreate()`
104: -  flg - `PETSC_TRUE` indicates you want the error generated

106:    Level: advanced

108:    Note:
109:    If this flag is not set then the matrix operation will note the error and continue. The error may cause a later `PC` or `KSP` error
110:    or result in a `KSPConvergedReason` indicating the method did not converge.

112: .seealso: `PCSetErrorIfFailure()`, `KSPConvergedReason`, `SNESConvergedReason`
113: @*/
114: PetscErrorCode MatSetErrorIfFailure(Mat mat, PetscBool flg)
115: {
118:   mat->erroriffailure = flg;
119:   return 0;
120: }

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

125:   Collective

127:   Input Parameters:
128: +  A - the matrix
129: .  m - number of local rows (or `PETSC_DECIDE`)
130: .  n - number of local columns (or `PETSC_DECIDE`)
131: .  M - number of global rows (or `PETSC_DETERMINE`)
132: -  N - number of global columns (or `PETSC_DETERMINE`)

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

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

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

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

150:   Level: beginner

152: .seealso: `MatGetSize()`, `PetscSplitOwnership()`
153: @*/
154: PetscErrorCode MatSetSizes(Mat A, PetscInt m, PetscInt n, PetscInt M, PetscInt N)
155: {
162:              A->rmap->n, A->rmap->N);
164:              A->cmap->n, A->cmap->N);
165:   A->rmap->n = m;
166:   A->cmap->n = n;
167:   A->rmap->N = M > -1 ? M : A->rmap->N;
168:   A->cmap->N = N > -1 ? N : A->cmap->N;
169:   return 0;
170: }

172: /*@
173:    MatSetFromOptions - Creates a matrix where the type is determined
174:    from the options database. Generates a parallel MPI matrix if the
175:    communicator has more than one processor.  The default matrix type is
176:    `MATAIJ`, using the routines `MatCreateSeqAIJ()` and `MatCreateAIJ()` if
177:    you do not select a type in the options database.

179:    Collective

181:    Input Parameter:
182: .  A - the matrix

184:    Options Database Keys:
185: +    -mat_type seqaij   - `MATSEQAIJ` type, uses `MatCreateSeqAIJ()`
186: .    -mat_type mpiaij   - `MATMPIAIJ` type, uses `MatCreateAIJ()`
187: .    -mat_type seqdense - `MATSEQDENSE` type, uses `MatCreateSeqDense()`
188: .    -mat_type mpidense - `MATMPIDENSE`, uses `MatCreateDense()`
189: .    -mat_type seqbaij  - `MATSEQBAIJ`, uses `MatCreateSeqBAIJ()`
190: -    -mat_type mpibaij  - `MATMPIBAIJ`, uses `MatCreateBAIJ()`

192:    Even More Options Database Keys:
193:    See the manpages for particular formats (e.g., `MatCreateSeqAIJ()`)
194:    for additional format-specific options.

196:    Level: beginner

198: .seealso: `MatCreateSeqAIJ(()`, `MatCreateAIJ()`,
199:           `MatCreateSeqDense()`, `MatCreateDense()`,
200:           `MatCreateSeqBAIJ()`, `MatCreateBAIJ()`,
201:           `MatCreateSeqSBAIJ()`, `MatCreateSBAIJ()`,
202:           `MatConvert()`
203: @*/
204: PetscErrorCode MatSetFromOptions(Mat B)
205: {
206:   const char *deft = MATAIJ;
207:   char        type[256];
208:   PetscBool   flg, set;
209:   PetscInt    bind_below = 0;


213:   PetscObjectOptionsBegin((PetscObject)B);

215:   if (B->rmap->bs < 0) {
216:     PetscInt newbs = -1;
217:     PetscOptionsInt("-mat_block_size", "Set the blocksize used to store the matrix", "MatSetBlockSize", newbs, &newbs, &flg);
218:     if (flg) {
219:       PetscLayoutSetBlockSize(B->rmap, newbs);
220:       PetscLayoutSetBlockSize(B->cmap, newbs);
221:     }
222:   }

224:   PetscOptionsFList("-mat_type", "Matrix type", "MatSetType", MatList, deft, type, 256, &flg);
225:   if (flg) {
226:     MatSetType(B, type);
227:   } else if (!((PetscObject)B)->type_name) {
228:     MatSetType(B, deft);
229:   }

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

236:   PetscTryTypeMethod(B, setfromoptions, PetscOptionsObject);

238:   flg = PETSC_FALSE;
239:   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);
240:   if (set) MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, flg);
241:   flg = PETSC_FALSE;
242:   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);
243:   if (set) MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, flg);
244:   flg = PETSC_FALSE;
245:   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);
246:   if (set) MatSetOption(B, MAT_IGNORE_ZERO_ENTRIES, flg);

248:   flg = PETSC_FALSE;
249:   PetscOptionsBool("-mat_form_explicit_transpose", "Hint to form an explicit transpose for operations like MatMultTranspose", "MatSetOption", flg, &flg, &set);
250:   if (set) MatSetOption(B, MAT_FORM_EXPLICIT_TRANSPOSE, flg);

252:   /* Bind to CPU if below a user-specified size threshold.
253:    * This perhaps belongs in the options for the GPU Mat types, but MatBindToCPU() does nothing when called on non-GPU types,
254:    * and putting it here makes is more maintainable than duplicating this for all. */
255:   PetscOptionsInt("-mat_bind_below", "Set the size threshold (in local rows) below which the Mat is bound to the CPU", "MatBindToCPU", bind_below, &bind_below, &flg);
256:   if (flg && B->rmap->n < bind_below) MatBindToCPU(B, PETSC_TRUE);

258:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
259:   PetscObjectProcessOptionsHandlers((PetscObject)B, PetscOptionsObject);
260:   PetscOptionsEnd();
261:   return 0;
262: }

264: /*@C
265:    MatXAIJSetPreallocation - set preallocation for serial and parallel `MATAIJ`, `MATBAIJ`, and `MATSBAIJ` matrices and their unassembled versions.

267:    Collective

269:    Input Parameters:
270: +  A - matrix being preallocated
271: .  bs - block size
272: .  dnnz - number of nonzero column blocks per block row of diagonal part of parallel matrix
273: .  onnz - number of nonzero column blocks per block row of off-diagonal part of parallel matrix
274: .  dnnzu - number of nonzero column blocks per block row of upper-triangular part of diagonal part of parallel matrix
275: -  onnzu - number of nonzero column blocks per block row of upper-triangular part of off-diagonal part of parallel matrix

277:    Level: beginner

279: .seealso: `MatSeqAIJSetPreallocation()`, `MatMPIAIJSetPreallocation()`, `MatSeqBAIJSetPreallocation()`, `MatMPIBAIJSetPreallocation()`, `MatSeqSBAIJSetPreallocation()`, `MatMPISBAIJSetPreallocation()`,
280:           `PetscSplitOwnership()`
281: @*/
282: PetscErrorCode MatXAIJSetPreallocation(Mat A, PetscInt bs, const PetscInt dnnz[], const PetscInt onnz[], const PetscInt dnnzu[], const PetscInt onnzu[])
283: {
284:   PetscInt cbs;
285:   void (*aij)(void);
286:   void (*is)(void);
287:   void (*hyp)(void) = NULL;

289:   if (bs != PETSC_DECIDE) { /* don't mess with an already set block size */
290:     MatSetBlockSize(A, bs);
291:   }
292:   PetscLayoutSetUp(A->rmap);
293:   PetscLayoutSetUp(A->cmap);
294:   MatGetBlockSizes(A, &bs, &cbs);
295:   /* these routines assumes bs == cbs, this should be checked somehow */
296:   MatSeqBAIJSetPreallocation(A, bs, 0, dnnz);
297:   MatMPIBAIJSetPreallocation(A, bs, 0, dnnz, 0, onnz);
298:   MatSeqSBAIJSetPreallocation(A, bs, 0, dnnzu);
299:   MatMPISBAIJSetPreallocation(A, bs, 0, dnnzu, 0, onnzu);
300:   /*
301:     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
302:     good before going on with it.
303:   */
304:   PetscObjectQueryFunction((PetscObject)A, "MatMPIAIJSetPreallocation_C", &aij);
305:   PetscObjectQueryFunction((PetscObject)A, "MatISSetPreallocation_C", &is);
306: #if defined(PETSC_HAVE_HYPRE)
307:   PetscObjectQueryFunction((PetscObject)A, "MatHYPRESetPreallocation_C", &hyp);
308: #endif
309:   if (!aij && !is && !hyp) PetscObjectQueryFunction((PetscObject)A, "MatSeqAIJSetPreallocation_C", &aij);
310:   if (aij || is || hyp) {
311:     if (bs == cbs && bs == 1) {
312:       MatSeqAIJSetPreallocation(A, 0, dnnz);
313:       MatMPIAIJSetPreallocation(A, 0, dnnz, 0, onnz);
314:       MatISSetPreallocation(A, 0, dnnz, 0, onnz);
315: #if defined(PETSC_HAVE_HYPRE)
316:       MatHYPRESetPreallocation(A, 0, dnnz, 0, onnz);
317: #endif
318:     } else { /* Convert block-row precallocation to scalar-row */
319:       PetscInt i, m, *sdnnz, *sonnz;
320:       MatGetLocalSize(A, &m, NULL);
321:       PetscMalloc2((!!dnnz) * m, &sdnnz, (!!onnz) * m, &sonnz);
322:       for (i = 0; i < m; i++) {
323:         if (dnnz) sdnnz[i] = dnnz[i / bs] * cbs;
324:         if (onnz) sonnz[i] = onnz[i / bs] * cbs;
325:       }
326:       MatSeqAIJSetPreallocation(A, 0, dnnz ? sdnnz : NULL);
327:       MatMPIAIJSetPreallocation(A, 0, dnnz ? sdnnz : NULL, 0, onnz ? sonnz : NULL);
328:       MatISSetPreallocation(A, 0, dnnz ? sdnnz : NULL, 0, onnz ? sonnz : NULL);
329: #if defined(PETSC_HAVE_HYPRE)
330:       MatHYPRESetPreallocation(A, 0, dnnz ? sdnnz : NULL, 0, onnz ? sonnz : NULL);
331: #endif
332:       PetscFree2(sdnnz, sonnz);
333:     }
334:   }
335:   return 0;
336: }

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

341:         This is somewhat different from MatHeaderReplace() it would be nice to merge the code
342: */
343: PetscErrorCode MatHeaderMerge(Mat A, Mat *C)
344: {
345:   PetscInt         refct;
346:   PetscOps         Abops;
347:   struct _MatOps   Aops;
348:   char            *mtype, *mname, *mprefix;
349:   Mat_Product     *product;
350:   Mat_Redundant   *redundant;
351:   PetscObjectState state;

355:   if (A == *C) return 0;
357:   /* save the parts of A we need */
358:   Abops     = ((PetscObject)A)->bops[0];
359:   Aops      = A->ops[0];
360:   refct     = ((PetscObject)A)->refct;
361:   mtype     = ((PetscObject)A)->type_name;
362:   mname     = ((PetscObject)A)->name;
363:   state     = ((PetscObject)A)->state;
364:   mprefix   = ((PetscObject)A)->prefix;
365:   product   = A->product;
366:   redundant = A->redundant;

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

372:   /*
373:      free all the interior data structures from mat
374:      cannot use PetscUseTypeMethod(A,destroy); because compiler
375:      thinks it may print NULL type_name and name
376:   */
377:   PetscTryTypeMethod(A, destroy);

379:   PetscFree(A->defaultvectype);
380:   PetscFree(A->defaultrandtype);
381:   PetscLayoutDestroy(&A->rmap);
382:   PetscLayoutDestroy(&A->cmap);
383:   PetscFunctionListDestroy(&((PetscObject)A)->qlist);
384:   PetscObjectListDestroy(&((PetscObject)A)->olist);
385:   PetscComposedQuantitiesDestroy((PetscObject)A);

387:   /* copy C over to A */
388:   PetscFree(A->factorprefix);
389:   PetscMemcpy(A, *C, sizeof(struct _p_Mat));

391:   /* return the parts of A we saved */
392:   ((PetscObject)A)->bops[0]   = Abops;
393:   A->ops[0]                   = Aops;
394:   ((PetscObject)A)->refct     = refct;
395:   ((PetscObject)A)->type_name = mtype;
396:   ((PetscObject)A)->name      = mname;
397:   ((PetscObject)A)->prefix    = mprefix;
398:   ((PetscObject)A)->state     = state + 1;
399:   A->product                  = product;
400:   A->redundant                = redundant;

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

406:   PetscHeaderDestroy(C);
407:   return 0;
408: }
409: /*
410:         Replace A's header with that of C; the C object is then destroyed

412:         This is essentially code moved from MatDestroy()

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

416:         Used in DM hence is declared PETSC_EXTERN
417: */
418: PETSC_EXTERN PetscErrorCode MatHeaderReplace(Mat A, Mat *C)
419: {
420:   PetscInt         refct;
421:   PetscObjectState state;
422:   struct _p_Mat    buffer;
423:   MatStencilInfo   stencil;

427:   if (A == *C) return 0;

431:   /* swap C and A */
432:   refct   = ((PetscObject)A)->refct;
433:   state   = ((PetscObject)A)->state;
434:   stencil = A->stencil;
435:   PetscMemcpy(&buffer, A, sizeof(struct _p_Mat));
436:   PetscMemcpy(A, *C, sizeof(struct _p_Mat));
437:   PetscMemcpy(*C, &buffer, sizeof(struct _p_Mat));
438:   ((PetscObject)A)->refct = refct;
439:   ((PetscObject)A)->state = state + 1;
440:   A->stencil              = stencil;

442:   ((PetscObject)*C)->refct = 1;
443:   MatShellSetOperation(*C, MATOP_DESTROY, (void (*)(void))NULL);
444:   MatDestroy(C);
445:   return 0;
446: }

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

451:    Logically collective

453:    Input Parameters:
454: +   A - the matrix
455: -   flg - bind to the CPU if value of `PETSC_TRUE`

457:    Level: intermediate

459: .seealso: `MatBoundToCPU()`
460: @*/
461: PetscErrorCode MatBindToCPU(Mat A, PetscBool flg)
462: {
465: #if defined(PETSC_HAVE_DEVICE)
466:   if (A->boundtocpu == flg) return 0;
467:   A->boundtocpu = flg;
468:   PetscTryTypeMethod(A, bindtocpu, flg);
469: #endif
470:   return 0;
471: }

473: /*@
474:      MatBoundToCPU - query if a matrix is bound to the CPU

476:    Input Parameter:
477: .   A - the matrix

479:    Output Parameter:
480: .   flg - the logical flag

482:    Level: intermediate

484: .seealso: `MatBindToCPU()`
485: @*/
486: PetscErrorCode MatBoundToCPU(Mat A, PetscBool *flg)
487: {
490: #if defined(PETSC_HAVE_DEVICE)
491:   *flg = A->boundtocpu;
492: #else
493:   *flg = PETSC_TRUE;
494: #endif
495:   return 0;
496: }

498: PetscErrorCode MatSetValuesCOO_Basic(Mat A, const PetscScalar coo_v[], InsertMode imode)
499: {
500:   IS              is_coo_i, is_coo_j;
501:   const PetscInt *coo_i, *coo_j;
502:   PetscInt        n, n_i, n_j;
503:   PetscScalar     zero = 0.;

505:   PetscObjectQuery((PetscObject)A, "__PETSc_coo_i", (PetscObject *)&is_coo_i);
506:   PetscObjectQuery((PetscObject)A, "__PETSc_coo_j", (PetscObject *)&is_coo_j);
509:   ISGetLocalSize(is_coo_i, &n_i);
510:   ISGetLocalSize(is_coo_j, &n_j);
512:   ISGetIndices(is_coo_i, &coo_i);
513:   ISGetIndices(is_coo_j, &coo_j);
514:   if (imode != ADD_VALUES) MatZeroEntries(A);
515:   for (n = 0; n < n_i; n++) MatSetValue(A, coo_i[n], coo_j[n], coo_v ? coo_v[n] : zero, ADD_VALUES);
516:   ISRestoreIndices(is_coo_i, &coo_i);
517:   ISRestoreIndices(is_coo_j, &coo_j);
518:   return 0;
519: }

521: PetscErrorCode MatSetPreallocationCOO_Basic(Mat A, PetscCount ncoo, const PetscInt coo_i[], const PetscInt coo_j[])
522: {
523:   Mat         preallocator;
524:   IS          is_coo_i, is_coo_j;
525:   PetscScalar zero = 0.0;

527:   PetscLayoutSetUp(A->rmap);
528:   PetscLayoutSetUp(A->cmap);
529:   MatCreate(PetscObjectComm((PetscObject)A), &preallocator);
530:   MatSetType(preallocator, MATPREALLOCATOR);
531:   MatSetSizes(preallocator, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N);
532:   MatSetLayouts(preallocator, A->rmap, A->cmap);
533:   MatSetUp(preallocator);
534:   for (PetscCount n = 0; n < ncoo; n++) MatSetValue(preallocator, coo_i[n], coo_j[n], zero, INSERT_VALUES);
535:   MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
536:   MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);
537:   MatPreallocatorPreallocate(preallocator, PETSC_TRUE, A);
538:   MatDestroy(&preallocator);
540:   ISCreateGeneral(PETSC_COMM_SELF, ncoo, coo_i, PETSC_COPY_VALUES, &is_coo_i);
541:   ISCreateGeneral(PETSC_COMM_SELF, ncoo, coo_j, PETSC_COPY_VALUES, &is_coo_j);
542:   PetscObjectCompose((PetscObject)A, "__PETSc_coo_i", (PetscObject)is_coo_i);
543:   PetscObjectCompose((PetscObject)A, "__PETSc_coo_j", (PetscObject)is_coo_j);
544:   ISDestroy(&is_coo_i);
545:   ISDestroy(&is_coo_j);
546:   return 0;
547: }

549: /*@C
550:    MatSetPreallocationCOO - set preallocation for matrices using a coordinate format of the entries with global indices

552:    Collective

554:    Input Parameters:
555: +  A - matrix being preallocated
556: .  ncoo - number of entries
557: .  coo_i - row indices
558: -  coo_j - column indices

560:    Level: beginner

562:    Notes:
563:    The indices coo_i and coo_j may be modified within this function. The caller should not rely on them
564:    having any specific value after this function returns. The arrays can be freed or reused immediately
565:    after this function returns.

567:    Entries can be repeated, see `MatSetValuesCOO()`. Entries with negative row or column indices are allowed
568:    but will be ignored. The corresponding entries in `MatSetValuesCOO()` will be ignored too. Remote entries
569:    are allowed and will be properly added or inserted to the matrix, unless the matrix option `MAT_IGNORE_OFF_PROC_ENTRIES`
570:    is set, in which case remote entries are ignored, or `MAT_NO_OFF_PROC_ENTRIES` is set, in which case an error will be generated.

572:    If you just want to create a sequential AIJ matrix (`MATSEQAIJ`), and your matrix entries in COO format are unique, you can also use
573:    `MatCreateSeqAIJFromTriple()`. But that is not recommended for iterative applications.

575: .seealso: `MatSetValuesCOO()`, `MatSeqAIJSetPreallocation()`, `MatMPIAIJSetPreallocation()`, `MatSeqBAIJSetPreallocation()`, `MatMPIBAIJSetPreallocation()`, `MatSeqSBAIJSetPreallocation()`, `MatMPISBAIJSetPreallocation()`, `MatSetPreallocationCOOLocal()`, `DMSetMatrixPreallocateSkip()`, `MatCreateSeqAIJFromTriple()`
576: @*/
577: PetscErrorCode MatSetPreallocationCOO(Mat A, PetscCount ncoo, PetscInt coo_i[], PetscInt coo_j[])
578: {
579:   PetscErrorCode (*f)(Mat, PetscCount, const PetscInt[], const PetscInt[]) = NULL;

585:   PetscLayoutSetUp(A->rmap);
586:   PetscLayoutSetUp(A->cmap);
587:   PetscObjectQueryFunction((PetscObject)A, "MatSetPreallocationCOO_C", &f);

589:   PetscLogEventBegin(MAT_PreallCOO, A, 0, 0, 0);
590:   if (f) {
591:     (*f)(A, ncoo, coo_i, coo_j);
592:   } else { /* allow fallback, very slow */
593:     MatSetPreallocationCOO_Basic(A, ncoo, coo_i, coo_j);
594:   }
595:   PetscLogEventEnd(MAT_PreallCOO, A, 0, 0, 0);
596:   A->preallocated = PETSC_TRUE;
597:   A->nonzerostate++;
598:   return 0;
599: }

601: /*@C
602:    MatSetPreallocationCOOLocal - set preallocation for matrices using a coordinate format of the entries with local indices

604:    Collective

606:    Input Parameters:
607: +  A - matrix being preallocated
608: .  ncoo - number of entries
609: .  coo_i - row indices (local numbering; may be modified)
610: -  coo_j - column indices (local numbering; may be modified)

612:    Level: beginner

614:    Notes:
615:    The local indices are translated using the local to global mapping, thus `MatSetLocalToGlobalMapping()` must have been
616:    called prior to this function. For matrices created with `DMCreateMatrix()` the local to global mapping is often already provided.

618:    The indices coo_i and coo_j may be modified within this function. They might be translated to corresponding global
619:    indices, but the caller should not rely on them having any specific value after this function returns. The arrays
620:    can be freed or reused immediately after this function returns.

622:    Entries can be repeated, see `MatSetValuesCOO()`. Entries with negative row or column indices are allowed
623:    but will be ignored. The corresponding entries in `MatSetValuesCOO()` will be ignored too. Remote entries
624:    are allowed and will be properly added or inserted to the matrix.

626: .seealso: `MatSetValuesCOO()`, `MatSeqAIJSetPreallocation()`, `MatMPIAIJSetPreallocation()`, `MatSeqBAIJSetPreallocation()`, `MatMPIBAIJSetPreallocation()`, `MatSeqSBAIJSetPreallocation()`, `MatMPISBAIJSetPreallocation()`, `MatSetPreallocationCOO()`, `DMSetMatrixPreallocateSkip()`
627: @*/
628: PetscErrorCode MatSetPreallocationCOOLocal(Mat A, PetscCount ncoo, PetscInt coo_i[], PetscInt coo_j[])
629: {
630:   PetscErrorCode (*f)(Mat, PetscCount, PetscInt[], PetscInt[]) = NULL;

637:   PetscLayoutSetUp(A->rmap);
638:   PetscLayoutSetUp(A->cmap);

640:   PetscObjectQueryFunction((PetscObject)A, "MatSetPreallocationCOOLocal_C", &f);
641:   if (f) {
642:     (*f)(A, ncoo, coo_i, coo_j);
643:     A->nonzerostate++;
644:   } else {
645:     ISLocalToGlobalMapping ltog_row, ltog_col;
646:     MatGetLocalToGlobalMapping(A, &ltog_row, &ltog_col);
647:     if (ltog_row) ISLocalToGlobalMappingApply(ltog_row, ncoo, coo_i, coo_i);
648:     if (ltog_col) ISLocalToGlobalMappingApply(ltog_col, ncoo, coo_j, coo_j);
649:     MatSetPreallocationCOO(A, ncoo, coo_i, coo_j);
650:   }
651:   A->preallocated = PETSC_TRUE;
652:   return 0;
653: }

655: /*@
656:    MatSetValuesCOO - set values at once in a matrix preallocated using `MatSetPreallocationCOO()`

658:    Collective

660:    Input Parameters:
661: +  A - matrix being preallocated
662: .  coo_v - the matrix values (can be NULL)
663: -  imode - the insert mode

665:    Level: beginner

667:    Notes:
668:    The values must follow the order of the indices prescribed with `MatSetPreallocationCOO()` or `MatSetPreallocationCOOLocal()`.

670:           When repeated entries are specified in the COO indices the coo_v values are first properly summed, regardless of the value of imode.
671:           The imode flag indicates if coo_v must be added to the current values of the matrix (`ADD_VALUES`) or overwritten (`INSERT_VALUES`).

673:           `MatAssemblyBegin()` and `MatAssemblyEnd()` do not need to be called after this routine. It automatically handles the assembly process.

675: .seealso: `MatSetPreallocationCOO()`, `MatSetPreallocationCOOLocal()`, `InsertMode`, `INSERT_VALUES`, `ADD_VALUES`
676: @*/
677: PetscErrorCode MatSetValuesCOO(Mat A, const PetscScalar coo_v[], InsertMode imode)
678: {
679:   PetscErrorCode (*f)(Mat, const PetscScalar[], InsertMode) = NULL;

683:   MatCheckPreallocated(A, 1);
685:   PetscObjectQueryFunction((PetscObject)A, "MatSetValuesCOO_C", &f);
686:   PetscLogEventBegin(MAT_SetVCOO, A, 0, 0, 0);
687:   if (f) {
688:     (*f)(A, coo_v, imode);
689:   } else { /* allow fallback */
690:     MatSetValuesCOO_Basic(A, coo_v, imode);
691:   }
692:   PetscLogEventEnd(MAT_SetVCOO, A, 0, 0, 0);
693:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
694:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
695:   return 0;
696: }

698: /*@
699:    MatSetBindingPropagates - Sets whether the state of being bound to the CPU for a GPU matrix type propagates to child and some other associated objects

701:    Input Parameters:
702: +  A - the matrix
703: -  flg - flag indicating whether the boundtocpu flag should be propagated

705:    Level: developer

707:    Notes:
708:    If the value of flg is set to true, the following will occur:

710:    `MatCreateSubMatrices()` and `MatCreateRedundantMatrix()` will bind created matrices to CPU if the input matrix is bound to the CPU.

712:    `MatCreateVecs()` will bind created vectors to CPU if the input matrix is bound to the CPU.
713:    The bindingpropagates flag itself is also propagated by the above routines.

715:    Developer Note:
716:    If the fine-scale `DMDA `has the -dm_bind_below option set to true, then `DMCreateInterpolationScale()` calls `MatSetBindingPropagates()`
717:    on the restriction/interpolation operator to set the bindingpropagates flag to true.

719: .seealso: `VecSetBindingPropagates()`, `MatGetBindingPropagates()`
720: @*/
721: PetscErrorCode MatSetBindingPropagates(Mat A, PetscBool flg)
722: {
724: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
725:   A->bindingpropagates = flg;
726: #endif
727:   return 0;
728: }

730: /*@
731:    MatGetBindingPropagates - Gets whether the state of being bound to the CPU for a GPU matrix type propagates to child and some other associated objects

733:    Input Parameter:
734: .  A - the matrix

736:    Output Parameter:
737: .  flg - flag indicating whether the boundtocpu flag will be propagated

739:    Level: developer

741: .seealso: `MatSetBindingPropagates()`
742: @*/
743: PetscErrorCode MatGetBindingPropagates(Mat A, PetscBool *flg)
744: {
747: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
748:   *flg = A->bindingpropagates;
749: #else
750:   *flg = PETSC_FALSE;
751: #endif
752:   return 0;
753: }