Actual source code: gcreate.c

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

  3: #include <../src/mat/impls/aij/seq/aij.h>
  4: #include <../src/mat/impls/aij/mpi/mpiaij.h>

  6: PetscErrorCode MatSetBlockSizes_Default(Mat mat, PetscInt rbs, PetscInt cbs)
  7: {
  8:   PetscFunctionBegin;
  9:   if (!mat->preallocated) PetscFunctionReturn(PETSC_SUCCESS);
 10:   PetscCheck(mat->rmap->bs <= 0 || mat->rmap->bs == rbs, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot change row block size %" PetscInt_FMT " to %" PetscInt_FMT, mat->rmap->bs, rbs);
 11:   PetscCheck(mat->cmap->bs <= 0 || mat->cmap->bs == cbs, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot change column block size %" PetscInt_FMT " to %" PetscInt_FMT, mat->cmap->bs, cbs);
 12:   PetscFunctionReturn(PETSC_SUCCESS);
 13: }

 15: PetscErrorCode MatShift_Basic(Mat Y, PetscScalar a)
 16: {
 17:   PetscInt    i, start, end, oldValA = 0, oldValB = 0;
 18:   PetscScalar alpha = a;
 19:   PetscBool   prevoption;
 20:   PetscBool   isSeqAIJDerived, isMPIAIJDerived; // all classes sharing SEQAIJHEADER or MPIAIJHEADER
 21:   Mat         A = NULL, B = NULL;

 23:   PetscFunctionBegin;
 24:   PetscCall(MatGetOption(Y, MAT_NO_OFF_PROC_ENTRIES, &prevoption));
 25:   PetscCall(MatSetOption(Y, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
 26:   PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)Y, &isSeqAIJDerived, MATSEQAIJ, MATSEQBAIJ, MATSEQSBAIJ, ""));
 27:   PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)Y, &isMPIAIJDerived, MATMPIAIJ, MATMPIBAIJ, MATMPISBAIJ, ""));

 29:   if (isSeqAIJDerived) A = Y;
 30:   else if (isMPIAIJDerived) {
 31:     Mat_MPIAIJ *mpiaij = (Mat_MPIAIJ *)Y->data;
 32:     A                  = mpiaij->A;
 33:     B                  = mpiaij->B;
 34:   }

 36:   if (A) {
 37:     oldValA                        = ((Mat_SeqAIJ *)A->data)->nonew;
 38:     ((Mat_SeqAIJ *)A->data)->nonew = 0; // so that new nonzero locations are allowed
 39:   }
 40:   if (B) {
 41:     oldValB                        = ((Mat_SeqAIJ *)B->data)->nonew;
 42:     ((Mat_SeqAIJ *)B->data)->nonew = 0;
 43:   }

 45:   PetscCall(MatGetOwnershipRange(Y, &start, &end));
 46:   for (i = start; i < end; i++) {
 47:     if (i < Y->cmap->N) PetscCall(MatSetValues(Y, 1, &i, 1, &i, &alpha, ADD_VALUES));
 48:   }
 49:   PetscCall(MatAssemblyBegin(Y, MAT_FINAL_ASSEMBLY));
 50:   PetscCall(MatAssemblyEnd(Y, MAT_FINAL_ASSEMBLY));
 51:   PetscCall(MatSetOption(Y, MAT_NO_OFF_PROC_ENTRIES, prevoption));
 52:   if (A) ((Mat_SeqAIJ *)A->data)->nonew = oldValA;
 53:   if (B) ((Mat_SeqAIJ *)B->data)->nonew = oldValB;
 54:   PetscFunctionReturn(PETSC_SUCCESS);
 55: }

 57: /*@
 58:   MatCreate - Creates a matrix where the type is determined
 59:   from either a call to `MatSetType()` or from the options database
 60:   with a call to `MatSetFromOptions()`.

 62:   Collective

 64:   Input Parameter:
 65: . comm - MPI communicator

 67:   Output Parameter:
 68: . A - the matrix

 70:   Options Database Keys:
 71: + -mat_type seqaij   - `MATSEQAIJ` type, uses `MatCreateSeqAIJ()`
 72: . -mat_type mpiaij   - `MATMPIAIJ` type, uses `MatCreateAIJ()`
 73: . -mat_type seqdense - `MATSEQDENSE`, uses `MatCreateSeqDense()`
 74: . -mat_type mpidense - `MATMPIDENSE` type, uses `MatCreateDense()`
 75: . -mat_type seqbaij  - `MATSEQBAIJ` type, uses `MatCreateSeqBAIJ()`
 76: - -mat_type mpibaij  - `MATMPIBAIJ` type, uses `MatCreateBAIJ()`

 78:    See the manpages for particular formats (e.g., `MATSEQAIJ`)
 79:    for additional format-specific options.

 81:   Level: beginner

 83:   Notes:
 84:   The default matrix type is `MATAIJ`, using the routines `MatCreateSeqAIJ()` or
 85:   `MatCreateAIJ()` if you do not set a type in the options database. If you never call
 86:   `MatSetType()` or `MatSetFromOptions()` it will generate an error when you try to use the
 87:   matrix.

 89: .seealso: [](ch_matrices), `Mat`, `MatCreateSeqAIJ()`, `MatCreateAIJ()`,
 90:           `MatCreateSeqDense()`, `MatCreateDense()`,
 91:           `MatCreateSeqBAIJ()`, `MatCreateBAIJ()`,
 92:           `MatCreateSeqSBAIJ()`, `MatCreateSBAIJ()`,
 93:           `MatConvert()`
 94: @*/
 95: PetscErrorCode MatCreate(MPI_Comm comm, Mat *A)
 96: {
 97:   Mat B;

 99:   PetscFunctionBegin;
100:   PetscAssertPointer(A, 2);
101:   PetscCall(MatInitializePackage());

103:   PetscCall(PetscHeaderCreate(B, MAT_CLASSID, "Mat", "Matrix", "Mat", comm, MatDestroy, MatView));
104:   PetscCall(PetscLayoutCreate(comm, &B->rmap));
105:   PetscCall(PetscLayoutCreate(comm, &B->cmap));
106:   PetscCall(PetscStrallocpy(VECSTANDARD, &B->defaultvectype));
107:   PetscCall(PetscStrallocpy(PETSCRANDER48, &B->defaultrandtype));

109:   B->symmetric                   = PETSC_BOOL3_UNKNOWN;
110:   B->hermitian                   = PETSC_BOOL3_UNKNOWN;
111:   B->structurally_symmetric      = PETSC_BOOL3_UNKNOWN;
112:   B->spd                         = PETSC_BOOL3_UNKNOWN;
113:   B->symmetry_eternal            = PETSC_FALSE;
114:   B->structural_symmetry_eternal = PETSC_FALSE;

116:   B->congruentlayouts = PETSC_DECIDE;
117:   B->preallocated     = PETSC_FALSE;
118: #if defined(PETSC_HAVE_DEVICE)
119:   B->boundtocpu = PETSC_TRUE;
120: #endif
121:   *A = B;
122:   PetscFunctionReturn(PETSC_SUCCESS);
123: }

