Actual source code: cdiagonal.c
2: #include <petsc/private/matimpl.h>
4: typedef struct {
5: PetscScalar diag;
6: } Mat_ConstantDiagonal;
8: static PetscErrorCode MatAXPY_ConstantDiagonal(Mat Y, PetscScalar a, Mat X, MatStructure str)
9: {
10: Mat_ConstantDiagonal *yctx = (Mat_ConstantDiagonal *)Y->data;
11: Mat_ConstantDiagonal *xctx = (Mat_ConstantDiagonal *)X->data;
13: PetscFunctionBegin;
14: yctx->diag += a * xctx->diag;
15: PetscFunctionReturn(PETSC_SUCCESS);
16: }
18: static PetscErrorCode MatGetRow_ConstantDiagonal(Mat A, PetscInt row, PetscInt *ncols, PetscInt *cols[], PetscScalar *vals[])
19: {
20: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;
22: PetscFunctionBegin;
23: if (ncols) *ncols = 1;
24: if (cols) {
25: PetscCall(PetscMalloc1(1, cols));
26: (*cols)[0] = row;
27: }
28: if (vals) {
29: PetscCall(PetscMalloc1(1, vals));
30: (*vals)[0] = ctx->diag;
31: }
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
35: static PetscErrorCode MatRestoreRow_ConstantDiagonal(Mat A, PetscInt row, PetscInt *ncols, PetscInt *cols[], PetscScalar *vals[])
36: {
37: PetscFunctionBegin;
38: if (ncols) *ncols = 0;
39: if (cols) PetscCall(PetscFree(*cols));
40: if (vals) PetscCall(PetscFree(*vals));
41: PetscFunctionReturn(PETSC_SUCCESS);
42: }
44: static PetscErrorCode MatMultTranspose_ConstantDiagonal(Mat A, Vec x, Vec y)
45: {
46: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;
48: PetscFunctionBegin;
49: PetscCall(VecAXPBY(y, ctx->diag, 0.0, x));
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: static PetscErrorCode MatMultAdd_ConstantDiagonal(Mat mat, Vec v1, Vec v2, Vec v3)
54: {
55: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data;
57: PetscFunctionBegin;
58: if (v2 == v3) {
59: PetscCall(VecAXPBY(v3, ctx->diag, 1.0, v1));
60: } else {
61: PetscCall(VecAXPBYPCZ(v3, ctx->diag, 1.0, 0.0, v1, v2));
62: }
63: PetscFunctionReturn(PETSC_SUCCESS);
64: }
66: static PetscErrorCode MatMultTransposeAdd_ConstantDiagonal(Mat mat, Vec v1, Vec v2, Vec v3)
67: {
68: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data;
70: PetscFunctionBegin;
71: if (v2 == v3) {
72: PetscCall(VecAXPBY(v3, ctx->diag, 1.0, v1));
73: } else {
74: PetscCall(VecAXPBYPCZ(v3, ctx->diag, 1.0, 0.0, v1, v2));
75: }
76: PetscFunctionReturn(PETSC_SUCCESS);
77: }
79: static PetscErrorCode MatNorm_ConstantDiagonal(Mat A, NormType type, PetscReal *nrm)
80: {
81: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;
83: PetscFunctionBegin;
84: if (type == NORM_FROBENIUS || type == NORM_2 || type == NORM_1 || type == NORM_INFINITY) *nrm = PetscAbsScalar(ctx->diag);
85: else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Unsupported norm");
86: PetscFunctionReturn(PETSC_SUCCESS);
87: }
89: static PetscErrorCode MatCreateSubMatrices_ConstantDiagonal(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[])
91: {
92: Mat B;
94: PetscFunctionBegin;
95: PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &B));
96: PetscCall(MatCreateSubMatrices(B, n, irow, icol, scall, submat));
97: PetscCall(MatDestroy(&B));
98: PetscFunctionReturn(PETSC_SUCCESS);
99: }
101: static PetscErrorCode MatDuplicate_ConstantDiagonal(Mat A, MatDuplicateOption op, Mat *B)
102: {
103: Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data;
105: PetscFunctionBegin;
106: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B));
107: PetscCall(MatSetSizes(*B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
108: PetscCall(MatSetBlockSizesFromMats(*B, A, A));
109: PetscCall(MatSetType(*B, MATCONSTANTDIAGONAL));
110: PetscCall(PetscLayoutReference(A->rmap, &(*B)->rmap));
111: PetscCall(PetscLayoutReference(A->cmap, &(*B)->cmap));
112: if (op == MAT_COPY_VALUES) {
113: Mat_ConstantDiagonal *bctx = (Mat_ConstantDiagonal *)(*B)->data;
114: bctx->diag = actx->diag;
115: }
116: PetscFunctionReturn(PETSC_SUCCESS);
117: }
119: static PetscErrorCode MatMissingDiagonal_ConstantDiagonal(Mat mat, PetscBool *missing, PetscInt *dd)
120: {
121: PetscFunctionBegin;
122: *missing = PETSC_FALSE;
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: static PetscErrorCode MatDestroy_ConstantDiagonal(Mat mat)
127: {
128: PetscFunctionBegin;
129: PetscCall(PetscFree(mat->data));
130: PetscFunctionReturn(PETSC_SUCCESS);
131: }
133: static PetscErrorCode MatView_ConstantDiagonal(Mat J, PetscViewer viewer)
134: {
135: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
136: PetscBool iascii;
138: PetscFunctionBegin;
139: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
140: if (iascii) {
141: PetscViewerFormat format;
143: PetscCall(PetscViewerGetFormat(viewer, &format));
144: if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS);
145: #if defined(PETSC_USE_COMPLEX)
146: PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g + i %g\n", (double)PetscRealPart(ctx->diag), (double)PetscImaginaryPart(ctx->diag)));
147: #else
148: PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g\n", (double)(ctx->diag)));
149: #endif
150: }
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: static PetscErrorCode MatAssemblyEnd_ConstantDiagonal(Mat J, MatAssemblyType mt)
155: {
156: PetscFunctionBegin;
157: PetscFunctionReturn(PETSC_SUCCESS);
158: }
160: static PetscErrorCode MatMult_ConstantDiagonal(Mat J, Vec x, Vec y)
161: {
162: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
164: PetscFunctionBegin;
165: PetscCall(VecAXPBY(y, ctx->diag, 0.0, x));
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
169: PetscErrorCode MatGetDiagonal_ConstantDiagonal(Mat J, Vec x)
170: {
171: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
173: PetscFunctionBegin;
174: PetscCall(VecSet(x, ctx->diag));
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: static PetscErrorCode MatShift_ConstantDiagonal(Mat Y, PetscScalar a)
179: {
180: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
182: PetscFunctionBegin;
183: ctx->diag += a;
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: static PetscErrorCode MatScale_ConstantDiagonal(Mat Y, PetscScalar a)
188: {
189: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
191: PetscFunctionBegin;
192: ctx->diag *= a;
193: PetscFunctionReturn(PETSC_SUCCESS);
194: }
196: static PetscErrorCode MatZeroEntries_ConstantDiagonal(Mat Y)
197: {
198: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
200: PetscFunctionBegin;
201: ctx->diag = 0.0;
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
205: PetscErrorCode MatSOR_ConstantDiagonal(Mat matin, Vec x, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec y)
206: {
207: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)matin->data;
209: PetscFunctionBegin;
210: if (ctx->diag == 0.0) matin->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
211: else matin->factorerrortype = MAT_FACTOR_NOERROR;
212: PetscCall(VecAXPBY(y, 1.0 / ctx->diag, 0.0, x));
213: PetscFunctionReturn(PETSC_SUCCESS);
214: }
216: PetscErrorCode MatGetInfo_ConstantDiagonal(Mat A, MatInfoType flag, MatInfo *info)
217: {
218: PetscFunctionBegin;
219: info->block_size = 1.0;
220: info->nz_allocated = 1.0;
221: info->nz_used = 1.0;
222: info->nz_unneeded = 0.0;
223: info->assemblies = A->num_ass;
224: info->mallocs = 0.0;
225: info->memory = 0; /* REVIEW ME */
226: if (A->factortype) {
227: info->fill_ratio_given = 1.0;
228: info->fill_ratio_needed = 1.0;
229: info->factor_mallocs = 0.0;
230: } else {
231: info->fill_ratio_given = 0;
232: info->fill_ratio_needed = 0;
233: info->factor_mallocs = 0;
234: }
235: PetscFunctionReturn(PETSC_SUCCESS);
236: }
238: /*@
239: MatCreateConstantDiagonal - Creates a matrix with a uniform value along the diagonal
241: Collective
243: Input Parameters:
244: + comm - MPI communicator
245: . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
246: This value should be the same as the local size used in creating the
247: y vector for the matrix-vector product y = Ax.
248: . n - This value should be the same as the local size used in creating the
249: x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have
250: calculated if `N` is given) For square matrices n is almost always `m`.
251: . M - number of global rows (or `PETSC_DETERMINE` to have calculated if m is given)
252: . N - number of global columns (or `PETSC_DETERMINE` to have calculated if n is given)
253: - diag - the diagonal value
255: Output Parameter:
256: . J - the diagonal matrix
258: Level: advanced
260: Notes:
261: Only supports square matrices with the same number of local rows and columns
263: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`, `MatScale()`, `MatShift()`, `MatMult()`, `MatGetDiagonal()`, `MatGetFactor()`, `MatSolve()`
264: @*/
265: PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscScalar diag, Mat *J)
266: {
267: PetscFunctionBegin;
268: PetscCall(MatCreate(comm, J));
269: PetscCall(MatSetSizes(*J, m, n, M, N));
270: PetscCall(MatSetType(*J, MATCONSTANTDIAGONAL));
271: PetscCall(MatShift(*J, diag));
272: PetscCall(MatSetUp(*J));
273: PetscFunctionReturn(PETSC_SUCCESS);
274: }
276: PETSC_EXTERN PetscErrorCode MatCreate_ConstantDiagonal(Mat A)
277: {
278: Mat_ConstantDiagonal *ctx;
280: PetscFunctionBegin;
281: PetscCall(PetscNew(&ctx));
282: ctx->diag = 0.0;
283: A->data = (void *)ctx;
285: A->assembled = PETSC_TRUE;
286: A->preallocated = PETSC_TRUE;
288: A->ops->mult = MatMult_ConstantDiagonal;
289: A->ops->multadd = MatMultAdd_ConstantDiagonal;
290: A->ops->multtranspose = MatMultTranspose_ConstantDiagonal;
291: A->ops->multtransposeadd = MatMultTransposeAdd_ConstantDiagonal;
292: A->ops->norm = MatNorm_ConstantDiagonal;
293: A->ops->createsubmatrices = MatCreateSubMatrices_ConstantDiagonal;
294: A->ops->duplicate = MatDuplicate_ConstantDiagonal;
295: A->ops->missingdiagonal = MatMissingDiagonal_ConstantDiagonal;
296: A->ops->getrow = MatGetRow_ConstantDiagonal;
297: A->ops->restorerow = MatRestoreRow_ConstantDiagonal;
298: A->ops->sor = MatSOR_ConstantDiagonal;
299: A->ops->shift = MatShift_ConstantDiagonal;
300: A->ops->scale = MatScale_ConstantDiagonal;
301: A->ops->getdiagonal = MatGetDiagonal_ConstantDiagonal;
302: A->ops->view = MatView_ConstantDiagonal;
303: A->ops->zeroentries = MatZeroEntries_ConstantDiagonal;
304: A->ops->assemblyend = MatAssemblyEnd_ConstantDiagonal;
305: A->ops->destroy = MatDestroy_ConstantDiagonal;
306: A->ops->getinfo = MatGetInfo_ConstantDiagonal;
307: A->ops->axpy = MatAXPY_ConstantDiagonal;
309: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATCONSTANTDIAGONAL));
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: static PetscErrorCode MatFactorNumeric_ConstantDiagonal(Mat fact, Mat A, const MatFactorInfo *info)
314: {
315: Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data, *fctx = (Mat_ConstantDiagonal *)fact->data;
317: PetscFunctionBegin;
318: if (actx->diag == 0.0) fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
319: else fact->factorerrortype = MAT_FACTOR_NOERROR;
320: fctx->diag = 1.0 / actx->diag;
321: fact->ops->solve = MatMult_ConstantDiagonal;
322: PetscFunctionReturn(PETSC_SUCCESS);
323: }
325: static PetscErrorCode MatFactorSymbolic_LU_ConstantDiagonal(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
326: {
327: PetscFunctionBegin;
328: fact->ops->lufactornumeric = MatFactorNumeric_ConstantDiagonal;
329: PetscFunctionReturn(PETSC_SUCCESS);
330: }
332: static PetscErrorCode MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact, Mat A, IS isrow, const MatFactorInfo *info)
333: {
334: PetscFunctionBegin;
335: fact->ops->choleskyfactornumeric = MatFactorNumeric_ConstantDiagonal;
336: PetscFunctionReturn(PETSC_SUCCESS);
337: }
339: PETSC_INTERN PetscErrorCode MatGetFactor_constantdiagonal_petsc(Mat A, MatFactorType ftype, Mat *B)
340: {
341: PetscInt n = A->rmap->n, N = A->rmap->N;
343: PetscFunctionBegin;
344: PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), n, n, N, N, 0, B));
346: (*B)->factortype = ftype;
347: (*B)->ops->ilufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal;
348: (*B)->ops->lufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal;
349: (*B)->ops->iccfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;
350: (*B)->ops->choleskyfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;
352: (*B)->ops->shift = NULL;
353: (*B)->ops->scale = NULL;
354: (*B)->ops->mult = NULL;
355: (*B)->ops->sor = NULL;
356: (*B)->ops->zeroentries = NULL;
358: PetscCall(PetscFree((*B)->solvertype));
359: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype));
360: PetscFunctionReturn(PETSC_SUCCESS);
361: }