Actual source code: cdiagonal.c
1: #include <petsc/private/matimpl.h>
3: typedef struct {
4: PetscScalar diag;
5: } Mat_ConstantDiagonal;
7: static PetscErrorCode MatAXPY_ConstantDiagonal(Mat Y, PetscScalar a, Mat X, MatStructure str)
8: {
9: Mat_ConstantDiagonal *yctx = (Mat_ConstantDiagonal *)Y->data;
10: Mat_ConstantDiagonal *xctx = (Mat_ConstantDiagonal *)X->data;
12: PetscFunctionBegin;
13: yctx->diag += a * xctx->diag;
14: PetscFunctionReturn(PETSC_SUCCESS);
15: }
17: static PetscErrorCode MatEqual_ConstantDiagonal(Mat Y, Mat X, PetscBool *equal)
18: {
19: Mat_ConstantDiagonal *yctx = (Mat_ConstantDiagonal *)Y->data;
20: Mat_ConstantDiagonal *xctx = (Mat_ConstantDiagonal *)X->data;
22: PetscFunctionBegin;
23: *equal = (yctx->diag == xctx->diag) ? PETSC_TRUE : PETSC_FALSE;
24: PetscFunctionReturn(PETSC_SUCCESS);
25: }
27: static PetscErrorCode MatGetRow_ConstantDiagonal(Mat A, PetscInt row, PetscInt *ncols, PetscInt *cols[], PetscScalar *vals[])
28: {
29: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;
31: PetscFunctionBegin;
32: if (ncols) *ncols = 1;
33: if (cols) {
34: PetscCall(PetscMalloc1(1, cols));
35: (*cols)[0] = row;
36: }
37: if (vals) {
38: PetscCall(PetscMalloc1(1, vals));
39: (*vals)[0] = ctx->diag;
40: }
41: PetscFunctionReturn(PETSC_SUCCESS);
42: }
44: static PetscErrorCode MatRestoreRow_ConstantDiagonal(Mat A, PetscInt row, PetscInt *ncols, PetscInt *cols[], PetscScalar *vals[])
45: {
46: PetscFunctionBegin;
47: if (cols) PetscCall(PetscFree(*cols));
48: if (vals) PetscCall(PetscFree(*vals));
49: PetscFunctionReturn(PETSC_SUCCESS);
50: }
52: static PetscErrorCode MatMultAdd_ConstantDiagonal(Mat mat, Vec v1, Vec v2, Vec v3)
53: {
54: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data;
56: PetscFunctionBegin;
57: if (v2 == v3) {
58: PetscCall(VecAXPBY(v3, ctx->diag, 1.0, v1));
59: } else {
60: PetscCall(VecAXPBYPCZ(v3, ctx->diag, 1.0, 0.0, v1, v2));
61: }
62: PetscFunctionReturn(PETSC_SUCCESS);
63: }
65: static PetscErrorCode MatMultHermitianTransposeAdd_ConstantDiagonal(Mat mat, Vec v1, Vec v2, Vec v3)
66: {
67: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data;
69: PetscFunctionBegin;
70: if (v2 == v3) {
71: PetscCall(VecAXPBY(v3, PetscConj(ctx->diag), 1.0, v1));
72: } else {
73: PetscCall(VecAXPBYPCZ(v3, PetscConj(ctx->diag), 1.0, 0.0, v1, v2));
74: }
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: static PetscErrorCode MatNorm_ConstantDiagonal(Mat A, NormType type, PetscReal *nrm)
79: {
80: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;
82: PetscFunctionBegin;
83: if (type == NORM_FROBENIUS || type == NORM_2 || type == NORM_1 || type == NORM_INFINITY) *nrm = PetscAbsScalar(ctx->diag);
84: else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Unsupported norm");
85: PetscFunctionReturn(PETSC_SUCCESS);
86: }
88: static PetscErrorCode MatCreateSubMatrices_ConstantDiagonal(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[])
90: {
91: Mat B;
93: PetscFunctionBegin;
94: PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &B));
95: PetscCall(MatCreateSubMatrices(B, n, irow, icol, scall, submat));
96: PetscCall(MatDestroy(&B));
97: PetscFunctionReturn(PETSC_SUCCESS);
98: }
100: static PetscErrorCode MatDuplicate_ConstantDiagonal(Mat A, MatDuplicateOption op, Mat *B)
101: {
102: Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data;
104: PetscFunctionBegin;
105: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B));
106: PetscCall(MatSetSizes(*B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
107: PetscCall(MatSetBlockSizesFromMats(*B, A, A));
108: PetscCall(MatSetType(*B, MATCONSTANTDIAGONAL));
109: PetscCall(PetscLayoutReference(A->rmap, &(*B)->rmap));
110: PetscCall(PetscLayoutReference(A->cmap, &(*B)->cmap));
111: if (op == MAT_COPY_VALUES) {
112: Mat_ConstantDiagonal *bctx = (Mat_ConstantDiagonal *)(*B)->data;
113: bctx->diag = actx->diag;
114: }
115: PetscFunctionReturn(PETSC_SUCCESS);
116: }
118: static PetscErrorCode MatMissingDiagonal_ConstantDiagonal(Mat mat, PetscBool *missing, PetscInt *dd)
119: {
120: PetscFunctionBegin;
121: *missing = PETSC_FALSE;
122: PetscFunctionReturn(PETSC_SUCCESS);
123: }
125: static PetscErrorCode MatDestroy_ConstantDiagonal(Mat mat)
126: {
127: PetscFunctionBegin;
128: PetscCall(PetscFree(mat->data));
129: mat->structural_symmetry_eternal = PETSC_FALSE;
130: mat->symmetry_eternal = PETSC_FALSE;
131: PetscFunctionReturn(PETSC_SUCCESS);
132: }
134: static PetscErrorCode MatView_ConstantDiagonal(Mat J, PetscViewer viewer)
135: {
136: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
137: PetscBool iascii;
139: PetscFunctionBegin;
140: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
141: if (iascii) {
142: PetscViewerFormat format;
144: PetscCall(PetscViewerGetFormat(viewer, &format));
145: if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS);
146: #if defined(PETSC_USE_COMPLEX)
147: PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g + i %g\n", (double)PetscRealPart(ctx->diag), (double)PetscImaginaryPart(ctx->diag)));
148: #else
149: PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g\n", (double)ctx->diag));
150: #endif
151: }
152: PetscFunctionReturn(PETSC_SUCCESS);
153: }
155: static PetscErrorCode MatAssemblyEnd_ConstantDiagonal(Mat J, MatAssemblyType mt)
156: {
157: PetscFunctionBegin;
158: PetscFunctionReturn(PETSC_SUCCESS);
159: }
161: static PetscErrorCode MatMult_ConstantDiagonal(Mat J, Vec x, Vec y)
162: {
163: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
165: PetscFunctionBegin;
166: PetscCall(VecAXPBY(y, ctx->diag, 0.0, x));
167: PetscFunctionReturn(PETSC_SUCCESS);
168: }
170: static PetscErrorCode MatMultHermitianTranspose_ConstantDiagonal(Mat J, Vec x, Vec y)
171: {
172: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
174: PetscFunctionBegin;
175: PetscCall(VecAXPBY(y, PetscConj(ctx->diag), 0.0, x));
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: static PetscErrorCode MatGetDiagonal_ConstantDiagonal(Mat J, Vec x)
180: {
181: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
183: PetscFunctionBegin;
184: PetscCall(VecSet(x, ctx->diag));
185: PetscFunctionReturn(PETSC_SUCCESS);
186: }
188: static PetscErrorCode MatShift_ConstantDiagonal(Mat Y, PetscScalar a)
189: {
190: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
192: PetscFunctionBegin;
193: ctx->diag += a;
194: PetscFunctionReturn(PETSC_SUCCESS);
195: }
197: static PetscErrorCode MatScale_ConstantDiagonal(Mat Y, PetscScalar a)
198: {
199: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
201: PetscFunctionBegin;
202: ctx->diag *= a;
203: PetscFunctionReturn(PETSC_SUCCESS);
204: }
206: static PetscErrorCode MatZeroEntries_ConstantDiagonal(Mat Y)
207: {
208: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
210: PetscFunctionBegin;
211: ctx->diag = 0.0;
212: PetscFunctionReturn(PETSC_SUCCESS);
213: }
215: static PetscErrorCode MatSolve_ConstantDiagonal(Mat matin, Vec b, Vec x)
216: {
217: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)matin->data;
219: PetscFunctionBegin;
220: if (ctx->diag == 0.0) matin->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
221: else matin->factorerrortype = MAT_FACTOR_NOERROR;
222: PetscCall(VecAXPBY(x, 1.0 / ctx->diag, 0.0, b));
223: PetscFunctionReturn(PETSC_SUCCESS);
224: }
226: static PetscErrorCode MatSOR_ConstantDiagonal(Mat matin, Vec x, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec y)
227: {
228: PetscFunctionBegin;
229: PetscCall(MatSolve_ConstantDiagonal(matin, x, y));
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }
233: static PetscErrorCode MatGetInfo_ConstantDiagonal(Mat A, MatInfoType flag, MatInfo *info)
234: {
235: PetscFunctionBegin;
236: info->block_size = 1.0;
237: info->nz_allocated = 1.0;
238: info->nz_used = 1.0;
239: info->nz_unneeded = 0.0;
240: info->assemblies = A->num_ass;
241: info->mallocs = 0.0;
242: info->memory = 0; /* REVIEW ME */
243: if (A->factortype) {
244: info->fill_ratio_given = 1.0;
245: info->fill_ratio_needed = 1.0;
246: info->factor_mallocs = 0.0;
247: } else {
248: info->fill_ratio_given = 0;
249: info->fill_ratio_needed = 0;
250: info->factor_mallocs = 0;
251: }
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: /*@
256: MatCreateConstantDiagonal - Creates a matrix with a uniform value along the diagonal
258: Collective
260: Input Parameters:
261: + comm - MPI communicator
262: . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
263: This value should be the same as the local size used in creating the
264: y vector for the matrix-vector product y = Ax.
265: . n - This value should be the same as the local size used in creating the
266: x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have
267: calculated if `N` is given) For square matrices n is almost always `m`.
268: . M - number of global rows (or `PETSC_DETERMINE` to have calculated if m is given)
269: . N - number of global columns (or `PETSC_DETERMINE` to have calculated if n is given)
270: - diag - the diagonal value
272: Output Parameter:
273: . J - the diagonal matrix
275: Level: advanced
277: Notes:
278: Only supports square matrices with the same number of local rows and columns
280: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`, `MatScale()`, `MatShift()`, `MatMult()`, `MatGetDiagonal()`, `MatGetFactor()`, `MatSolve()`
281: @*/
282: PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscScalar diag, Mat *J)
283: {
284: PetscFunctionBegin;
285: PetscCall(MatCreate(comm, J));
286: PetscCall(MatSetSizes(*J, m, n, M, N));
287: PetscCall(MatSetType(*J, MATCONSTANTDIAGONAL));
288: PetscCall(MatShift(*J, diag));
289: PetscCall(MatSetUp(*J));
290: PetscFunctionReturn(PETSC_SUCCESS);
291: }
293: PETSC_EXTERN PetscErrorCode MatCreate_ConstantDiagonal(Mat A)
294: {
295: Mat_ConstantDiagonal *ctx;
297: PetscFunctionBegin;
298: PetscCall(PetscNew(&ctx));
299: ctx->diag = 0.0;
300: A->data = (void *)ctx;
302: A->assembled = PETSC_TRUE;
303: A->preallocated = PETSC_TRUE;
304: A->structurally_symmetric = PETSC_BOOL3_TRUE;
305: A->structural_symmetry_eternal = PETSC_TRUE;
306: A->symmetric = PETSC_BOOL3_TRUE;
307: if (!PetscDefined(USE_COMPLEX)) A->hermitian = PETSC_BOOL3_TRUE;
308: A->symmetry_eternal = PETSC_TRUE;
310: A->ops->mult = MatMult_ConstantDiagonal;
311: A->ops->multadd = MatMultAdd_ConstantDiagonal;
312: A->ops->multtranspose = MatMult_ConstantDiagonal;
313: A->ops->multtransposeadd = MatMultAdd_ConstantDiagonal;
314: A->ops->multhermitiantranspose = MatMultHermitianTranspose_ConstantDiagonal;
315: A->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_ConstantDiagonal;
316: A->ops->solve = MatSolve_ConstantDiagonal;
317: A->ops->solvetranspose = MatSolve_ConstantDiagonal;
318: A->ops->norm = MatNorm_ConstantDiagonal;
319: A->ops->createsubmatrices = MatCreateSubMatrices_ConstantDiagonal;
320: A->ops->duplicate = MatDuplicate_ConstantDiagonal;
321: A->ops->missingdiagonal = MatMissingDiagonal_ConstantDiagonal;
322: A->ops->getrow = MatGetRow_ConstantDiagonal;
323: A->ops->restorerow = MatRestoreRow_ConstantDiagonal;
324: A->ops->sor = MatSOR_ConstantDiagonal;
325: A->ops->shift = MatShift_ConstantDiagonal;
326: A->ops->scale = MatScale_ConstantDiagonal;
327: A->ops->getdiagonal = MatGetDiagonal_ConstantDiagonal;
328: A->ops->view = MatView_ConstantDiagonal;
329: A->ops->zeroentries = MatZeroEntries_ConstantDiagonal;
330: A->ops->assemblyend = MatAssemblyEnd_ConstantDiagonal;
331: A->ops->destroy = MatDestroy_ConstantDiagonal;
332: A->ops->getinfo = MatGetInfo_ConstantDiagonal;
333: A->ops->equal = MatEqual_ConstantDiagonal;
334: A->ops->axpy = MatAXPY_ConstantDiagonal;
336: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATCONSTANTDIAGONAL));
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
340: static PetscErrorCode MatFactorNumeric_ConstantDiagonal(Mat fact, Mat A, const MatFactorInfo *info)
341: {
342: Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data, *fctx = (Mat_ConstantDiagonal *)fact->data;
344: PetscFunctionBegin;
345: if (actx->diag == 0.0) fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
346: else fact->factorerrortype = MAT_FACTOR_NOERROR;
347: fctx->diag = 1.0 / actx->diag;
348: fact->ops->solve = MatMult_ConstantDiagonal;
349: PetscFunctionReturn(PETSC_SUCCESS);
350: }
352: static PetscErrorCode MatFactorSymbolic_LU_ConstantDiagonal(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
353: {
354: PetscFunctionBegin;
355: fact->ops->lufactornumeric = MatFactorNumeric_ConstantDiagonal;
356: PetscFunctionReturn(PETSC_SUCCESS);
357: }
359: static PetscErrorCode MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact, Mat A, IS isrow, const MatFactorInfo *info)
360: {
361: PetscFunctionBegin;
362: fact->ops->choleskyfactornumeric = MatFactorNumeric_ConstantDiagonal;
363: PetscFunctionReturn(PETSC_SUCCESS);
364: }
366: PETSC_INTERN PetscErrorCode MatGetFactor_constantdiagonal_petsc(Mat A, MatFactorType ftype, Mat *B)
367: {
368: PetscInt n = A->rmap->n, N = A->rmap->N;
370: PetscFunctionBegin;
371: PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), n, n, N, N, 0, B));
373: (*B)->factortype = ftype;
374: (*B)->ops->ilufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal;
375: (*B)->ops->lufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal;
376: (*B)->ops->iccfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;
377: (*B)->ops->choleskyfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;
379: (*B)->ops->shift = NULL;
380: (*B)->ops->scale = NULL;
381: (*B)->ops->mult = NULL;
382: (*B)->ops->sor = NULL;
383: (*B)->ops->zeroentries = NULL;
385: PetscCall(PetscFree((*B)->solvertype));
386: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype));
387: PetscFunctionReturn(PETSC_SUCCESS);
388: }