125: /*@
126:   MatCreateFromOptions - Creates a matrix whose type is set from the options database

128:   Collective

130:   Input Parameters:
131: + comm   - MPI communicator
132: . prefix - [optional] prefix for the options database
133: . bs     - the blocksize (commonly 1)
134: . m      - the local number of rows (or `PETSC_DECIDE`)
135: . n      - the local number of columns (or `PETSC_DECIDE` or `PETSC_DETERMINE`)
136: . M      - the global number of rows (or `PETSC_DETERMINE`)
137: - N      - the global number of columns (or `PETSC_DETERMINE`)

139:   Output Parameter:
140: . A - the matrix

142:   Options Database Key:
143: . -mat_type - see `MatType`, for example `aij`, `aijcusparse`, `baij`, `sbaij`, dense, defaults to `aij`

145:   Level: beginner

147: .seealso: [](ch_matrices), `Mat`, `MatCreateSeqAIJ()`, `MatCreateAIJ()`,
148:           `MatCreateSeqDense()`, `MatCreateDense()`,
149:           `MatCreateSeqBAIJ()`, `MatCreateBAIJ()`,
150:           `MatCreateSeqSBAIJ()`, `MatCreateSBAIJ()`,
151:           `MatConvert()`, `MatCreate()`
152: @*/
153: PetscErrorCode MatCreateFromOptions(MPI_Comm comm, const char *prefix, PetscInt bs, PetscInt m, PetscInt n, PetscInt M, PetscInt N, Mat *A)
154: {
155:   PetscFunctionBegin;
156:   PetscAssertPointer(A, 8);
157:   PetscCall(MatCreate(comm, A));
158:   if (prefix) PetscCall(MatSetOptionsPrefix(*A, prefix));
159:   PetscCall(MatSetBlockSize(*A, bs));
160:   PetscCall(MatSetSizes(*A, m, n, M, N));
161:   PetscCall(MatSetFromOptions(*A));
162:   PetscFunctionReturn(PETSC_SUCCESS);
163: }

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

168:   Logically Collective

170:   Input Parameters:
171: + mat - matrix obtained from `MatCreate()`
172: - flg - `PETSC_TRUE` indicates you want the error generated

174:   Level: advanced

176:   Note:
177:   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
178:   or result in a `KSPConvergedReason` indicating the method did not converge.

180: .seealso: [](ch_matrices), `Mat`, `PCSetErrorIfFailure()`, `KSPConvergedReason`, `SNESConvergedReason`
181: @*/
182: PetscErrorCode MatSetErrorIfFailure(Mat mat, PetscBool flg)
183: {
184:   PetscFunctionBegin;
187:   mat->erroriffailure = flg;
188:   PetscFunctionReturn(PETSC_SUCCESS);
189: }

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

194:   Collective

196:   Input Parameters:
197: + A - the matrix
198: . m - number of local rows (or `PETSC_DECIDE`)
199: . n - number of local columns (or `PETSC_DECIDE`)
200: . M - number of global rows (or `PETSC_DETERMINE`)
201: - N - number of global columns (or `PETSC_DETERMINE`)

203:   Level: beginner

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

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

217:   If `m` and `n` are not `PETSC_DECIDE`, then the values determine the `PetscLayout` of the matrix and the ranges returned by
218:   `MatGetOwnershipRange()`,  `MatGetOwnershipRanges()`, `MatGetOwnershipRangeColumn()`, and `MatGetOwnershipRangesColumn()`.

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

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

