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: yctx->diag += a * xctx->diag;
14: return 0;
15: }
17: static PetscErrorCode MatGetRow_ConstantDiagonal(Mat A, PetscInt row, PetscInt *ncols, PetscInt *cols[], PetscScalar *vals[])
18: {
19: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;
21: if (ncols) *ncols = 1;
22: if (cols) {
23: PetscMalloc1(1, cols);
24: (*cols)[0] = row;
25: }
26: if (vals) {
27: PetscMalloc1(1, vals);
28: (*vals)[0] = ctx->diag;
29: }
30: return 0;
31: }
33: static PetscErrorCode MatRestoreRow_ConstantDiagonal(Mat A, PetscInt row, PetscInt *ncols, PetscInt *cols[], PetscScalar *vals[])
34: {
35: if (ncols) *ncols = 0;
36: if (cols) PetscFree(*cols);
37: if (vals) PetscFree(*vals);
38: return 0;
39: }
41: static PetscErrorCode MatMultTranspose_ConstantDiagonal(Mat A, Vec x, Vec y)
42: {
43: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;
45: VecAXPBY(y, ctx->diag, 0.0, x);
46: return 0;
47: }
49: static PetscErrorCode MatMultAdd_ConstantDiagonal(Mat mat, Vec v1, Vec v2, Vec v3)
50: {
51: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data;
53: if (v2 == v3) {
54: VecAXPBY(v3, ctx->diag, 1.0, v1);
55: } else {
56: VecAXPBYPCZ(v3, ctx->diag, 1.0, 0.0, v1, v2);
57: }
58: return 0;
59: }
61: static PetscErrorCode MatMultTransposeAdd_ConstantDiagonal(Mat mat, Vec v1, Vec v2, Vec v3)
62: {
63: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data;
65: if (v2 == v3) {
66: VecAXPBY(v3, ctx->diag, 1.0, v1);
67: } else {
68: VecAXPBYPCZ(v3, ctx->diag, 1.0, 0.0, v1, v2);
69: }
70: return 0;
71: }
73: static PetscErrorCode MatNorm_ConstantDiagonal(Mat A, NormType type, PetscReal *nrm)
74: {
75: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;
77: if (type == NORM_FROBENIUS || type == NORM_2 || type == NORM_1 || type == NORM_INFINITY) *nrm = PetscAbsScalar(ctx->diag);
78: else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Unsupported norm");
79: return 0;
80: }
82: static PetscErrorCode MatCreateSubMatrices_ConstantDiagonal(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[])
84: {
85: Mat B;
87: MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &B);
88: MatCreateSubMatrices(B, n, irow, icol, scall, submat);
89: MatDestroy(&B);
90: return 0;
91: }
93: static PetscErrorCode MatDuplicate_ConstantDiagonal(Mat A, MatDuplicateOption op, Mat *B)
94: {
95: Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data;
97: MatCreate(PetscObjectComm((PetscObject)A), B);
98: MatSetSizes(*B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N);
99: MatSetBlockSizesFromMats(*B, A, A);
100: MatSetType(*B, MATCONSTANTDIAGONAL);
101: PetscLayoutReference(A->rmap, &(*B)->rmap);
102: PetscLayoutReference(A->cmap, &(*B)->cmap);
103: if (op == MAT_COPY_VALUES) {
104: Mat_ConstantDiagonal *bctx = (Mat_ConstantDiagonal *)(*B)->data;
105: bctx->diag = actx->diag;
106: }
107: return 0;
108: }
110: static PetscErrorCode MatMissingDiagonal_ConstantDiagonal(Mat mat, PetscBool *missing, PetscInt *dd)
111: {
112: *missing = PETSC_FALSE;
113: return 0;
114: }
116: static PetscErrorCode MatDestroy_ConstantDiagonal(Mat mat)
117: {
118: PetscFree(mat->data);
119: return 0;
120: }
122: static PetscErrorCode MatView_ConstantDiagonal(Mat J, PetscViewer viewer)
123: {
124: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
125: PetscBool iascii;
127: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
128: if (iascii) {
129: PetscViewerFormat format;
131: PetscViewerGetFormat(viewer, &format);
132: if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) return 0;
133: #if defined(PETSC_USE_COMPLEX)
134: PetscViewerASCIIPrintf(viewer, "Diagonal value: %g + i %g\n", (double)PetscRealPart(ctx->diag), (double)PetscImaginaryPart(ctx->diag));
135: #else
136: PetscViewerASCIIPrintf(viewer, "Diagonal value: %g\n", (double)(ctx->diag));
137: #endif
138: }
139: return 0;
140: }
142: static PetscErrorCode MatAssemblyEnd_ConstantDiagonal(Mat J, MatAssemblyType mt)
143: {
144: return 0;
145: }
147: static PetscErrorCode MatMult_ConstantDiagonal(Mat J, Vec x, Vec y)
148: {
149: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
151: VecAXPBY(y, ctx->diag, 0.0, x);
152: return 0;
153: }
155: PetscErrorCode MatGetDiagonal_ConstantDiagonal(Mat J, Vec x)
156: {
157: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
159: VecSet(x, ctx->diag);
160: return 0;
161: }
163: static PetscErrorCode MatShift_ConstantDiagonal(Mat Y, PetscScalar a)
164: {
165: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
167: ctx->diag += a;
168: return 0;
169: }
171: static PetscErrorCode MatScale_ConstantDiagonal(Mat Y, PetscScalar a)
172: {
173: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
175: ctx->diag *= a;
176: return 0;
177: }
179: static PetscErrorCode MatZeroEntries_ConstantDiagonal(Mat Y)
180: {
181: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
183: ctx->diag = 0.0;
184: return 0;
185: }
187: PetscErrorCode MatSOR_ConstantDiagonal(Mat matin, Vec x, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec y)
188: {
189: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)matin->data;
191: if (ctx->diag == 0.0) matin->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
192: else matin->factorerrortype = MAT_FACTOR_NOERROR;
193: VecAXPBY(y, 1.0 / ctx->diag, 0.0, x);
194: return 0;
195: }
197: PetscErrorCode MatGetInfo_ConstantDiagonal(Mat A, MatInfoType flag, MatInfo *info)
198: {
199: info->block_size = 1.0;
200: info->nz_allocated = 1.0;
201: info->nz_used = 1.0;
202: info->nz_unneeded = 0.0;
203: info->assemblies = A->num_ass;
204: info->mallocs = 0.0;
205: info->memory = 0; /* REVIEW ME */
206: if (A->factortype) {
207: info->fill_ratio_given = 1.0;
208: info->fill_ratio_needed = 1.0;
209: info->factor_mallocs = 0.0;
210: } else {
211: info->fill_ratio_given = 0;
212: info->fill_ratio_needed = 0;
213: info->factor_mallocs = 0;
214: }
215: return 0;
216: }
218: /*@
219: MatCreateConstantDiagonal - Creates a matrix with a uniform value along the diagonal
221: Collective
223: Input Parameters:
224: + comm - MPI communicator
225: . m - number of local rows (or `PETSC_DECIDE` to have calculated if M is given)
226: This value should be the same as the local size used in creating the
227: y vector for the matrix-vector product y = Ax.
228: . n - This value should be the same as the local size used in creating the
229: x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
230: calculated if N is given) For square matrices n is almost always m.
231: . M - number of global rows (or `PETSC_DETERMINE` to have calculated if m is given)
232: . N - number of global columns (or `PETSC_DETERMINE` to have calculated if n is given)
233: - diag - the diagonal value
235: Output Parameter:
236: . J - the diagonal matrix
238: Level: advanced
240: Notes:
241: Only supports square matrices with the same number of local rows and columns
243: .seealso: `MatDestroy()`, `MATCONSTANTDIAGONAL`, `MatScale()`, `MatShift()`, `MatMult()`, `MatGetDiagonal()`, `MatGetFactor()`, `MatSolve()`
244: @*/
245: PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscScalar diag, Mat *J)
246: {
247: MatCreate(comm, J);
248: MatSetSizes(*J, m, n, M, N);
249: MatSetType(*J, MATCONSTANTDIAGONAL);
250: MatShift(*J, diag);
251: MatSetUp(*J);
252: return 0;
253: }
255: PETSC_EXTERN PetscErrorCode MatCreate_ConstantDiagonal(Mat A)
256: {
257: Mat_ConstantDiagonal *ctx;
259: PetscNew(&ctx);
260: ctx->diag = 0.0;
261: A->data = (void *)ctx;
263: A->assembled = PETSC_TRUE;
264: A->preallocated = PETSC_TRUE;
266: A->ops->mult = MatMult_ConstantDiagonal;
267: A->ops->multadd = MatMultAdd_ConstantDiagonal;
268: A->ops->multtranspose = MatMultTranspose_ConstantDiagonal;
269: A->ops->multtransposeadd = MatMultTransposeAdd_ConstantDiagonal;
270: A->ops->norm = MatNorm_ConstantDiagonal;
271: A->ops->createsubmatrices = MatCreateSubMatrices_ConstantDiagonal;
272: A->ops->duplicate = MatDuplicate_ConstantDiagonal;
273: A->ops->missingdiagonal = MatMissingDiagonal_ConstantDiagonal;
274: A->ops->getrow = MatGetRow_ConstantDiagonal;
275: A->ops->restorerow = MatRestoreRow_ConstantDiagonal;
276: A->ops->sor = MatSOR_ConstantDiagonal;
277: A->ops->shift = MatShift_ConstantDiagonal;
278: A->ops->scale = MatScale_ConstantDiagonal;
279: A->ops->getdiagonal = MatGetDiagonal_ConstantDiagonal;
280: A->ops->view = MatView_ConstantDiagonal;
281: A->ops->zeroentries = MatZeroEntries_ConstantDiagonal;
282: A->ops->assemblyend = MatAssemblyEnd_ConstantDiagonal;
283: A->ops->destroy = MatDestroy_ConstantDiagonal;
284: A->ops->getinfo = MatGetInfo_ConstantDiagonal;
285: A->ops->axpy = MatAXPY_ConstantDiagonal;
287: PetscObjectChangeTypeName((PetscObject)A, MATCONSTANTDIAGONAL);
288: return 0;
289: }
291: static PetscErrorCode MatFactorNumeric_ConstantDiagonal(Mat fact, Mat A, const MatFactorInfo *info)
292: {
293: Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data, *fctx = (Mat_ConstantDiagonal *)fact->data;
295: if (actx->diag == 0.0) fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
296: else fact->factorerrortype = MAT_FACTOR_NOERROR;
297: fctx->diag = 1.0 / actx->diag;
298: fact->ops->solve = MatMult_ConstantDiagonal;
299: return 0;
300: }
302: static PetscErrorCode MatFactorSymbolic_LU_ConstantDiagonal(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
303: {
304: fact->ops->lufactornumeric = MatFactorNumeric_ConstantDiagonal;
305: return 0;
306: }
308: static PetscErrorCode MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact, Mat A, IS isrow, const MatFactorInfo *info)
309: {
310: fact->ops->choleskyfactornumeric = MatFactorNumeric_ConstantDiagonal;
311: return 0;
312: }
314: PETSC_INTERN PetscErrorCode MatGetFactor_constantdiagonal_petsc(Mat A, MatFactorType ftype, Mat *B)
315: {
316: PetscInt n = A->rmap->n, N = A->rmap->N;
318: MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), n, n, N, N, 0, B);
320: (*B)->factortype = ftype;
321: (*B)->ops->ilufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal;
322: (*B)->ops->lufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal;
323: (*B)->ops->iccfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;
324: (*B)->ops->choleskyfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;
326: (*B)->ops->shift = NULL;
327: (*B)->ops->scale = NULL;
328: (*B)->ops->mult = NULL;
329: (*B)->ops->sor = NULL;
330: (*B)->ops->zeroentries = NULL;
332: PetscFree((*B)->solvertype);
333: PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype);
334: return 0;
335: }