Actual source code: cdiagonal.c

petsc-3.14.6 2021-03-30
Report Typos and Errors

  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: }