224: .seealso: [](ch_matrices), `Mat`, `MatGetSize()`, `PetscSplitOwnership()`, `MatGetOwnershipRange()`, `MatGetOwnershipRanges()`,
225:           `MatGetOwnershipRangeColumn()`, `MatGetOwnershipRangesColumn()`, `PetscLayout`, `VecSetSizes()`
226: @*/
227: PetscErrorCode MatSetSizes(Mat A, PetscInt m, PetscInt n, PetscInt M, PetscInt N)
228: {
229:   PetscFunctionBegin;
233:   PetscCheck(M <= 0 || m <= M, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local row size %" PetscInt_FMT " cannot be larger than global row size %" PetscInt_FMT, m, M);
234:   PetscCheck(N <= 0 || n <= N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local column size %" PetscInt_FMT " cannot be larger than global column size %" PetscInt_FMT, n, N);
235:   PetscCheck((A->rmap->n < 0 || A->rmap->N < 0) || (A->rmap->n == m && (M <= 0 || A->rmap->N == M)), PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot change/reset row sizes to %" PetscInt_FMT " local %" PetscInt_FMT " global after previously setting them to %" PetscInt_FMT " local %" PetscInt_FMT " global", m, M,
236:              A->rmap->n, A->rmap->N);
237:   PetscCheck((A->cmap->n < 0 || A->cmap->N < 0) || (A->cmap->n == n && (N <= 0 || A->cmap->N == N)), PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot change/reset column sizes to %" PetscInt_FMT " local %" PetscInt_FMT " global after previously setting them to %" PetscInt_FMT " local %" PetscInt_FMT " global", n, N,
238:              A->cmap->n, A->cmap->N);
239:   A->rmap->n = m;
240:   A->cmap->n = n;
241:   A->rmap->N = M > -1 ? M : A->rmap->N;
242:   A->cmap->N = N > -1 ? N : A->cmap->N;
243:   PetscFunctionReturn(PETSC_SUCCESS);
244: }

246: /*@
247:   MatSetFromOptions - Creates a matrix where the type is determined
248:   from the options database.

250:   Collective

252:   Input Parameter:
253: . B - the matrix

255:   Options Database Keys:
256: + -mat_type seqaij   - `MATSEQAIJ` type, uses `MatCreateSeqAIJ()`
257: . -mat_type mpiaij   - `MATMPIAIJ` type, uses `MatCreateAIJ()`
258: . -mat_type seqdense - `MATSEQDENSE` type, uses `MatCreateSeqDense()`
259: . -mat_type mpidense - `MATMPIDENSE`, uses `MatCreateDense()`
260: . -mat_type seqbaij  - `MATSEQBAIJ`, uses `MatCreateSeqBAIJ()`
261: - -mat_type mpibaij  - `MATMPIBAIJ`, uses `MatCreateBAIJ()`

263:    See the manpages for particular formats (e.g., `MATSEQAIJ`)
264:    for additional format-specific options.

266:   Level: beginner

268:   Notes:
269:   Generates a parallel MPI matrix if the communicator has more than one processor.  The default
270:   matrix type is `MATAIJ`, using the routines `MatCreateSeqAIJ()` and `MatCreateAIJ()` if you
271:   do not select a type in the options database.

273: .seealso: [](ch_matrices), `Mat`, `MatCreateSeqAIJ()`, `MatCreateAIJ()`,
274:           `MatCreateSeqDense()`, `MatCreateDense()`,
275:           `MatCreateSeqBAIJ()`, `MatCreateBAIJ()`,
276:           `MatCreateSeqSBAIJ()`, `MatCreateSBAIJ()`,
277:           `MatConvert()`
278: @*/
279: PetscErrorCode MatSetFromOptions(Mat B)
280: {
281:   const char *deft = MATAIJ;
282:   char        type[256];
283:   PetscBool   flg, set;
284:   PetscInt    bind_below = 0;

286:   PetscFunctionBegin;

289:   PetscObjectOptionsBegin((PetscObject)B);

291:   if (B->rmap->bs < 0) {
292:     PetscInt newbs = -1;
293:     PetscCall(PetscOptionsInt("-mat_block_size", "Set the blocksize used to store the matrix", "MatSetBlockSize", newbs, &newbs, &flg));
294:     if (flg) {
295:       PetscCall(PetscLayoutSetBlockSize(B->rmap, newbs));
296:       PetscCall(PetscLayoutSetBlockSize(B->cmap, newbs));
297:     }
298:   }

300:   PetscCall(PetscOptionsFList("-mat_type", "Matrix type", "MatSetType", MatList, deft, type, 256, &flg));
301:   if (flg) {
302:     PetscCall(MatSetType(B, type));
303:   } else if (!((PetscObject)B)->type_name) {
304:     PetscCall(MatSetType(B, deft));
305:   }

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

312:   PetscTryTypeMethod(B, setfromoptions, PetscOptionsObject);

314:   flg = PETSC_FALSE;
315:   PetscCall(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));
316:   if (set) PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, flg));
317:   flg = PETSC_FALSE;
318:   PetscCall(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));
319:   if (set) PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, flg));
320:   flg = PETSC_FALSE;
321:   PetscCall(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));
322:   if (set) PetscCall(MatSetOption(B, MAT_IGNORE_ZERO_ENTRIES, flg));

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

