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: PetscFunctionBegin;
6: if (!mat->preallocated) PetscFunctionReturn(PETSC_SUCCESS);
7: 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);
8: 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);
9: PetscFunctionReturn(PETSC_SUCCESS);
10: }
12: PETSC_INTERN PetscErrorCode MatShift_Basic(Mat Y, PetscScalar a)
13: {
14: PetscInt i, start, end;
15: PetscScalar alpha = a;
16: PetscBool prevoption;
18: PetscFunctionBegin;
19: PetscCall(MatGetOption(Y, MAT_NO_OFF_PROC_ENTRIES, &prevoption));
20: PetscCall(MatSetOption(Y, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
21: PetscCall(MatGetOwnershipRange(Y, &start, &end));
22: for (i = start; i < end; i++) {
23: if (i < Y->cmap->N) PetscCall(MatSetValues(Y, 1, &i, 1, &i, &alpha, ADD_VALUES));
24: }
25: PetscCall(MatAssemblyBegin(Y, MAT_FINAL_ASSEMBLY));
26: PetscCall(MatAssemblyEnd(Y, MAT_FINAL_ASSEMBLY));
27: PetscCall(MatSetOption(Y, MAT_NO_OFF_PROC_ENTRIES, prevoption));
28: PetscFunctionReturn(PETSC_SUCCESS);
29: }
31: /*@
32: MatCreate - Creates a matrix where the type is determined
33: from either a call to `MatSetType()` or from the options database
34: with a call to `MatSetFromOptions()`.
36: Collective
38: Input Parameter:
39: . comm - MPI communicator
41: Output Parameter:
42: . A - the matrix
44: Options Database Keys:
45: + -mat_type seqaij - `MATSEQAIJ` type, uses `MatCreateSeqAIJ()`
46: . -mat_type mpiaij - `MATMPIAIJ` type, uses `MatCreateAIJ()`
47: . -mat_type seqdense - `MATSEQDENSE`, uses `MatCreateSeqDense()`
48: . -mat_type mpidense - `MATMPIDENSE` type, uses `MatCreateDense()`
49: . -mat_type seqbaij - `MATSEQBAIJ` type, uses `MatCreateSeqBAIJ()`
50: - -mat_type mpibaij - `MATMPIBAIJ` type, uses `MatCreateBAIJ()`
52: See the manpages for particular formats (e.g., `MATSEQAIJ`)
53: for additional format-specific options.
55: Level: beginner
57: Notes:
58: The default matrix type is `MATAIJ`, using the routines `MatCreateSeqAIJ()` or
59: `MatCreateAIJ()` if you do not set a type in the options database. If you never call
60: `MatSetType()` or `MatSetFromOptions()` it will generate an error when you try to use the
61: matrix.
63: .seealso: [](ch_matrices), `Mat`, `MatCreateSeqAIJ()`, `MatCreateAIJ()`,
64: `MatCreateSeqDense()`, `MatCreateDense()`,
65: `MatCreateSeqBAIJ()`, `MatCreateBAIJ()`,
66: `MatCreateSeqSBAIJ()`, `MatCreateSBAIJ()`,
67: `MatConvert()`
68: @*/
69: PetscErrorCode MatCreate(MPI_Comm comm, Mat *A)
70: {
71: Mat B;
73: PetscFunctionBegin;
74: PetscAssertPointer(A, 2);
76: *A = NULL;
77: PetscCall(MatInitializePackage());
79: PetscCall(PetscHeaderCreate(B, MAT_CLASSID, "Mat", "Matrix", "Mat", comm, MatDestroy, MatView));
80: PetscCall(PetscLayoutCreate(comm, &B->rmap));
81: PetscCall(PetscLayoutCreate(comm, &B->cmap));
82: PetscCall(PetscStrallocpy(VECSTANDARD, &B->defaultvectype));
83: PetscCall(PetscStrallocpy(PETSCRANDER48, &B->defaultrandtype));
85: B->symmetric = PETSC_BOOL3_UNKNOWN;
86: B->hermitian = PETSC_BOOL3_UNKNOWN;
87: B->structurally_symmetric = PETSC_BOOL3_UNKNOWN;
88: B->spd = PETSC_BOOL3_UNKNOWN;
89: B->symmetry_eternal = PETSC_FALSE;
90: B->structural_symmetry_eternal = PETSC_FALSE;
92: B->congruentlayouts = PETSC_DECIDE;
93: B->preallocated = PETSC_FALSE;
94: #if defined(PETSC_HAVE_DEVICE)
95: B->boundtocpu = PETSC_TRUE;
96: #endif
97: *A = B;
98: PetscFunctionReturn(PETSC_SUCCESS);
99: }
101: /*@C
102: MatCreateFromOptions - Creates a matrix whose type is set from the options database
104: Collective
106: Input Parameters:
107: + comm - MPI communicator
108: . prefix - [optional] prefix for the options database
109: . bs - the blocksize (commonly 1)
110: . m - the local number of rows (or `PETSC_DECIDE`)
111: . n - the local number of columns (or `PETSC_DECIDE` or `PETSC_DETERMINE`)
112: . M - the global number of rows (or `PETSC_DETERMINE`)
113: - N - the global number of columns (or `PETSC_DETERMINE`)
115: Output Parameter:
116: . A - the matrix
118: Options Database Key:
119: . -mat_type - see `MatType`, for example `aij`, `aijcusparse`, `baij`, `sbaij`, dense, defaults to `aij`
121: Level: beginner
123: .seealso: [](ch_matrices), `Mat`, `MatCreateSeqAIJ()`, `MatCreateAIJ()`,
124: `MatCreateSeqDense()`, `MatCreateDense()`,
125: `MatCreateSeqBAIJ()`, `MatCreateBAIJ()`,
126: `MatCreateSeqSBAIJ()`, `MatCreateSBAIJ()`,
127: `MatConvert()`, `MatCreate()`
128: @*/
129: PetscErrorCode MatCreateFromOptions(MPI_Comm comm, const char *prefix, PetscInt bs, PetscInt m, PetscInt n, PetscInt M, PetscInt N, Mat *A)
130: {
131: PetscFunctionBegin;
132: PetscAssertPointer(A, 8);
133: PetscCall(MatCreate(comm, A));
134: if (prefix) PetscCall(MatSetOptionsPrefix(*A, prefix));
135: PetscCall(MatSetBlockSize(*A, bs));
136: PetscCall(MatSetSizes(*A, m, n, M, N));
137: PetscCall(MatSetFromOptions(*A));
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
141: /*@
142: MatSetErrorIfFailure - Causes `Mat` to generate an immediate error, for example a zero pivot, is detected.
144: Logically Collective
146: Input Parameters:
147: + mat - matrix obtained from `MatCreate()`
148: - flg - `PETSC_TRUE` indicates you want the error generated
150: Level: advanced
152: Note:
153: 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
154: or result in a `KSPConvergedReason` indicating the method did not converge.
156: .seealso: [](ch_matrices), `Mat`, `PCSetErrorIfFailure()`, `KSPConvergedReason`, `SNESConvergedReason`
157: @*/
158: PetscErrorCode MatSetErrorIfFailure(Mat mat, PetscBool flg)
159: {
160: PetscFunctionBegin;
163: mat->erroriffailure = flg;
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
167: /*@
168: MatSetSizes - Sets the local and global sizes, and checks to determine compatibility
170: Collective
172: Input Parameters:
173: + A - the matrix
174: . m - number of local rows (or `PETSC_DECIDE`)
175: . n - number of local columns (or `PETSC_DECIDE`)
176: . M - number of global rows (or `PETSC_DETERMINE`)
177: - N - number of global columns (or `PETSC_DETERMINE`)
179: Level: beginner
181: Notes:
182: `m` (`n`) and `M` (`N`) cannot be both `PETSC_DECIDE`
183: If one processor calls this with `M` (`N`) of `PETSC_DECIDE` then all processors must, otherwise the program will hang.
185: If `PETSC_DECIDE` is not used for the arguments 'm' and 'n', then the
186: user must ensure that they are chosen to be compatible with the
187: vectors. To do this, one first considers the matrix-vector product
188: 'y = A x'. The `m` that is used in the above routine must match the
189: local size used in the vector creation routine `VecCreateMPI()` for 'y'.
190: Likewise, the `n` used must match that used as the local size in
191: `VecCreateMPI()` for 'x'.
193: You cannot change the sizes once they have been set.
195: The sizes must be set before `MatSetUp()` or MatXXXSetPreallocation() is called.
197: .seealso: [](ch_matrices), `Mat`, `MatGetSize()`, `PetscSplitOwnership()`
198: @*/
199: PetscErrorCode MatSetSizes(Mat A, PetscInt m, PetscInt n, PetscInt M, PetscInt N)
200: {
201: PetscFunctionBegin;
205: 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);
206: 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);
207: 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,
208: A->rmap->n, A->rmap->N);
209: 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,
210: A->cmap->n, A->cmap->N);
211: A->rmap->n = m;
212: A->cmap->n = n;
213: A->rmap->N = M > -1 ? M : A->rmap->N;
214: A->cmap->N = N > -1 ? N : A->cmap->N;
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: /*@
219: MatSetFromOptions - Creates a matrix where the type is determined
220: from the options database.
222: Collective
224: Input Parameter:
225: . B - the matrix
227: Options Database Keys:
228: + -mat_type seqaij - `MATSEQAIJ` type, uses `MatCreateSeqAIJ()`
229: . -mat_type mpiaij - `MATMPIAIJ` type, uses `MatCreateAIJ()`
230: . -mat_type seqdense - `MATSEQDENSE` type, uses `MatCreateSeqDense()`
231: . -mat_type mpidense - `MATMPIDENSE`, uses `MatCreateDense()`
232: . -mat_type seqbaij - `MATSEQBAIJ`, uses `MatCreateSeqBAIJ()`
233: - -mat_type mpibaij - `MATMPIBAIJ`, uses `MatCreateBAIJ()`
235: See the manpages for particular formats (e.g., `MATSEQAIJ`)
236: for additional format-specific options.
238: Level: beginner
240: Notes:
241: Generates a parallel MPI matrix if the communicator has more than one processor. The default
242: matrix type is `MATAIJ`, using the routines `MatCreateSeqAIJ()` and `MatCreateAIJ()` if you
243: do not select a type in the options database.
245: .seealso: [](ch_matrices), `Mat`, `MatCreateSeqAIJ()`, `MatCreateAIJ()`,
246: `MatCreateSeqDense()`, `MatCreateDense()`,
247: `MatCreateSeqBAIJ()`, `MatCreateBAIJ()`,
248: `MatCreateSeqSBAIJ()`, `MatCreateSBAIJ()`,
249: `MatConvert()`
250: @*/
251: PetscErrorCode MatSetFromOptions(Mat B)
252: {
253: const char *deft = MATAIJ;
254: char type[256];
255: PetscBool flg, set;
256: PetscInt bind_below = 0;
258: PetscFunctionBegin;
261: PetscObjectOptionsBegin((PetscObject)B);
263: if (B->rmap->bs < 0) {
264: PetscInt newbs = -1;
265: PetscCall(PetscOptionsInt("-mat_block_size", "Set the blocksize used to store the matrix", "MatSetBlockSize", newbs, &newbs, &flg));
266: if (flg) {
267: PetscCall(PetscLayoutSetBlockSize(B->rmap, newbs));
268: PetscCall(PetscLayoutSetBlockSize(B->cmap, newbs));
269: }
270: }
272: PetscCall(PetscOptionsFList("-mat_type", "Matrix type", "MatSetType", MatList, deft, type, 256, &flg));
273: if (flg) {
274: PetscCall(MatSetType(B, type));
275: } else if (!((PetscObject)B)->type_name) {
276: PetscCall(MatSetType(B, deft));
277: }
279: PetscCall(PetscOptionsName("-mat_is_symmetric", "Checks if mat is symmetric on MatAssemblyEnd()", "MatIsSymmetric", &B->checksymmetryonassembly));
280: PetscCall(PetscOptionsReal("-mat_is_symmetric", "Checks if mat is symmetric on MatAssemblyEnd()", "MatIsSymmetric", B->checksymmetrytol, &B->checksymmetrytol, NULL));
281: PetscCall(PetscOptionsBool("-mat_null_space_test", "Checks if provided null space is correct in MatAssemblyEnd()", "MatSetNullSpaceTest", B->checknullspaceonassembly, &B->checknullspaceonassembly, NULL));
282: PetscCall(PetscOptionsBool("-mat_error_if_failure", "Generate an error if an error occurs when factoring the matrix", "MatSetErrorIfFailure", B->erroriffailure, &B->erroriffailure, NULL));
284: PetscTryTypeMethod(B, setfromoptions, PetscOptionsObject);
286: flg = PETSC_FALSE;
287: 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));
288: if (set) PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, flg));
289: flg = PETSC_FALSE;
290: 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));
291: if (set) PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, flg));
292: flg = PETSC_FALSE;
293: 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));
294: if (set) PetscCall(MatSetOption(B, MAT_IGNORE_ZERO_ENTRIES, flg));
296: flg = PETSC_FALSE;
297: PetscCall(PetscOptionsBool("-mat_form_explicit_transpose", "Hint to form an explicit transpose for operations like MatMultTranspose", "MatSetOption", flg, &flg, &set));
298: if (set) PetscCall(MatSetOption(B, MAT_FORM_EXPLICIT_TRANSPOSE, flg));
300: /* Bind to CPU if below a user-specified size threshold.
301: * This perhaps belongs in the options for the GPU Mat types, but MatBindToCPU() does nothing when called on non-GPU types,
302: * and putting it here makes is more maintainable than duplicating this for all. */
303: 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));
304: if (flg && B->rmap->n < bind_below) PetscCall(MatBindToCPU(B, PETSC_TRUE));
306: /* process any options handlers added with PetscObjectAddOptionsHandler() */
307: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)B, PetscOptionsObject));
308: PetscOptionsEnd();
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: /*@C
313: MatXAIJSetPreallocation - set preallocation for serial and parallel `MATAIJ`, `MATBAIJ`, and `MATSBAIJ` matrices and their unassembled versions.
315: Collective
317: Input Parameters:
318: + A - matrix being preallocated
319: . bs - block size
320: . dnnz - number of nonzero column blocks per block row of diagonal part of parallel matrix
321: . onnz - number of nonzero column blocks per block row of off-diagonal part of parallel matrix
322: . dnnzu - number of nonzero column blocks per block row of upper-triangular part of diagonal part of parallel matrix
323: - onnzu - number of nonzero column blocks per block row of upper-triangular part of off-diagonal part of parallel matrix
325: Level: beginner
327: .seealso: [](ch_matrices), `Mat`, `MatSeqAIJSetPreallocation()`, `MatMPIAIJSetPreallocation()`, `MatSeqBAIJSetPreallocation()`, `MatMPIBAIJSetPreallocation()`,
328: `MatSeqSBAIJSetPreallocation()`, `MatMPISBAIJSetPreallocation()`,
329: `PetscSplitOwnership()`
330: @*/
331: PetscErrorCode MatXAIJSetPreallocation(Mat A, PetscInt bs, const PetscInt dnnz[], const PetscInt onnz[], const PetscInt dnnzu[], const PetscInt onnzu[])
332: {
333: PetscInt cbs;
334: void (*aij)(void);
335: void (*is)(void);
336: void (*hyp)(void) = NULL;
338: PetscFunctionBegin;
339: if (bs != PETSC_DECIDE) { /* don't mess with an already set block size */
340: PetscCall(MatSetBlockSize(A, bs));
341: }
342: PetscCall(PetscLayoutSetUp(A->rmap));
343: PetscCall(PetscLayoutSetUp(A->cmap));
344: PetscCall(MatGetBlockSizes(A, &bs, &cbs));
345: /* these routines assumes bs == cbs, this should be checked somehow */
346: PetscCall(MatSeqBAIJSetPreallocation(A, bs, 0, dnnz));
347: PetscCall(MatMPIBAIJSetPreallocation(A, bs, 0, dnnz, 0, onnz));
348: PetscCall(MatSeqSBAIJSetPreallocation(A, bs, 0, dnnzu));
349: PetscCall(MatMPISBAIJSetPreallocation(A, bs, 0, dnnzu, 0, onnzu));
350: /*
351: 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
352: good before going on with it.
353: */
354: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatMPIAIJSetPreallocation_C", &aij));
355: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatISSetPreallocation_C", &is));
356: #if defined(PETSC_HAVE_HYPRE)
357: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatHYPRESetPreallocation_C", &hyp));
358: #endif
359: if (!aij && !is && !hyp) PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatSeqAIJSetPreallocation_C", &aij));
360: if (aij || is || hyp) {
361: if (bs == cbs && bs == 1) {
362: PetscCall(MatSeqAIJSetPreallocation(A, 0, dnnz));
363: PetscCall(MatMPIAIJSetPreallocation(A, 0, dnnz, 0, onnz));
364: PetscCall(MatISSetPreallocation(A, 0, dnnz, 0, onnz));
365: #if defined(PETSC_HAVE_HYPRE)
366: PetscCall(MatHYPRESetPreallocation(A, 0, dnnz, 0, onnz));
367: #endif
368: } else { /* Convert block-row precallocation to scalar-row */
369: PetscInt i, m, *sdnnz, *sonnz;
370: PetscCall(MatGetLocalSize(A, &m, NULL));
371: PetscCall(PetscMalloc2((!!dnnz) * m, &sdnnz, (!!onnz) * m, &sonnz));
372: for (i = 0; i < m; i++) {
373: if (dnnz) sdnnz[i] = dnnz[i / bs] * cbs;
374: if (onnz) sonnz[i] = onnz[i / bs] * cbs;
375: }
376: PetscCall(MatSeqAIJSetPreallocation(A, 0, dnnz ? sdnnz : NULL));
377: PetscCall(MatMPIAIJSetPreallocation(A, 0, dnnz ? sdnnz : NULL, 0, onnz ? sonnz : NULL));
378: PetscCall(MatISSetPreallocation(A, 0, dnnz ? sdnnz : NULL, 0, onnz ? sonnz : NULL));
379: #if defined(PETSC_HAVE_HYPRE)
380: PetscCall(MatHYPRESetPreallocation(A, 0, dnnz ? sdnnz : NULL, 0, onnz ? sonnz : NULL));
381: #endif
382: PetscCall(PetscFree2(sdnnz, sonnz));
383: }
384: }
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }
388: /*@C
389: MatHeaderMerge - Merges some information from the header of `C` to `A`; the `C` object is then destroyed
391: Collective, No Fortran Support
393: Input Parameters:
394: + A - a `Mat` being merged into
395: - C - the `Mat` providing the merge information
397: Level: developer
399: Developer Note:
400: This is somewhat different from `MatHeaderReplace()`, it would be nice to merge the code
402: .seealso: `Mat`, `MatHeaderReplace()`
403: @*/
404: PetscErrorCode MatHeaderMerge(Mat A, Mat *C)
405: {
406: PetscInt refct;
407: PetscOps Abops;
408: struct _MatOps Aops;
409: char *mtype, *mname, *mprefix;
410: Mat_Product *product;
411: Mat_Redundant *redundant;
412: PetscObjectState state;
414: PetscFunctionBegin;
417: if (A == *C) PetscFunctionReturn(PETSC_SUCCESS);
418: PetscCheckSameComm(A, 1, *C, 2);
419: /* save the parts of A we need */
420: Abops = ((PetscObject)A)->bops[0];
421: Aops = A->ops[0];
422: refct = ((PetscObject)A)->refct;
423: mtype = ((PetscObject)A)->type_name;
424: mname = ((PetscObject)A)->name;
425: state = ((PetscObject)A)->state;
426: mprefix = ((PetscObject)A)->prefix;
427: product = A->product;
428: redundant = A->redundant;
430: /* zero these so the destroy below does not free them */
431: ((PetscObject)A)->type_name = NULL;
432: ((PetscObject)A)->name = NULL;
434: /*
435: free all the interior data structures from mat
436: cannot use PetscUseTypeMethod(A,destroy); because compiler
437: thinks it may print NULL type_name and name
438: */
439: PetscTryTypeMethod(A, destroy);
441: PetscCall(PetscFree(A->defaultvectype));
442: PetscCall(PetscFree(A->defaultrandtype));
443: PetscCall(PetscLayoutDestroy(&A->rmap));
444: PetscCall(PetscLayoutDestroy(&A->cmap));
445: PetscCall(PetscFunctionListDestroy(&((PetscObject)A)->qlist));
446: PetscCall(PetscObjectListDestroy(&((PetscObject)A)->olist));
447: PetscCall(PetscComposedQuantitiesDestroy((PetscObject)A));
449: /* copy C over to A */
450: PetscCall(PetscFree(A->factorprefix));
451: PetscCall(PetscMemcpy(A, *C, sizeof(struct _p_Mat)));
453: /* return the parts of A we saved */
454: ((PetscObject)A)->bops[0] = Abops;
455: A->ops[0] = Aops;
456: ((PetscObject)A)->refct = refct;
457: ((PetscObject)A)->type_name = mtype;
458: ((PetscObject)A)->name = mname;
459: ((PetscObject)A)->prefix = mprefix;
460: ((PetscObject)A)->state = state + 1;
461: A->product = product;
462: A->redundant = redundant;
464: /* since these two are copied into A we do not want them destroyed in C */
465: ((PetscObject)*C)->qlist = NULL;
466: ((PetscObject)*C)->olist = NULL;
468: PetscCall(PetscHeaderDestroy(C));
469: PetscFunctionReturn(PETSC_SUCCESS);
470: }
472: /*@
473: MatHeaderReplace - Replaces the internal data of matrix `A` by the internal data of matrix `C` while deleting the outer wrapper of `C`
475: Input Parameters:
476: + A - a `Mat` whose internal data is to be replaced
477: - C - the `Mat` providing new internal data for `A`
479: Level: advanced
481: Example Usage\:
482: .vb
483: Mat C;
484: MatCreateSeqAIJWithArrays(..., &C);
485: MatHeaderReplace(A, &C);
486: // C has been destroyed and A contains the matrix entries of C
487: .ve
489: Note:
490: 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
491: (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.
493: Developer Note:
494: This is somewhat different from `MatHeaderMerge()` it would be nice to merge the code
496: .seealso: `Mat`, `MatHeaderMerge()`
497: @*/
498: PetscErrorCode MatHeaderReplace(Mat A, Mat *C)
499: {
500: PetscInt refct;
501: PetscObjectState state;
502: struct _p_Mat buffer;
503: MatStencilInfo stencil;
505: PetscFunctionBegin;
508: if (A == *C) PetscFunctionReturn(PETSC_SUCCESS);
509: PetscCheckSameComm(A, 1, *C, 2);
510: 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);
512: /* swap C and A */
513: refct = ((PetscObject)A)->refct;
514: state = ((PetscObject)A)->state;
515: stencil = A->stencil;
516: PetscCall(PetscMemcpy(&buffer, A, sizeof(struct _p_Mat)));
517: PetscCall(PetscMemcpy(A, *C, sizeof(struct _p_Mat)));
518: PetscCall(PetscMemcpy(*C, &buffer, sizeof(struct _p_Mat)));
519: ((PetscObject)A)->refct = refct;
520: ((PetscObject)A)->state = state + 1;
521: A->stencil = stencil;
523: ((PetscObject)*C)->refct = 1;
524: PetscCall(MatShellSetOperation(*C, MATOP_DESTROY, (void (*)(void))NULL));
525: PetscCall(MatDestroy(C));
526: PetscFunctionReturn(PETSC_SUCCESS);
527: }
529: /*@
530: MatBindToCPU - marks a matrix to temporarily stay on the CPU and perform computations on the CPU
532: Logically Collective
534: Input Parameters:
535: + A - the matrix
536: - flg - bind to the CPU if value of `PETSC_TRUE`
538: Level: intermediate
540: .seealso: [](ch_matrices), `Mat`, `MatBoundToCPU()`
541: @*/
542: PetscErrorCode MatBindToCPU(Mat A, PetscBool flg)
543: {
544: PetscFunctionBegin;
547: #if defined(PETSC_HAVE_DEVICE)
548: if (A->boundtocpu == flg) PetscFunctionReturn(PETSC_SUCCESS);
549: A->boundtocpu = flg;
550: PetscTryTypeMethod(A, bindtocpu, flg);
551: #endif
552: PetscFunctionReturn(PETSC_SUCCESS);
553: }
555: /*@
556: MatBoundToCPU - query if a matrix is bound to the CPU
558: Input Parameter:
559: . A - the matrix
561: Output Parameter:
562: . flg - the logical flag
564: Level: intermediate
566: .seealso: [](ch_matrices), `Mat`, `MatBindToCPU()`
567: @*/
568: PetscErrorCode MatBoundToCPU(Mat A, PetscBool *flg)
569: {
570: PetscFunctionBegin;
572: PetscAssertPointer(flg, 2);
573: #if defined(PETSC_HAVE_DEVICE)
574: *flg = A->boundtocpu;
575: #else
576: *flg = PETSC_TRUE;
577: #endif
578: PetscFunctionReturn(PETSC_SUCCESS);
579: }
581: PetscErrorCode MatSetValuesCOO_Basic(Mat A, const PetscScalar coo_v[], InsertMode imode)
582: {
583: IS is_coo_i, is_coo_j;
584: const PetscInt *coo_i, *coo_j;
585: PetscInt n, n_i, n_j;
586: PetscScalar zero = 0.;
588: PetscFunctionBegin;
589: PetscCall(PetscObjectQuery((PetscObject)A, "__PETSc_coo_i", (PetscObject *)&is_coo_i));
590: PetscCall(PetscObjectQuery((PetscObject)A, "__PETSc_coo_j", (PetscObject *)&is_coo_j));
591: PetscCheck(is_coo_i, PetscObjectComm((PetscObject)A), PETSC_ERR_COR, "Missing coo_i IS");
592: PetscCheck(is_coo_j, PetscObjectComm((PetscObject)A), PETSC_ERR_COR, "Missing coo_j IS");
593: PetscCall(ISGetLocalSize(is_coo_i, &n_i));
594: PetscCall(ISGetLocalSize(is_coo_j, &n_j));
595: PetscCheck(n_i == n_j, PETSC_COMM_SELF, PETSC_ERR_COR, "Wrong local size %" PetscInt_FMT " != %" PetscInt_FMT, n_i, n_j);
596: PetscCall(ISGetIndices(is_coo_i, &coo_i));
597: PetscCall(ISGetIndices(is_coo_j, &coo_j));
598: if (imode != ADD_VALUES) PetscCall(MatZeroEntries(A));
599: for (n = 0; n < n_i; n++) PetscCall(MatSetValue(A, coo_i[n], coo_j[n], coo_v ? coo_v[n] : zero, ADD_VALUES));
600: PetscCall(ISRestoreIndices(is_coo_i, &coo_i));
601: PetscCall(ISRestoreIndices(is_coo_j, &coo_j));
602: PetscFunctionReturn(PETSC_SUCCESS);
603: }
605: PetscErrorCode MatSetPreallocationCOO_Basic(Mat A, PetscCount ncoo, PetscInt coo_i[], PetscInt coo_j[])
606: {
607: Mat preallocator;
608: IS is_coo_i, is_coo_j;
609: PetscScalar zero = 0.0;
611: PetscFunctionBegin;
612: PetscCall(PetscLayoutSetUp(A->rmap));
613: PetscCall(PetscLayoutSetUp(A->cmap));
614: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &preallocator));
615: PetscCall(MatSetType(preallocator, MATPREALLOCATOR));
616: PetscCall(MatSetSizes(preallocator, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
617: PetscCall(MatSetLayouts(preallocator, A->rmap, A->cmap));
618: PetscCall(MatSetUp(preallocator));
619: for (PetscCount n = 0; n < ncoo; n++) PetscCall(MatSetValue(preallocator, coo_i[n], coo_j[n], zero, INSERT_VALUES));
620: PetscCall(MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY));
621: PetscCall(MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY));
622: PetscCall(MatPreallocatorPreallocate(preallocator, PETSC_TRUE, A));
623: PetscCall(MatDestroy(&preallocator));
624: PetscCheck(ncoo <= PETSC_MAX_INT, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ncoo %" PetscCount_FMT " overflowed PetscInt; configure --with-64-bit-indices or request support", ncoo);
625: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncoo, coo_i, PETSC_COPY_VALUES, &is_coo_i));
626: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncoo, coo_j, PETSC_COPY_VALUES, &is_coo_j));
627: PetscCall(PetscObjectCompose((PetscObject)A, "__PETSc_coo_i", (PetscObject)is_coo_i));
628: PetscCall(PetscObjectCompose((PetscObject)A, "__PETSc_coo_j", (PetscObject)is_coo_j));
629: PetscCall(ISDestroy(&is_coo_i));
630: PetscCall(ISDestroy(&is_coo_j));
631: PetscFunctionReturn(PETSC_SUCCESS);
632: }
634: /*@C
635: MatSetPreallocationCOO - set preallocation for matrices using a coordinate format of the entries with global indices
637: Collective
639: Input Parameters:
640: + A - matrix being preallocated
641: . ncoo - number of entries
642: . coo_i - row indices
643: - coo_j - column indices
645: Level: beginner
647: Notes:
648: The indices `coo_i` and `coo_j` may be modified within this function. The caller should not rely on them
649: having any specific value after this function returns. The arrays can be freed or reused immediately
650: after this function returns.
652: Entries can be repeated, see `MatSetValuesCOO()`. Entries with negative row or column indices are allowed
653: but will be ignored. The corresponding entries in `MatSetValuesCOO()` will be ignored too. Remote entries
654: are allowed and will be properly added or inserted to the matrix, unless the matrix option `MAT_IGNORE_OFF_PROC_ENTRIES`
655: 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.
657: If you just want to create a sequential AIJ matrix (`MATSEQAIJ`), and your matrix entries in COO format are unique, you can also use
658: `MatCreateSeqAIJFromTriple()`. But that is not recommended for iterative applications.
660: .seealso: [](ch_matrices), `Mat`, `MatSetValuesCOO()`, `MatSeqAIJSetPreallocation()`, `MatMPIAIJSetPreallocation()`, `MatSeqBAIJSetPreallocation()`,
661: `MatMPIBAIJSetPreallocation()`, `MatSeqSBAIJSetPreallocation()`, `MatMPISBAIJSetPreallocation()`, `MatSetPreallocationCOOLocal()`,
662: `DMSetMatrixPreallocateSkip()`, `MatCreateSeqAIJFromTriple()`
663: @*/
664: PetscErrorCode MatSetPreallocationCOO(Mat A, PetscCount ncoo, PetscInt coo_i[], PetscInt coo_j[])
665: {
666: PetscErrorCode (*f)(Mat, PetscCount, PetscInt[], PetscInt[]) = NULL;
668: PetscFunctionBegin;
671: if (ncoo) PetscAssertPointer(coo_i, 3);
672: if (ncoo) PetscAssertPointer(coo_j, 4);
673: PetscCall(PetscLayoutSetUp(A->rmap));
674: PetscCall(PetscLayoutSetUp(A->cmap));
675: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatSetPreallocationCOO_C", &f));
677: PetscCall(PetscLogEventBegin(MAT_PreallCOO, A, 0, 0, 0));
678: if (f) {
679: PetscCall((*f)(A, ncoo, coo_i, coo_j));
680: } else { /* allow fallback, very slow */
681: PetscCall(MatSetPreallocationCOO_Basic(A, ncoo, coo_i, coo_j));
682: }
683: PetscCall(PetscLogEventEnd(MAT_PreallCOO, A, 0, 0, 0));
684: A->preallocated = PETSC_TRUE;
685: A->nonzerostate++;
686: PetscFunctionReturn(PETSC_SUCCESS);
687: }
689: /*@C
690: MatSetPreallocationCOOLocal - set preallocation for matrices using a coordinate format of the entries with local indices
692: Collective
694: Input Parameters:
695: + A - matrix being preallocated
696: . ncoo - number of entries
697: . coo_i - row indices (local numbering; may be modified)
698: - coo_j - column indices (local numbering; may be modified)
700: Level: beginner
702: Notes:
703: The local indices are translated using the local to global mapping, thus `MatSetLocalToGlobalMapping()` must have been
704: called prior to this function. For matrices created with `DMCreateMatrix()` the local to global mapping is often already provided.
706: The indices `coo_i` and `coo_j` may be modified within this function. They might be translated to corresponding global
707: indices, but the caller should not rely on them having any specific value after this function returns. The arrays
708: can be freed or reused immediately after this function returns.
710: Entries can be repeated, see `MatSetValuesCOO()`. Entries with negative row or column indices are allowed
711: but will be ignored. The corresponding entries in `MatSetValuesCOO()` will be ignored too. Remote entries
712: are allowed and will be properly added or inserted to the matrix.
714: .seealso: [](ch_matrices), `Mat`, `MatSetValuesCOO()`, `MatSeqAIJSetPreallocation()`, `MatMPIAIJSetPreallocation()`, `MatSeqBAIJSetPreallocation()`,
715: `MatMPIBAIJSetPreallocation()`, `MatSeqSBAIJSetPreallocation()`, `MatMPISBAIJSetPreallocation()`, `MatSetPreallocationCOO()`,
716: `DMSetMatrixPreallocateSkip()`
717: @*/
718: PetscErrorCode MatSetPreallocationCOOLocal(Mat A, PetscCount ncoo, PetscInt coo_i[], PetscInt coo_j[])
719: {
720: PetscErrorCode (*f)(Mat, PetscCount, PetscInt[], PetscInt[]) = NULL;
722: PetscFunctionBegin;
725: if (ncoo) PetscAssertPointer(coo_i, 3);
726: if (ncoo) PetscAssertPointer(coo_j, 4);
727: PetscCheck(ncoo <= PETSC_MAX_INT, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ncoo %" PetscCount_FMT " overflowed PetscInt; configure --with-64-bit-indices or request support", ncoo);
728: PetscCall(PetscLayoutSetUp(A->rmap));
729: PetscCall(PetscLayoutSetUp(A->cmap));
731: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatSetPreallocationCOOLocal_C", &f));
732: if (f) {
733: PetscCall((*f)(A, ncoo, coo_i, coo_j));
734: A->nonzerostate++;
735: } else {
736: ISLocalToGlobalMapping ltog_row, ltog_col;
737: PetscCall(MatGetLocalToGlobalMapping(A, <og_row, <og_col));
738: if (ltog_row) PetscCall(ISLocalToGlobalMappingApply(ltog_row, ncoo, coo_i, coo_i));
739: if (ltog_col) PetscCall(ISLocalToGlobalMappingApply(ltog_col, ncoo, coo_j, coo_j));
740: PetscCall(MatSetPreallocationCOO(A, ncoo, coo_i, coo_j));
741: }
742: A->preallocated = PETSC_TRUE;
743: PetscFunctionReturn(PETSC_SUCCESS);
744: }
746: /*@
747: MatSetValuesCOO - set values at once in a matrix preallocated using `MatSetPreallocationCOO()`
749: Collective
751: Input Parameters:
752: + A - matrix being preallocated
753: . coo_v - the matrix values (can be `NULL`)
754: - imode - the insert mode
756: Level: beginner
758: Notes:
759: The values must follow the order of the indices prescribed with `MatSetPreallocationCOO()` or `MatSetPreallocationCOOLocal()`.
761: When repeated entries are specified in the COO indices the `coo_v` values are first properly summed, regardless of the value of imode.
762: The imode flag indicates if coo_v must be added to the current values of the matrix (`ADD_VALUES`) or overwritten (`INSERT_VALUES`).
764: `MatAssemblyBegin()` and `MatAssemblyEnd()` do not need to be called after this routine. It automatically handles the assembly process.
766: .seealso: [](ch_matrices), `Mat`, `MatSetPreallocationCOO()`, `MatSetPreallocationCOOLocal()`, `InsertMode`, `INSERT_VALUES`, `ADD_VALUES`
767: @*/
768: PetscErrorCode MatSetValuesCOO(Mat A, const PetscScalar coo_v[], InsertMode imode)
769: {
770: PetscErrorCode (*f)(Mat, const PetscScalar[], InsertMode) = NULL;
771: PetscBool oldFlg;
773: PetscFunctionBegin;
776: MatCheckPreallocated(A, 1);
778: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatSetValuesCOO_C", &f));
779: PetscCall(PetscLogEventBegin(MAT_SetVCOO, A, 0, 0, 0));
780: if (f) {
781: PetscCall((*f)(A, coo_v, imode)); // all known COO implementations do not use MatStash. They do their own off-proc communication
782: PetscCall(MatGetOption(A, MAT_NO_OFF_PROC_ENTRIES, &oldFlg));
783: PetscCall(MatSetOption(A, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE)); // set A->nooffprocentries to avoid costly MatStash scatter in MatAssembly
784: } else {
785: PetscCall(MatSetValuesCOO_Basic(A, coo_v, imode)); // fall back to MatSetValues, which might use MatStash
786: }
787: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
788: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
789: if (f) PetscCall(MatSetOption(A, MAT_NO_OFF_PROC_ENTRIES, oldFlg));
790: PetscCall(PetscLogEventEnd(MAT_SetVCOO, A, 0, 0, 0));
791: PetscFunctionReturn(PETSC_SUCCESS);
792: }
794: /*@
795: MatSetBindingPropagates - Sets whether the state of being bound to the CPU for a GPU matrix type propagates to child and some other associated objects
797: Input Parameters:
798: + A - the matrix
799: - flg - flag indicating whether the boundtocpu flag should be propagated
801: Level: developer
803: Notes:
804: If the value of flg is set to true, the following will occur
805: + `MatCreateSubMatrices()` and `MatCreateRedundantMatrix()` - bind created matrices to CPU if the input matrix is bound to the CPU.
806: - `MatCreateVecs()` - bind created vectors to CPU if the input matrix is bound to the CPU.
808: The bindingpropagates flag itself is also propagated by the above routines.
810: Developer Notes:
811: If the fine-scale `DMDA` has the `-dm_bind_below` option set to true, then `DMCreateInterpolationScale()` calls `MatSetBindingPropagates()`
812: on the restriction/interpolation operator to set the bindingpropagates flag to true.
814: .seealso: [](ch_matrices), `Mat`, `VecSetBindingPropagates()`, `MatGetBindingPropagates()`
815: @*/
816: PetscErrorCode MatSetBindingPropagates(Mat A, PetscBool flg)
817: {
818: PetscFunctionBegin;
820: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
821: A->bindingpropagates = flg;
822: #endif
823: PetscFunctionReturn(PETSC_SUCCESS);
824: }
826: /*@
827: MatGetBindingPropagates - Gets whether the state of being bound to the CPU for a GPU matrix type propagates to child and some other associated objects
829: Input Parameter:
830: . A - the matrix
832: Output Parameter:
833: . flg - flag indicating whether the boundtocpu flag will be propagated
835: Level: developer
837: .seealso: [](ch_matrices), `Mat`, `MatSetBindingPropagates()`
838: @*/
839: PetscErrorCode MatGetBindingPropagates(Mat A, PetscBool *flg)
840: {
841: PetscFunctionBegin;
843: PetscAssertPointer(flg, 2);
844: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
845: *flg = A->bindingpropagates;
846: #else
847: *flg = PETSC_FALSE;
848: #endif
849: PetscFunctionReturn(PETSC_SUCCESS);
850: }