Actual source code: cdiagonal.c
petsc-3.12.5 2020-03-29
2: #include <petsc/private/matimpl.h>
4: typedef struct {
5: PetscScalar diag;
6: } Mat_ConstantDiagonal;
9: /* ----------------------------------------------------------------------------------------*/
10: static PetscErrorCode MatDestroy_ConstantDiagonal(Mat mat)
11: {
12: PetscErrorCode ierr;
15: PetscFree(mat->data);
16: return(0);
17: }
19: static PetscErrorCode MatView_ConstantDiagonal(Mat J,PetscViewer viewer)
20: {
21: PetscErrorCode ierr;
22: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)J->data;
23: PetscBool iascii;
26: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
27: if (iascii) {
28: PetscViewerFormat format;
30: PetscViewerGetFormat(viewer,&format);
31: if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) return(0);
32: #if !defined(PETSC_USE_COMPLEX)
33: PetscViewerASCIIPrintf(viewer,"Diagonal value: %g\n",ctx->diag);
34: #else
35: PetscViewerASCIIPrintf(viewer,"Diagonal value: %g + i %g\n",PetscRealPart(ctx->diag),PetscImaginaryPart(ctx->diag));
36: #endif
37: }
38: return(0);
39: }
41: static PetscErrorCode MatAssemblyEnd_ConstantDiagonal(Mat J,MatAssemblyType mt)
42: {
44: return(0);
45: }
47: static PetscErrorCode MatMult_ConstantDiagonal(Mat J,Vec x,Vec y)
48: {
49: PetscErrorCode ierr;
50: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)J->data;
53: VecAXPBY(y,ctx->diag,0.0,x);
54: return(0);
55: }
57: PetscErrorCode MatGetDiagonal_ConstantDiagonal(Mat J,Vec x)
58: {
59: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)J->data;
60: PetscErrorCode ierr;
63: VecSet(x,ctx->diag);
64: return(0);
65: }
67: static PetscErrorCode MatShift_ConstantDiagonal(Mat Y,PetscScalar a)
68: {
69: PetscErrorCode ierr;
70: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)Y->data;
73: if (a != 0.) {PetscObjectStateIncrease((PetscObject)Y);}
74: ctx->diag += a;
75: return(0);
76: }
78: static PetscErrorCode MatScale_ConstantDiagonal(Mat Y,PetscScalar a)
79: {
80: PetscErrorCode ierr;
81: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)Y->data;
84: if (a != 1.) {PetscObjectStateIncrease((PetscObject)Y);}
85: ctx->diag *= a;
86: return(0);
87: }
89: static PetscErrorCode MatZeroEntries_ConstantDiagonal(Mat Y)
90: {
91: PetscErrorCode ierr;
92: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)Y->data;
95: if (ctx->diag != 0.0) {PetscObjectStateIncrease((PetscObject)Y);}
96: ctx->diag = 0.0;
97: return(0);
98: }
100: PetscErrorCode MatSOR_ConstantDiagonal(Mat matin,Vec x,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec y)
101: {
102: PetscErrorCode ierr;
103: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)matin->data;
106: if (ctx->diag == 0.0) matin->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
107: else matin->factorerrortype = MAT_FACTOR_NOERROR;
108: VecAXPBY(y,1.0/ctx->diag,0.0,x);
109: return(0);
110: }
112: PetscErrorCode MatGetInfo_ConstantDiagonal(Mat A,MatInfoType flag,MatInfo *info)
113: {
115: info->block_size = 1.0;
116: info->nz_allocated = 1.0;
117: info->nz_used = 1.0;
118: info->nz_unneeded = 0.0;
119: info->assemblies = A->num_ass;
120: info->mallocs = 0.0;
121: info->memory = ((PetscObject)A)->mem;
122: if (A->factortype) {
123: info->fill_ratio_given = 1.0;
124: info->fill_ratio_needed = 1.0;
125: info->factor_mallocs = 0.0;
126: } else {
127: info->fill_ratio_given = 0;
128: info->fill_ratio_needed = 0;
129: info->factor_mallocs = 0;
130: }
131: return(0);
132: }
134: /*@
135: MatCreateConstantDiagonal - Creates a matrix with a uniform value along the diagonal
137: Collective
139: Input Parameters:
140: + comm - MPI communicator
141: . m - number of local rows (or PETSC_DECIDE to have calculated if M is given)
142: This value should be the same as the local size used in creating the
143: y vector for the matrix-vector product y = Ax.
144: . n - This value should be the same as the local size used in creating the
145: x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
146: calculated if N is given) For square matrices n is almost always m.
147: . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given)
148: . N - number of global columns (or PETSC_DETERMINE to have calculated if n is given)
149: - diag - the diagonal value
151: Output Parameter:
152: . J - the diagonal matrix
154: Level: advanced
156: Notes:
157: Only supports square matrices with the same number of local rows and columns
159: .seealso: MatDestroy(), MATCONSTANTDIAGONAL, MatScale(), MatShift(), MatMult(), MatGetDiagonal(), MatGetFactor(), MatSolve()
161: @*/
162: PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscScalar diag,Mat *J)
163: {
167: MatCreate(comm,J);
168: MatSetSizes(*J,m,n,M,N);
169: MatSetType(*J,MATCONSTANTDIAGONAL);
170: MatShift(*J,diag);
171: MatSetUp(*J);
172: return(0);
173: }
175: PETSC_EXTERN PetscErrorCode MatCreate_ConstantDiagonal(Mat A)
176: {
177: PetscErrorCode ierr;
178: Mat_ConstantDiagonal *ctx;
181: PetscNew(&ctx);
182: ctx->diag = 0.0;
183: A->data = (void*)ctx;
185: A->assembled = PETSC_TRUE;
186: A->preallocated = PETSC_TRUE;
187: A->ops->mult = MatMult_ConstantDiagonal;
188: A->ops->sor = MatSOR_ConstantDiagonal;
189: A->ops->shift = MatShift_ConstantDiagonal;
190: A->ops->scale = MatScale_ConstantDiagonal;
191: A->ops->getdiagonal = MatGetDiagonal_ConstantDiagonal;
192: A->ops->view = MatView_ConstantDiagonal;
193: A->ops->zeroentries = MatZeroEntries_ConstantDiagonal;
194: A->ops->assemblyend = MatAssemblyEnd_ConstantDiagonal;
195: A->ops->destroy = MatDestroy_ConstantDiagonal;
196: A->ops->getinfo = MatGetInfo_ConstantDiagonal;
197: PetscObjectChangeTypeName((PetscObject)A,MATCONSTANTDIAGONAL);
198: return(0);
199: }
201: static PetscErrorCode MatFactorNumeric_ConstantDiagonal(Mat fact,Mat A,const MatFactorInfo *info)
202: {
203: PetscErrorCode ierr;
204: Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal*)A->data,*fctx = (Mat_ConstantDiagonal*)fact->data;
207: if (actx->diag == 0.0) fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
208: else fact->factorerrortype = MAT_FACTOR_NOERROR;
209: fctx->diag = 1.0/actx->diag;
210: PetscObjectStateIncrease((PetscObject)fact);
211: fact->ops->solve = MatMult_ConstantDiagonal;
212: return(0);
213: }
215: static PetscErrorCode MatFactorSymbolic_LU_ConstantDiagonal(Mat fact,Mat A,IS isrow,IS iscol,const MatFactorInfo *info)
216: {
218: fact->ops->lufactornumeric = MatFactorNumeric_ConstantDiagonal;
219: return(0);
220: }
222: static PetscErrorCode MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact,Mat A,IS isrow,const MatFactorInfo *info)
223: {
225: fact->ops->choleskyfactornumeric = MatFactorNumeric_ConstantDiagonal; return(0);
226: return(0);
227: }
229: PETSC_INTERN PetscErrorCode MatGetFactor_constantdiagonal_petsc(Mat A,MatFactorType ftype,Mat *B)
230: {
231: PetscInt n = A->rmap->n, N = A->rmap->N;
235: MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A),n,n,N,N,0,B);
237: (*B)->factortype = ftype;
238: (*B)->ops->ilufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal;
239: (*B)->ops->lufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal;
240: (*B)->ops->iccfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;
241: (*B)->ops->choleskyfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;
243: (*B)->ops->shift = NULL;
244: (*B)->ops->scale = NULL;
245: (*B)->ops->mult = NULL;
246: (*B)->ops->sor = NULL;
247: (*B)->ops->zeroentries = NULL;
249: PetscFree((*B)->solvertype);
250: PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype);
251: return(0);
252: }