328:   /* Bind to CPU if below a user-specified size threshold.
329:    * This perhaps belongs in the options for the GPU Mat types, but MatBindToCPU() does nothing when called on non-GPU types,
330:    * and putting it here makes is more maintainable than duplicating this for all. */
331:   PetscCall(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));
332:   if (flg && B->rmap->n < bind_below) PetscCall(MatBindToCPU(B, PETSC_TRUE));

334:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
335:   PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)B, PetscOptionsObject));
336:   PetscOptionsEnd();
337:   PetscFunctionReturn(PETSC_SUCCESS);
338: }

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

343:   Collective

345:   Input Parameters:
346: + A     - matrix being preallocated
347: . bs    - block size
348: . dnnz  - number of nonzero column blocks per block row of diagonal part of parallel matrix
349: . onnz  - number of nonzero column blocks per block row of off-diagonal part of parallel matrix
350: . dnnzu - number of nonzero column blocks per block row of upper-triangular part of diagonal part of parallel matrix
351: - onnzu - number of nonzero column blocks per block row of upper-triangular part of off-diagonal part of parallel matrix

353:   Level: beginner

355: .seealso: [](ch_matrices), `Mat`, `MatSeqAIJSetPreallocation()`, `MatMPIAIJSetPreallocation()`, `MatSeqBAIJSetPreallocation()`, `MatMPIBAIJSetPreallocation()`,
356:           `MatSeqSBAIJSetPreallocation()`, `MatMPISBAIJSetPreallocation()`,
357:           `PetscSplitOwnership()`
358: @*/
359: PetscErrorCode MatXAIJSetPreallocation(Mat A, PetscInt bs, const PetscInt dnnz[], const PetscInt onnz[], const PetscInt dnnzu[], const PetscInt onnzu[])
360: {
361:   PetscInt cbs;
362:   void (*aij)(void);
363:   void (*is)(void);
364:   void (*hyp)(void) = NULL;

366:   PetscFunctionBegin;
367:   if (bs != PETSC_DECIDE) { /* don't mess with an already set block size */
368:     PetscCall(MatSetBlockSize(A, bs));
369:   }
370:   PetscCall(PetscLayoutSetUp(A->rmap));
371:   PetscCall(PetscLayoutSetUp(A->cmap));
372:   PetscCall(MatGetBlockSizes(A, &bs, &cbs));
373:   /* these routines assumes bs == cbs, this should be checked somehow */
374:   PetscCall(MatSeqBAIJSetPreallocation(A, bs, 0, dnnz));
375:   PetscCall(MatMPIBAIJSetPreallocation(A, bs, 0, dnnz, 0, onnz));
376:   PetscCall(MatSeqSBAIJSetPreallocation(A, bs, 0, dnnzu));
377:   PetscCall(MatMPISBAIJSetPreallocation(A, bs, 0, dnnzu, 0, onnzu));
378:   /*
379:     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
380:     good before going on with it.
381:   */
382:   PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatMPIAIJSetPreallocation_C", &aij));
383:   PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatISSetPreallocation_C", &is));
384: #if defined(PETSC_HAVE_HYPRE)
385:   PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatHYPRESetPreallocation_C", &hyp));
386: #endif
387:   if (!aij && !is && !hyp) PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatSeqAIJSetPreallocation_C", &aij));
388:   if (aij || is || hyp) {
389:     if (bs == cbs && bs == 1) {
390:       PetscCall(MatSeqAIJSetPreallocation(A, 0, dnnz));
391:       PetscCall(MatMPIAIJSetPreallocation(A, 0, dnnz, 0, onnz));
392:       PetscCall(MatISSetPreallocation(A, 0, dnnz, 0, onnz));
393: #if defined(PETSC_HAVE_HYPRE)
394:       PetscCall(MatHYPRESetPreallocation(A, 0, dnnz, 0, onnz));
395: #endif
396:     } else { /* Convert block-row precallocation to scalar-row */
397:       PetscInt i, m, *sdnnz, *sonnz;
398:       PetscCall(MatGetLocalSize(A, &m, NULL));
399:       PetscCall(PetscMalloc2((!!dnnz) * m, &sdnnz, (!!onnz) * m, &sonnz));
400:       for (i = 0; i < m; i++) {
401:         if (dnnz) sdnnz[i] = dnnz[i / bs] * cbs;
402:         if (onnz) sonnz[i] = onnz[i / bs] * cbs;
403:       }
404:       PetscCall(MatSeqAIJSetPreallocation(A, 0, dnnz ? sdnnz : NULL));
405:       PetscCall(MatMPIAIJSetPreallocation(A, 0, dnnz ? sdnnz : NULL, 0, onnz ? sonnz : NULL));
406:       PetscCall(MatISSetPreallocation(A, 0, dnnz ? sdnnz : NULL, 0, onnz ? sonnz : NULL));
407: #if defined(PETSC_HAVE_HYPRE)
408:       PetscCall(MatHYPRESetPreallocation(A, 0, dnnz ? sdnnz : NULL, 0, onnz ? sonnz : NULL));
409: #endif
410:       PetscCall(PetscFree2(sdnnz, sonnz));
411:     }
412:   }
413:   PetscFunctionReturn(PETSC_SUCCESS);
414: }

