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