1 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
2
3 typedef struct {
4 PetscScalar diag;
5 } Mat_ConstantDiagonal;
6
MatAXPY_ConstantDiagonal(Mat Y,PetscScalar a,Mat X,MatStructure str)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;
11
12 PetscFunctionBegin;
13 yctx->diag += a * xctx->diag;
14 PetscFunctionReturn(PETSC_SUCCESS);
15 }
16
MatEqual_ConstantDiagonal(Mat Y,Mat X,PetscBool * equal)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;
21
22 PetscFunctionBegin;
23 *equal = (yctx->diag == xctx->diag) ? PETSC_TRUE : PETSC_FALSE;
24 PetscFunctionReturn(PETSC_SUCCESS);
25 }
26
MatGetRow_ConstantDiagonal(Mat A,PetscInt row,PetscInt * ncols,PetscInt * cols[],PetscScalar * vals[])27 static PetscErrorCode MatGetRow_ConstantDiagonal(Mat A, PetscInt row, PetscInt *ncols, PetscInt *cols[], PetscScalar *vals[])
28 {
29 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;
30
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 }
43
MatRestoreRow_ConstantDiagonal(Mat A,PetscInt row,PetscInt * ncols,PetscInt * cols[],PetscScalar * vals[])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 }
51
MatMultAdd_ConstantDiagonal(Mat mat,Vec v1,Vec v2,Vec v3)52 static PetscErrorCode MatMultAdd_ConstantDiagonal(Mat mat, Vec v1, Vec v2, Vec v3)
53 {
54 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data;
55
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 }
64
MatMultHermitianTransposeAdd_ConstantDiagonal(Mat mat,Vec v1,Vec v2,Vec v3)65 static PetscErrorCode MatMultHermitianTransposeAdd_ConstantDiagonal(Mat mat, Vec v1, Vec v2, Vec v3)
66 {
67 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data;
68
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 }
77
MatNorm_ConstantDiagonal(Mat A,NormType type,PetscReal * nrm)78 static PetscErrorCode MatNorm_ConstantDiagonal(Mat A, NormType type, PetscReal *nrm)
79 {
80 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;
81
82 PetscFunctionBegin;
83 PetscCheck(type == NORM_FROBENIUS || type == NORM_2 || type == NORM_1 || type == NORM_INFINITY, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Unsupported norm");
84 *nrm = PetscAbsScalar(ctx->diag);
85 PetscFunctionReturn(PETSC_SUCCESS);
86 }
87
MatCreateSubMatrices_ConstantDiagonal(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat * submat[])88 static PetscErrorCode MatCreateSubMatrices_ConstantDiagonal(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[])
89 {
90 Mat B;
91
92 PetscFunctionBegin;
93 PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &B));
94 PetscCall(MatCreateSubMatrices(B, n, irow, icol, scall, submat));
95 PetscCall(MatDestroy(&B));
96 PetscFunctionReturn(PETSC_SUCCESS);
97 }
98
MatDuplicate_ConstantDiagonal(Mat A,MatDuplicateOption op,Mat * B)99 static PetscErrorCode MatDuplicate_ConstantDiagonal(Mat A, MatDuplicateOption op, Mat *B)
100 {
101 Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data;
102
103 PetscFunctionBegin;
104 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B));
105 PetscCall(MatSetSizes(*B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
106 PetscCall(MatSetBlockSizesFromMats(*B, A, A));
107 PetscCall(MatSetType(*B, MATCONSTANTDIAGONAL));
108 PetscCall(PetscLayoutReference(A->rmap, &(*B)->rmap));
109 PetscCall(PetscLayoutReference(A->cmap, &(*B)->cmap));
110 if (op == MAT_COPY_VALUES) {
111 Mat_ConstantDiagonal *bctx = (Mat_ConstantDiagonal *)(*B)->data;
112 bctx->diag = actx->diag;
113 }
114 PetscFunctionReturn(PETSC_SUCCESS);
115 }
116
MatDestroy_ConstantDiagonal(Mat mat)117 static PetscErrorCode MatDestroy_ConstantDiagonal(Mat mat)
118 {
119 PetscFunctionBegin;
120 PetscCall(PetscFree(mat->data));
121 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConstantDiagonalGetConstant_C", NULL));
122 PetscFunctionReturn(PETSC_SUCCESS);
123 }
124
MatView_ConstantDiagonal(Mat J,PetscViewer viewer)125 static PetscErrorCode MatView_ConstantDiagonal(Mat J, PetscViewer viewer)
126 {
127 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
128 PetscBool isascii;
129
130 PetscFunctionBegin;
131 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
132 if (isascii) {
133 PetscViewerFormat format;
134
135 PetscCall(PetscViewerGetFormat(viewer, &format));
136 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS);
137 if (PetscImaginaryPart(ctx->diag) == 0) {
138 PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g\n", (double)PetscRealPart(ctx->diag)));
139 } else {
140 PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g + i %g\n", (double)PetscRealPart(ctx->diag), (double)PetscImaginaryPart(ctx->diag)));
141 }
142 }
143 PetscFunctionReturn(PETSC_SUCCESS);
144 }
145
MatMult_ConstantDiagonal(Mat J,Vec x,Vec y)146 static PetscErrorCode MatMult_ConstantDiagonal(Mat J, Vec x, Vec y)
147 {
148 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
149
150 PetscFunctionBegin;
151 PetscCall(VecAXPBY(y, ctx->diag, 0.0, x));
152 PetscFunctionReturn(PETSC_SUCCESS);
153 }
154
MatMultHermitianTranspose_ConstantDiagonal(Mat J,Vec x,Vec y)155 static PetscErrorCode MatMultHermitianTranspose_ConstantDiagonal(Mat J, Vec x, Vec y)
156 {
157 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
158
159 PetscFunctionBegin;
160 PetscCall(VecAXPBY(y, PetscConj(ctx->diag), 0.0, x));
161 PetscFunctionReturn(PETSC_SUCCESS);
162 }
163
MatGetDiagonal_ConstantDiagonal(Mat J,Vec x)164 static PetscErrorCode MatGetDiagonal_ConstantDiagonal(Mat J, Vec x)
165 {
166 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
167
168 PetscFunctionBegin;
169 PetscCall(VecSet(x, ctx->diag));
170 PetscFunctionReturn(PETSC_SUCCESS);
171 }
172
MatShift_ConstantDiagonal(Mat Y,PetscScalar a)173 static PetscErrorCode MatShift_ConstantDiagonal(Mat Y, PetscScalar a)
174 {
175 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
176
177 PetscFunctionBegin;
178 ctx->diag += a;
179 PetscFunctionReturn(PETSC_SUCCESS);
180 }
181
MatScale_ConstantDiagonal(Mat Y,PetscScalar a)182 static PetscErrorCode MatScale_ConstantDiagonal(Mat Y, PetscScalar a)
183 {
184 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
185
186 PetscFunctionBegin;
187 ctx->diag *= a;
188 PetscFunctionReturn(PETSC_SUCCESS);
189 }
190
MatZeroEntries_ConstantDiagonal(Mat Y)191 static PetscErrorCode MatZeroEntries_ConstantDiagonal(Mat Y)
192 {
193 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
194
195 PetscFunctionBegin;
196 ctx->diag = 0.0;
197 PetscFunctionReturn(PETSC_SUCCESS);
198 }
199
MatConjugate_ConstantDiagonal(Mat Y)200 static PetscErrorCode MatConjugate_ConstantDiagonal(Mat Y)
201 {
202 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
203
204 PetscFunctionBegin;
205 ctx->diag = PetscConj(ctx->diag);
206 PetscFunctionReturn(PETSC_SUCCESS);
207 }
208
MatTranspose_ConstantDiagonal(Mat A,MatReuse reuse,Mat * matout)209 static PetscErrorCode MatTranspose_ConstantDiagonal(Mat A, MatReuse reuse, Mat *matout)
210 {
211 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;
212
213 PetscFunctionBegin;
214 if (reuse == MAT_INPLACE_MATRIX) {
215 PetscLayout tmplayout = A->rmap;
216
217 A->rmap = A->cmap;
218 A->cmap = tmplayout;
219 } else {
220 if (reuse == MAT_INITIAL_MATRIX) {
221 PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), A->cmap->n, A->rmap->n, A->cmap->N, A->rmap->N, ctx->diag, matout));
222 } else {
223 PetscCall(MatZeroEntries(*matout));
224 PetscCall(MatShift(*matout, ctx->diag));
225 }
226 }
227 PetscFunctionReturn(PETSC_SUCCESS);
228 }
229
MatSetRandom_ConstantDiagonal(Mat A,PetscRandom rand)230 static PetscErrorCode MatSetRandom_ConstantDiagonal(Mat A, PetscRandom rand)
231 {
232 PetscMPIInt rank;
233 MPI_Comm comm;
234 PetscScalar v = 0.0;
235 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;
236
237 PetscFunctionBegin;
238 PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
239 PetscCallMPI(MPI_Comm_rank(comm, &rank));
240 if (!rank) PetscCall(PetscRandomGetValue(rand, &v));
241 PetscCallMPI(MPI_Bcast(&v, 1, MPIU_SCALAR, 0, comm));
242 ctx->diag = v;
243 PetscFunctionReturn(PETSC_SUCCESS);
244 }
245
MatSolve_ConstantDiagonal(Mat matin,Vec b,Vec x)246 static PetscErrorCode MatSolve_ConstantDiagonal(Mat matin, Vec b, Vec x)
247 {
248 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)matin->data;
249
250 PetscFunctionBegin;
251 if (ctx->diag == 0.0) matin->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
252 else matin->factorerrortype = MAT_FACTOR_NOERROR;
253 PetscCall(VecAXPBY(x, 1.0 / ctx->diag, 0.0, b));
254 PetscFunctionReturn(PETSC_SUCCESS);
255 }
256
MatSOR_ConstantDiagonal(Mat matin,Vec x,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec y)257 static PetscErrorCode MatSOR_ConstantDiagonal(Mat matin, Vec x, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec y)
258 {
259 PetscFunctionBegin;
260 PetscCall(MatSolve_ConstantDiagonal(matin, x, y));
261 PetscFunctionReturn(PETSC_SUCCESS);
262 }
263
MatGetInfo_ConstantDiagonal(Mat A,MatInfoType flag,MatInfo * info)264 static PetscErrorCode MatGetInfo_ConstantDiagonal(Mat A, MatInfoType flag, MatInfo *info)
265 {
266 PetscFunctionBegin;
267 info->block_size = 1.0;
268 info->nz_allocated = 1.0;
269 info->nz_used = 1.0;
270 info->nz_unneeded = 0.0;
271 info->assemblies = A->num_ass;
272 info->mallocs = 0.0;
273 info->memory = 0; /* REVIEW ME */
274 if (A->factortype) {
275 info->fill_ratio_given = 1.0;
276 info->fill_ratio_needed = 1.0;
277 info->factor_mallocs = 0.0;
278 } else {
279 info->fill_ratio_given = 0;
280 info->fill_ratio_needed = 0;
281 info->factor_mallocs = 0;
282 }
283 PetscFunctionReturn(PETSC_SUCCESS);
284 }
285
286 /*@
287 MatCreateConstantDiagonal - Creates a matrix with a uniform value along the diagonal
288
289 Collective
290
291 Input Parameters:
292 + comm - MPI communicator
293 . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
294 This value should be the same as the local size used in creating the
295 y vector for the matrix-vector product y = Ax.
296 . n - This value should be the same as the local size used in creating the
297 x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have
298 calculated if `N` is given) For square matrices n is almost always `m`.
299 . M - number of global rows (or `PETSC_DETERMINE` to have calculated if m is given)
300 . N - number of global columns (or `PETSC_DETERMINE` to have calculated if n is given)
301 - diag - the diagonal value
302
303 Output Parameter:
304 . J - the diagonal matrix
305
306 Level: advanced
307
308 Notes:
309 Only supports square matrices with the same number of local rows and columns
310
311 .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`, `MatScale()`, `MatShift()`, `MatMult()`, `MatGetDiagonal()`, `MatGetFactor()`, `MatSolve()`
312 @*/
MatCreateConstantDiagonal(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscScalar diag,Mat * J)313 PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscScalar diag, Mat *J)
314 {
315 PetscFunctionBegin;
316 PetscCall(MatCreate(comm, J));
317 PetscCall(MatSetSizes(*J, m, n, M, N));
318 PetscCall(MatSetType(*J, MATCONSTANTDIAGONAL));
319 PetscCall(MatShift(*J, diag));
320 PetscCall(MatSetUp(*J));
321 PetscFunctionReturn(PETSC_SUCCESS);
322 }
323
324 /*@
325 MatConstantDiagonalGetConstant - Get the scalar constant of a constant diagonal matrix
326
327 Not collective
328
329 Input Parameter:
330 . mat - a `MATCONSTANTDIAGONAL`
331
332 Output Parameter:
333 . value - the scalar value
334
335 Level: developer
336
337 .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`
338 @*/
MatConstantDiagonalGetConstant(Mat mat,PetscScalar * value)339 PetscErrorCode MatConstantDiagonalGetConstant(Mat mat, PetscScalar *value)
340 {
341 PetscFunctionBegin;
342 PetscUseMethod(mat, "MatConstantDiagonalGetConstant_C", (Mat, PetscScalar *), (mat, value));
343 PetscFunctionReturn(PETSC_SUCCESS);
344 }
345
MatConstantDiagonalGetConstant_ConstantDiagonal(Mat mat,PetscScalar * value)346 static PetscErrorCode MatConstantDiagonalGetConstant_ConstantDiagonal(Mat mat, PetscScalar *value)
347 {
348 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data;
349
350 PetscFunctionBegin;
351 *value = ctx->diag;
352 PetscFunctionReturn(PETSC_SUCCESS);
353 }
354
355 /*MC
356 MATCONSTANTDIAGONAL - "constant-diagonal" - A diagonal matrix type with a uniform value
357 along the diagonal.
358
359 Level: advanced
360
361 .seealso: [](ch_matrices), `Mat`, `MatCreateConstantDiagonal()`
362 M*/
MatCreate_ConstantDiagonal(Mat A)363 PETSC_EXTERN PetscErrorCode MatCreate_ConstantDiagonal(Mat A)
364 {
365 Mat_ConstantDiagonal *ctx;
366
367 PetscFunctionBegin;
368 PetscCall(PetscNew(&ctx));
369 ctx->diag = 0.0;
370 A->data = (void *)ctx;
371
372 A->assembled = PETSC_TRUE;
373 A->preallocated = PETSC_TRUE;
374 A->structurally_symmetric = PETSC_BOOL3_TRUE;
375 A->structural_symmetry_eternal = PETSC_TRUE;
376 A->symmetric = PETSC_BOOL3_TRUE;
377 if (!PetscDefined(USE_COMPLEX)) A->hermitian = PETSC_BOOL3_TRUE;
378 A->symmetry_eternal = PETSC_TRUE;
379
380 A->ops->mult = MatMult_ConstantDiagonal;
381 A->ops->multadd = MatMultAdd_ConstantDiagonal;
382 A->ops->multtranspose = MatMult_ConstantDiagonal;
383 A->ops->multtransposeadd = MatMultAdd_ConstantDiagonal;
384 A->ops->multhermitiantranspose = MatMultHermitianTranspose_ConstantDiagonal;
385 A->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_ConstantDiagonal;
386 A->ops->solve = MatSolve_ConstantDiagonal;
387 A->ops->solvetranspose = MatSolve_ConstantDiagonal;
388 A->ops->norm = MatNorm_ConstantDiagonal;
389 A->ops->createsubmatrices = MatCreateSubMatrices_ConstantDiagonal;
390 A->ops->duplicate = MatDuplicate_ConstantDiagonal;
391 A->ops->getrow = MatGetRow_ConstantDiagonal;
392 A->ops->restorerow = MatRestoreRow_ConstantDiagonal;
393 A->ops->sor = MatSOR_ConstantDiagonal;
394 A->ops->shift = MatShift_ConstantDiagonal;
395 A->ops->scale = MatScale_ConstantDiagonal;
396 A->ops->getdiagonal = MatGetDiagonal_ConstantDiagonal;
397 A->ops->view = MatView_ConstantDiagonal;
398 A->ops->zeroentries = MatZeroEntries_ConstantDiagonal;
399 A->ops->destroy = MatDestroy_ConstantDiagonal;
400 A->ops->getinfo = MatGetInfo_ConstantDiagonal;
401 A->ops->equal = MatEqual_ConstantDiagonal;
402 A->ops->axpy = MatAXPY_ConstantDiagonal;
403 A->ops->setrandom = MatSetRandom_ConstantDiagonal;
404 A->ops->conjugate = MatConjugate_ConstantDiagonal;
405 A->ops->transpose = MatTranspose_ConstantDiagonal;
406
407 PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATCONSTANTDIAGONAL));
408 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConstantDiagonalGetConstant_C", MatConstantDiagonalGetConstant_ConstantDiagonal));
409 PetscFunctionReturn(PETSC_SUCCESS);
410 }
411
MatFactorNumeric_ConstantDiagonal(Mat fact,Mat A,const MatFactorInfo * info)412 static PetscErrorCode MatFactorNumeric_ConstantDiagonal(Mat fact, Mat A, const MatFactorInfo *info)
413 {
414 Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data, *fctx = (Mat_ConstantDiagonal *)fact->data;
415
416 PetscFunctionBegin;
417 if (actx->diag == 0.0) fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
418 else fact->factorerrortype = MAT_FACTOR_NOERROR;
419 fctx->diag = 1.0 / actx->diag;
420 fact->ops->solve = MatMult_ConstantDiagonal;
421 PetscFunctionReturn(PETSC_SUCCESS);
422 }
423
MatFactorSymbolic_LU_ConstantDiagonal(Mat fact,Mat A,IS isrow,IS iscol,const MatFactorInfo * info)424 static PetscErrorCode MatFactorSymbolic_LU_ConstantDiagonal(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
425 {
426 PetscFunctionBegin;
427 fact->ops->lufactornumeric = MatFactorNumeric_ConstantDiagonal;
428 PetscFunctionReturn(PETSC_SUCCESS);
429 }
430
MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact,Mat A,IS isrow,const MatFactorInfo * info)431 static PetscErrorCode MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact, Mat A, IS isrow, const MatFactorInfo *info)
432 {
433 PetscFunctionBegin;
434 fact->ops->choleskyfactornumeric = MatFactorNumeric_ConstantDiagonal;
435 PetscFunctionReturn(PETSC_SUCCESS);
436 }
437
MatGetFactor_constantdiagonal_petsc(Mat A,MatFactorType ftype,Mat * B)438 PETSC_INTERN PetscErrorCode MatGetFactor_constantdiagonal_petsc(Mat A, MatFactorType ftype, Mat *B)
439 {
440 PetscInt n = A->rmap->n, N = A->rmap->N;
441
442 PetscFunctionBegin;
443 PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), n, n, N, N, 0, B));
444
445 (*B)->factortype = ftype;
446 (*B)->ops->ilufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal;
447 (*B)->ops->lufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal;
448 (*B)->ops->iccfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;
449 (*B)->ops->choleskyfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;
450
451 (*B)->ops->shift = NULL;
452 (*B)->ops->scale = NULL;
453 (*B)->ops->mult = NULL;
454 (*B)->ops->sor = NULL;
455 (*B)->ops->zeroentries = NULL;
456
457 PetscCall(PetscFree((*B)->solvertype));
458 PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype));
459 PetscFunctionReturn(PETSC_SUCCESS);
460 }
461