416: /*@C
417:   MatHeaderMerge - Merges some information from the header of `C` to `A`; the `C` object is then destroyed

419:   Collective, No Fortran Support

421:   Input Parameters:
422: + A - a `Mat` being merged into
423: - C - the `Mat` providing the merge information

425:   Level: developer

427:   Notes:
428:   `A` and `C` must be of the same type.
429:   The object list and query function list in `A` are retained, as well as the object name, and prefix.
430:   The object state of `A` is increased by 1.

432:   Developer Note:
433:   This is somewhat different from `MatHeaderReplace()`, it would be nice to merge the code

435: .seealso: `Mat`, `MatHeaderReplace()`
436:  @*/
437: PetscErrorCode MatHeaderMerge(Mat A, Mat *C)
438: {
439:   PetscInt          refct;
440:   PetscOps          Abops;
441:   struct _MatOps    Aops;
442:   char             *mtype, *mname, *mprefix;
443:   Mat_Product      *product;
444:   Mat_Redundant    *redundant;
445:   PetscObjectState  state;
446:   PetscObjectList   olist;
447:   PetscFunctionList qlist;

449:   PetscFunctionBegin;
452:   if (A == *C) PetscFunctionReturn(PETSC_SUCCESS);
453:   PetscCheckSameTypeAndComm(A, 1, *C, 2);
454:   /* save the parts of A we need */
455:   Abops     = ((PetscObject)A)->bops[0];
456:   Aops      = A->ops[0];
457:   refct     = ((PetscObject)A)->refct;
458:   mtype     = ((PetscObject)A)->type_name;
459:   mname     = ((PetscObject)A)->name;
460:   state     = ((PetscObject)A)->state;
461:   mprefix   = ((PetscObject)A)->prefix;
462:   product   = A->product;
463:   redundant = A->redundant;
464:   qlist     = ((PetscObject)A)->qlist;
465:   olist     = ((PetscObject)A)->olist;

467:   /* zero these so the destroy below does not free them */
468:   ((PetscObject)A)->type_name = NULL;
469:   ((PetscObject)A)->name      = NULL;
470:   ((PetscObject)A)->qlist     = NULL;
471:   ((PetscObject)A)->olist     = NULL;

473:   /*
474:      free all the interior data structures from mat
475:      cannot use PetscUseTypeMethod(A,destroy); because compiler
476:      thinks it may print NULL type_name and name
477:   */
478:   PetscTryTypeMethod(A, destroy);

480:   PetscCall(PetscFree(A->defaultvectype));
481:   PetscCall(PetscFree(A->defaultrandtype));
482:   PetscCall(PetscLayoutDestroy(&A->rmap));
483:   PetscCall(PetscLayoutDestroy(&A->cmap));
484:   PetscCall(PetscComposedQuantitiesDestroy((PetscObject)A));

486:   /* copy C over to A */
487:   PetscCall(PetscFree(A->factorprefix));
488:   PetscCall(PetscMemcpy(A, *C, sizeof(struct _p_Mat)));

490:   /* return the parts of A we saved */
491:   ((PetscObject)A)->bops[0]   = Abops;
492:   A->ops[0]                   = Aops;
493:   ((PetscObject)A)->refct     = refct;
494:   ((PetscObject)A)->type_name = mtype;
495:   ((PetscObject)A)->name      = mname;
496:   ((PetscObject)A)->prefix    = mprefix;
497:   ((PetscObject)A)->state     = state + 1;
498:   A->product                  = product;
499:   A->redundant                = redundant;

501:   /* Append the saved lists */
502:   PetscCall(PetscFunctionListDuplicate(qlist, &((PetscObject)A)->qlist));
503:   PetscCall(PetscObjectListDuplicate(olist, &((PetscObject)A)->olist));
504:   PetscCall(PetscFunctionListDestroy(&qlist));
505:   PetscCall(PetscObjectListDestroy(&olist));

507:   /* since these two are copied into A we do not want them destroyed in C */
508:   ((PetscObject)*C)->qlist = NULL;
509:   ((PetscObject)*C)->olist = NULL;
510:   PetscCall(PetscHeaderDestroy(C));
511:   PetscFunctionReturn(PETSC_SUCCESS);
512: }

