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