Actual source code: cdiagonal.c

petsc-master 2019-07-16
Report Typos and Errors

  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   = (double)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: }