514: /*@
515:   MatHeaderReplace - Replaces the internal data of matrix `A` by the internal data of matrix `C` while deleting the outer wrapper of `C`

517:   Input Parameters:
518: + A - a `Mat` whose internal data is to be replaced
519: - C - the `Mat` providing new internal data for `A`

521:   Level: advanced

523:   Example Usage\:
524: .vb
525:   Mat C;
526:   MatCreateSeqAIJWithArrays(..., &C);
527:   MatHeaderReplace(A, &C);
528:   // C has been destroyed and A contains the matrix entries of C
529: .ve

531:   Note:
532:   This can be used inside a function provided to `SNESSetJacobian()`, `TSSetRHSJacobian()`, or `TSSetIJacobian()` in cases where the user code computes an entirely new sparse matrix
533:   (generally with a different nonzero pattern) for each Newton update. It is usually better to reuse the matrix structure of `A` instead of constructing an entirely new one.

535:   Developer Note:
536:   This is somewhat different from `MatHeaderMerge()` it would be nice to merge the code

538: .seealso: `Mat`, `MatHeaderMerge()`
539:  @*/
540: PetscErrorCode MatHeaderReplace(Mat A, Mat *C)
541: {
542:   PetscInt         refct;
543:   PetscObjectState state;
544:   struct _p_Mat    buffer;
545:   MatStencilInfo   stencil;

547:   PetscFunctionBegin;
550:   if (A == *C) PetscFunctionReturn(PETSC_SUCCESS);
551:   PetscCheckSameComm(A, 1, *C, 2);
552:   PetscCheck(((PetscObject)*C)->refct == 1, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "Object C has refct %" PetscInt_FMT " > 1, would leave hanging reference", ((PetscObject)*C)->refct);

554:   /* swap C and A */
555:   refct   = ((PetscObject)A)->refct;
556:   state   = ((PetscObject)A)->state;
557:   stencil = A->stencil;
558:   PetscCall(PetscMemcpy(&buffer, A, sizeof(struct _p_Mat)));
559:   PetscCall(PetscMemcpy(A, *C, sizeof(struct _p_Mat)));
560:   PetscCall(PetscMemcpy(*C, &buffer, sizeof(struct _p_Mat)));
561:   ((PetscObject)A)->refct = refct;
562:   ((PetscObject)A)->state = state + 1;
563:   A->stencil              = stencil;

565:   ((PetscObject)*C)->refct = 1;
566:   PetscCall(MatDestroy(C));
567:   PetscFunctionReturn(PETSC_SUCCESS);
568: }

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

573:   Logically Collective

575:   Input Parameters:
576: + A   - the matrix
577: - flg - bind to the CPU if value of `PETSC_TRUE`

579:   Level: intermediate

581: .seealso: [](ch_matrices), `Mat`, `MatBoundToCPU()`
582: @*/
583: PetscErrorCode MatBindToCPU(Mat A, PetscBool flg)
584: {
585:   PetscFunctionBegin;
588: #if defined(PETSC_HAVE_DEVICE)
589:   if (A->boundtocpu == flg) PetscFunctionReturn(PETSC_SUCCESS);
590:   A->boundtocpu = flg;
591:   PetscTryTypeMethod(A, bindtocpu, flg);
592: #endif
593:   PetscFunctionReturn(PETSC_SUCCESS);
594: }

596: /*@
597:   MatBoundToCPU - query if a matrix is bound to the CPU

599:   Input Parameter:
600: . A - the matrix

602:   Output Parameter:
603: . flg - the logical flag

605:   Level: intermediate

607: .seealso: [](ch_matrices), `Mat`, `MatBindToCPU()`
608: @*/
609: PetscErrorCode MatBoundToCPU(Mat A, PetscBool *flg)
610: {
611:   PetscFunctionBegin;
613:   PetscAssertPointer(flg, 2);
614: #if defined(PETSC_HAVE_DEVICE)
615:   *flg = A->boundtocpu;
616: #else
617:   *flg = PETSC_TRUE;
618: #endif
619:   PetscFunctionReturn(PETSC_SUCCESS);
620: }

622: PetscErrorCode MatSetValuesCOO_Basic(Mat A, const PetscScalar coo_v[], InsertMode imode)
623: {
624:   IS              is_coo_i, is_coo_j;
625:   const PetscInt *coo_i, *coo_j;
626:   PetscInt        n, n_i, n_j;
627:   PetscScalar     zero = 0.;

629:   PetscFunctionBegin;
630:   PetscCall(PetscObjectQuery((PetscObject)A, "__PETSc_coo_i", (PetscObject *)&is_coo_i));
631:   PetscCall(PetscObjectQuery((PetscObject)A, "__PETSc_coo_j", (PetscObject *)&is_coo_j));
632:   PetscCheck(is_coo_i, PetscObjectComm((PetscObject)A), PETSC_ERR_COR, "Missing coo_i IS");
633:   PetscCheck(is_coo_j, PetscObjectComm((PetscObject)A), PETSC_ERR_COR, "Missing coo_j IS");
634:   PetscCall(ISGetLocalSize(is_coo_i, &n_i));
635:   PetscCall(ISGetLocalSize(is_coo_j, &n_j));
636:   PetscCheck(n_i == n_j, PETSC_COMM_SELF, PETSC_ERR_COR, "Wrong local size %" PetscInt_FMT " != %" PetscInt_FMT, n_i, n_j);
637:   PetscCall(ISGetIndices(is_coo_i, &coo_i));
638:   PetscCall(ISGetIndices(is_coo_j, &coo_j));
639:   if (imode != ADD_VALUES) PetscCall(MatZeroEntries(A));
640:   for (n = 0; n < n_i; n++) PetscCall(MatSetValue(A, coo_i[n], coo_j[n], coo_v ? coo_v[n] : zero, ADD_VALUES));
641:   PetscCall(ISRestoreIndices(is_coo_i, &coo_i));
642:   PetscCall(ISRestoreIndices(is_coo_j, &coo_j));
643:   PetscFunctionReturn(PETSC_SUCCESS);
644: }

