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;
14: yctx->diag += a*xctx->diag;
15: return(0);
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;
21: PetscErrorCode ierr;
24: if (ncols) *ncols = 1;
25: if (cols) {
26: PetscMalloc1(1,cols);
27: (*cols)[0] = row;
28: }
29: if (vals) {
30: PetscMalloc1(1,vals);
31: (*vals)[0] = ctx->diag;
32: }
33: return(0);
34: }
36: static PetscErrorCode MatRestoreRow_ConstantDiagonal(Mat A, PetscInt row, PetscInt *ncols, PetscInt *cols[], PetscScalar *vals[])
37: {
41: if (ncols) *ncols = 0;
42: if (cols) {
43: PetscFree(*cols);
44: }
45: if (vals) {
46: PetscFree(*vals);
47: }
48: return(0);
49: }
51: static PetscErrorCode MatMultTranspose_ConstantDiagonal(Mat A, Vec x, Vec y)
52: {
53: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)A->data;
54: PetscErrorCode ierr;
57: VecAXPBY(y,ctx->diag,0.0,x);
58: return(0);
59: }
61: static PetscErrorCode MatMultAdd_ConstantDiagonal(Mat mat,Vec v1,Vec v2,Vec v3)
62: {
63: PetscErrorCode ierr;
64: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)mat->data;
67: if (v2 == v3) {
68: VecAXPBY(v3,ctx->diag,1.0,v1);
69: } else {
70: VecAXPBYPCZ(v3,ctx->diag,1.0,0.0,v1,v2);
71: }
72: return(0);
73: }
75: static PetscErrorCode MatMultTransposeAdd_ConstantDiagonal(Mat mat,Vec v1,Vec v2,Vec v3)
76: {
77: PetscErrorCode ierr;
78: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)mat->data;
81: if (v2 == v3) {
82: VecAXPBY(v3,ctx->diag,1.0,v1);
83: } else {
84: VecAXPBYPCZ(v3,ctx->diag,1.0,0.0,v1,v2);
85: }
86: return(0);
87: }
89: static PetscErrorCode MatDuplicate_ConstantDiagonal(Mat A, MatDuplicateOption op, Mat *B)
90: {
91: PetscErrorCode ierr;
92: Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal*)A->data;
95: MatCreate(PetscObjectComm((PetscObject)A),B);
96: MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
97: MatSetBlockSizesFromMats(*B,A,A);
98: MatSetType(*B,MATCONSTANTDIAGONAL);
99: PetscLayoutReference(A->rmap,&(*B)->rmap);
100: PetscLayoutReference(A->cmap,&(*B)->cmap);
101: if (op == MAT_COPY_VALUES) {
102: Mat_ConstantDiagonal *bctx = (Mat_ConstantDiagonal*)(*B)->data;
103: bctx->diag = actx->diag;
104: }
105: return(0);
106: }
108: static PetscErrorCode MatMissingDiagonal_ConstantDiagonal(Mat mat,PetscBool *missing,PetscInt *dd)
109: {
111: *missing = PETSC_FALSE;
112: return(0);
113: }
115: static PetscErrorCode MatDestroy_ConstantDiagonal(Mat mat)
116: {
117: PetscErrorCode ierr;
120: PetscFree(mat->data);
121: return(0);
122: }
124: static PetscErrorCode MatView_ConstantDiagonal(Mat J,PetscViewer viewer)
125: {
126: PetscErrorCode ierr;
127: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)J->data;
128: 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\n",ctx->diag);
139: #else
140: PetscViewerASCIIPrintf(viewer,"Diagonal value: %g + i %g\n",PetscRealPart(ctx->diag),PetscImaginaryPart(ctx->diag));
141: #endif
142: }
143: return(0);
144: }
146: static PetscErrorCode MatAssemblyEnd_ConstantDiagonal(Mat J,MatAssemblyType mt)
147: {
149: return(0);
150: }
152: static PetscErrorCode MatMult_ConstantDiagonal(Mat J,Vec x,Vec y)
153: {
154: PetscErrorCode ierr;
155: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)J->data;
158: VecAXPBY(y,ctx->diag,0.0,x);
159: return(0);
160: }
162: PetscErrorCode MatGetDiagonal_ConstantDiagonal(Mat J,Vec x)
163: {
164: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)J->data;
165: PetscErrorCode ierr;
168: VecSet(x,ctx->diag);
169: return(0);
170: }
172: static PetscErrorCode MatShift_ConstantDiagonal(Mat Y,PetscScalar a)
173: {
174: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)Y->data;
177: ctx->diag += a;
178: return(0);
179: }
181: static PetscErrorCode MatScale_ConstantDiagonal(Mat Y,PetscScalar a)
182: {
183: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)Y->data;
186: ctx->diag *= a;
187: return(0);
188: }
190: static PetscErrorCode MatZeroEntries_ConstantDiagonal(Mat Y)
191: {
192: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)Y->data;
195: ctx->diag = 0.0;
196: return(0);
197: }
199: PetscErrorCode MatSOR_ConstantDiagonal(Mat matin,Vec x,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec y)
200: {
201: PetscErrorCode ierr;
202: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)matin->data;
205: if (ctx->diag == 0.0) matin->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
206: else matin->factorerrortype = MAT_FACTOR_NOERROR;
207: VecAXPBY(y,1.0/ctx->diag,0.0,x);
208: return(0);
209: }
211: PetscErrorCode MatGetInfo_ConstantDiagonal(Mat A,MatInfoType flag,MatInfo *info)
212: {
214: info->block_size = 1.0;
215: info->nz_allocated = 1.0;
216: info->nz_used = 1.0;
217: info->nz_unneeded = 0.0;
218: info->assemblies = A->num_ass;
219: info->mallocs = 0.0;
220: info->memory = ((PetscObject)A)->mem;
221: if (A->factortype) {
222: info->fill_ratio_given = 1.0;
223: info->fill_ratio_needed = 1.0;
224: info->factor_mallocs = 0.0;
225: } else {
226: info->fill_ratio_given = 0;
227: info->fill_ratio_needed = 0;
228: info->factor_mallocs = 0;
229: }
230: return(0);
231: }
233: /*@
234: MatCreateConstantDiagonal - Creates a matrix with a uniform value along the diagonal
236: Collective
238: Input Parameters:
239: + comm - MPI communicator
240: . m - number of local rows (or PETSC_DECIDE to have calculated if M is given)
241: This value should be the same as the local size used in creating the
242: y vector for the matrix-vector product y = Ax.
243: . n - This value should be the same as the local size used in creating the
244: x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
245: calculated if N is given) For square matrices n is almost always m.
246: . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given)
247: . N - number of global columns (or PETSC_DETERMINE to have calculated if n is given)
248: - diag - the diagonal value
250: Output Parameter:
251: . J - the diagonal matrix
253: Level: advanced
255: Notes:
256: Only supports square matrices with the same number of local rows and columns
258: .seealso: MatDestroy(), MATCONSTANTDIAGONAL, MatScale(), MatShift(), MatMult(), MatGetDiagonal(), MatGetFactor(), MatSolve()
260: @*/
261: PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscScalar diag,Mat *J)
262: {
266: MatCreate(comm,J);
267: MatSetSizes(*J,m,n,M,N);
268: MatSetType(*J,MATCONSTANTDIAGONAL);
269: MatShift(*J,diag);
270: MatSetUp(*J);
271: return(0);
272: }
274: PETSC_EXTERN PetscErrorCode MatCreate_ConstantDiagonal(Mat A)
275: {
276: PetscErrorCode ierr;
277: Mat_ConstantDiagonal *ctx;
280: PetscNew(&ctx);
281: ctx->diag = 0.0;
282: A->data = (void*)ctx;
284: A->assembled = PETSC_TRUE;
285: A->preallocated = PETSC_TRUE;
287: A->ops->mult = MatMult_ConstantDiagonal;
288: A->ops->multadd = MatMultAdd_ConstantDiagonal;
289: A->ops->multtranspose = MatMultTranspose_ConstantDiagonal;
290: A->ops->multtransposeadd = MatMultTransposeAdd_ConstantDiagonal;
291: A->ops->duplicate = MatDuplicate_ConstantDiagonal;
292: A->ops->missingdiagonal = MatMissingDiagonal_ConstantDiagonal;
293: A->ops->getrow = MatGetRow_ConstantDiagonal;
294: A->ops->restorerow = MatRestoreRow_ConstantDiagonal;
295: A->ops->sor = MatSOR_ConstantDiagonal;
296: A->ops->shift = MatShift_ConstantDiagonal;
297: A->ops->scale = MatScale_ConstantDiagonal;
298: A->ops->getdiagonal = MatGetDiagonal_ConstantDiagonal;
299: A->ops->view = MatView_ConstantDiagonal;
300: A->ops->zeroentries = MatZeroEntries_ConstantDiagonal;
301: A->ops->assemblyend = MatAssemblyEnd_ConstantDiagonal;
302: A->ops->destroy = MatDestroy_ConstantDiagonal;
303: A->ops->getinfo = MatGetInfo_ConstantDiagonal;
304: A->ops->axpy = MatAXPY_ConstantDiagonal;
306: PetscObjectChangeTypeName((PetscObject)A,MATCONSTANTDIAGONAL);
307: return(0);
308: }
310: static PetscErrorCode MatFactorNumeric_ConstantDiagonal(Mat fact,Mat A,const MatFactorInfo *info)
311: {
312: Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal*)A->data,*fctx = (Mat_ConstantDiagonal*)fact->data;
315: if (actx->diag == 0.0) fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
316: else fact->factorerrortype = MAT_FACTOR_NOERROR;
317: fctx->diag = 1.0/actx->diag;
318: fact->ops->solve = MatMult_ConstantDiagonal;
319: return(0);
320: }
322: static PetscErrorCode MatFactorSymbolic_LU_ConstantDiagonal(Mat fact,Mat A,IS isrow,IS iscol,const MatFactorInfo *info)
323: {
325: fact->ops->lufactornumeric = MatFactorNumeric_ConstantDiagonal;
326: return(0);
327: }
329: static PetscErrorCode MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact,Mat A,IS isrow,const MatFactorInfo *info)
330: {
332: fact->ops->choleskyfactornumeric = MatFactorNumeric_ConstantDiagonal;
333: return(0);
334: }
336: PETSC_INTERN PetscErrorCode MatGetFactor_constantdiagonal_petsc(Mat A,MatFactorType ftype,Mat *B)
337: {
338: PetscInt n = A->rmap->n, N = A->rmap->N;
342: MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A),n,n,N,N,0,B);
344: (*B)->factortype = ftype;
345: (*B)->ops->ilufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal;
346: (*B)->ops->lufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal;
347: (*B)->ops->iccfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;
348: (*B)->ops->choleskyfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;
350: (*B)->ops->shift = NULL;
351: (*B)->ops->scale = NULL;
352: (*B)->ops->mult = NULL;
353: (*B)->ops->sor = NULL;
354: (*B)->ops->zeroentries = NULL;
356: PetscFree((*B)->solvertype);
357: PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype);
358: return(0);
359: }