646: PetscErrorCode MatSetPreallocationCOO_Basic(Mat A, PetscCount ncoo, PetscInt coo_i[], PetscInt coo_j[])
647: {
648:   Mat         preallocator;
649:   IS          is_coo_i, is_coo_j;
650:   PetscScalar zero = 0.0;

652:   PetscFunctionBegin;
653:   PetscCheck(ncoo <= PETSC_INT_MAX, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ncoo %" PetscCount_FMT " overflowed PetscInt; configure --with-64-bit-indices or request support", ncoo);
654:   PetscCall(PetscLayoutSetUp(A->rmap));
655:   PetscCall(PetscLayoutSetUp(A->cmap));
656:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &preallocator));
657:   PetscCall(MatSetType(preallocator, MATPREALLOCATOR));
658:   PetscCall(MatSetSizes(preallocator, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
659:   PetscCall(MatSetLayouts(preallocator, A->rmap, A->cmap));
660:   PetscCall(MatSetUp(preallocator));
661:   for (PetscCount n = 0; n < ncoo; n++) PetscCall(MatSetValue(preallocator, coo_i[n], coo_j[n], zero, INSERT_VALUES));
662:   PetscCall(MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY));
663:   PetscCall(MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY));
664:   PetscCall(MatPreallocatorPreallocate(preallocator, PETSC_TRUE, A));
665:   PetscCall(MatDestroy(&preallocator));
666:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, (PetscInt)ncoo, coo_i, PETSC_COPY_VALUES, &is_coo_i));
667:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, (PetscInt)ncoo, coo_j, PETSC_COPY_VALUES, &is_coo_j));
668:   PetscCall(PetscObjectCompose((PetscObject)A, "__PETSc_coo_i", (PetscObject)is_coo_i));
669:   PetscCall(PetscObjectCompose((PetscObject)A, "__PETSc_coo_j", (PetscObject)is_coo_j));
670:   PetscCall(ISDestroy(&is_coo_i));
671:   PetscCall(ISDestroy(&is_coo_j));
672:   PetscFunctionReturn(PETSC_SUCCESS);
673: }

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

678:   Collective

680:   Input Parameters:
681: + A     - matrix being preallocated
682: . ncoo  - number of entries
683: . coo_i - row indices
684: - coo_j - column indices

686:   Level: beginner

688:   Notes:
689:   The indices within `coo_i` and `coo_j` may be modified within this function. The caller should not rely on them
690:   having any specific value after this function returns. The arrays can be freed or reused immediately
691:   after this function returns.

693:   Entries can be repeated, see `MatSetValuesCOO()`. Entries with negative row or column indices are allowed
694:   but will be ignored. The corresponding entries in `MatSetValuesCOO()` will be ignored too. Remote entries
695:   are allowed and will be properly added or inserted to the matrix, unless the matrix option `MAT_IGNORE_OFF_PROC_ENTRIES`
696:   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.

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

701: .seealso: [](ch_matrices), `Mat`, `MatSetValuesCOO()`, `MatSeqAIJSetPreallocation()`, `MatMPIAIJSetPreallocation()`, `MatSeqBAIJSetPreallocation()`,
702:           `MatMPIBAIJSetPreallocation()`, `MatSeqSBAIJSetPreallocation()`, `MatMPISBAIJSetPreallocation()`, `MatSetPreallocationCOOLocal()`,
703:           `DMSetMatrixPreallocateSkip()`, `MatCreateSeqAIJFromTriple()`
704: @*/
705: PetscErrorCode MatSetPreallocationCOO(Mat A, PetscCount ncoo, PetscInt coo_i[], PetscInt coo_j[])
706: {
707:   PetscErrorCode (*f)(Mat, PetscCount, PetscInt[], PetscInt[]) = NULL;

709:   PetscFunctionBegin;
712:   if (ncoo) PetscAssertPointer(coo_i, 3);
713:   if (ncoo) PetscAssertPointer(coo_j, 4);
714:   PetscCall(PetscLayoutSetUp(A->rmap));
715:   PetscCall(PetscLayoutSetUp(A->cmap));
716:   PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatSetPreallocationCOO_C", &f));

718:   PetscCall(PetscLogEventBegin(MAT_PreallCOO, A, 0, 0, 0));
719:   if (f) {
720:     PetscCall((*f)(A, ncoo, coo_i, coo_j));
721:   } else { /* allow fallback, very slow */
722:     PetscCall(MatSetPreallocationCOO_Basic(A, ncoo, coo_i, coo_j));
723:   }
724:   PetscCall(PetscLogEventEnd(MAT_PreallCOO, A, 0, 0, 0));
725:   A->preallocated = PETSC_TRUE;
726:   A->nonzerostate++;
727:   PetscFunctionReturn(PETSC_SUCCESS);
728: }

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

733:   Collective

735:   Input Parameters:
736: + A     - matrix being preallocated
737: . ncoo  - number of entries
738: . coo_i - row indices (local numbering; may be modified)
739: - coo_j - column indices (local numbering; may be modified)

741:   Level: beginner

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

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

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

755: .seealso: [](ch_matrices), `Mat`, `MatSetValuesCOO()`, `MatSeqAIJSetPreallocation()`, `MatMPIAIJSetPreallocation()`, `MatSeqBAIJSetPreallocation()`,
756:           `MatMPIBAIJSetPreallocation()`, `MatSeqSBAIJSetPreallocation()`, `MatMPISBAIJSetPreallocation()`, `MatSetPreallocationCOO()`,
757:           `DMSetMatrixPreallocateSkip()`
758: @*/
759: PetscErrorCode MatSetPreallocationCOOLocal(Mat A, PetscCount ncoo, PetscInt coo_i[], PetscInt coo_j[])
760: {
761:   PetscErrorCode (*f)(Mat, PetscCount, PetscInt[], PetscInt[]) = NULL;

763:   PetscFunctionBegin;
766:   if (ncoo) PetscAssertPointer(coo_i, 3);
767:   if (ncoo) PetscAssertPointer(coo_j, 4);
768:   PetscCheck(ncoo <= PETSC_INT_MAX, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ncoo %" PetscCount_FMT " overflowed PetscInt; configure --with-64-bit-indices or request support", ncoo);
769:   PetscCall(PetscLayoutSetUp(A->rmap));
770:   PetscCall(PetscLayoutSetUp(A->cmap));

772:   PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatSetPreallocationCOOLocal_C", &f));
773:   if (f) {
774:     PetscCall((*f)(A, ncoo, coo_i, coo_j));
775:     A->nonzerostate++;
776:   } else {
777:     ISLocalToGlobalMapping ltog_row, ltog_col;

779:     PetscCall(MatGetLocalToGlobalMapping(A, &ltog_row, &ltog_col));
780:     if (ltog_row) PetscCall(ISLocalToGlobalMappingApply(ltog_row, (PetscInt)ncoo, coo_i, coo_i));
781:     if (ltog_col) PetscCall(ISLocalToGlobalMappingApply(ltog_col, (PetscInt)ncoo, coo_j, coo_j));
782:     PetscCall(MatSetPreallocationCOO(A, ncoo, coo_i, coo_j));
783:   }
784:   A->preallocated = PETSC_TRUE;
785:   PetscFunctionReturn(PETSC_SUCCESS);
786: }

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

791:   Collective

793:   Input Parameters:
794: + A     - matrix being preallocated
795: . coo_v - the matrix values (can be `NULL`)
796: - imode - the insert mode

798:   Level: beginner

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

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

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

808: .seealso: [](ch_matrices), `Mat`, `MatSetPreallocationCOO()`, `MatSetPreallocationCOOLocal()`, `InsertMode`, `INSERT_VALUES`, `ADD_VALUES`
809: @*/
810: PetscErrorCode MatSetValuesCOO(Mat A, const PetscScalar coo_v[], InsertMode imode)
811: {
812:   PetscErrorCode (*f)(Mat, const PetscScalar[], InsertMode) = NULL;
813:   PetscBool oldFlg;

815:   PetscFunctionBegin;
818:   MatCheckPreallocated(A, 1);
820:   PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatSetValuesCOO_C", &f));
821:   PetscCall(PetscLogEventBegin(MAT_SetVCOO, A, 0, 0, 0));
822:   if (f) {
823:     PetscCall((*f)(A, coo_v, imode)); // all known COO implementations do not use MatStash. They do their own off-proc communication
824:     PetscCall(MatGetOption(A, MAT_NO_OFF_PROC_ENTRIES, &oldFlg));
825:     PetscCall(MatSetOption(A, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE)); // set A->nooffprocentries to avoid costly MatStash scatter in MatAssembly
826:   } else {
827:     PetscCall(MatSetValuesCOO_Basic(A, coo_v, imode)); // fall back to MatSetValues, which might use MatStash
828:   }
829:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
830:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
831:   if (f) PetscCall(MatSetOption(A, MAT_NO_OFF_PROC_ENTRIES, oldFlg));
832:   PetscCall(PetscLogEventEnd(MAT_SetVCOO, A, 0, 0, 0));
833:   PetscFunctionReturn(PETSC_SUCCESS);
834: }

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

839:   Input Parameters:
840: + A   - the matrix
841: - flg - flag indicating whether the boundtocpu flag should be propagated

843:   Level: developer

845:   Notes:
846:   If the value of flg is set to true, the following will occur
847: +   `MatCreateSubMatrices()` and `MatCreateRedundantMatrix()` - bind created matrices to CPU if the input matrix is bound to the CPU.
848: -   `MatCreateVecs()` - bind created vectors to CPU if the input matrix is bound to the CPU.

850:   The bindingpropagates flag itself is also propagated by the above routines.

852:   Developer Notes:
853:   If the fine-scale `DMDA` has the `-dm_bind_below` option set to true, then `DMCreateInterpolationScale()` calls `MatSetBindingPropagates()`
854:   on the restriction/interpolation operator to set the bindingpropagates flag to true.

856: .seealso: [](ch_matrices), `Mat`, `VecSetBindingPropagates()`, `MatGetBindingPropagates()`
857: @*/
858: PetscErrorCode MatSetBindingPropagates(Mat A, PetscBool flg)
859: {
860:   PetscFunctionBegin;
862: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
863:   A->bindingpropagates = flg;
864: #endif
865:   PetscFunctionReturn(PETSC_SUCCESS);
866: }

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

871:   Input Parameter:
872: . A - the matrix

874:   Output Parameter:
875: . flg - flag indicating whether the boundtocpu flag will be propagated

877:   Level: developer

879: .seealso: [](ch_matrices), `Mat`, `MatSetBindingPropagates()`
880: @*/
881: PetscErrorCode MatGetBindingPropagates(Mat A, PetscBool *flg)
882: {
883:   PetscFunctionBegin;
885:   PetscAssertPointer(flg, 2);
886: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
887:   *flg = A->bindingpropagates;
888: #else
889:   *flg = PETSC_FALSE;
890: #endif
891:   PetscFunctionReturn(PETSC_SUCCESS);
892: }