Actual source code: pcis.c

      2: #include <../src/ksp/pc/impls/is/pcis.h>

      4: static PetscErrorCode PCISSetUseStiffnessScaling_IS(PC pc, PetscBool use)
      5: {
      6:   PC_IS *pcis = (PC_IS*)pc->data;

      9:   pcis->use_stiffness_scaling = use;
     10:   return(0);
     11: }

     13: /*@
     14:  PCISSetUseStiffnessScaling - Tells PCIS to construct partition of unity using
     15:                               local matrices' diagonal.

     17:    Not collective

     19:    Input Parameters:
     20: +  pc - the preconditioning context
     21: -  use - whether or not pcis use matrix diagonal to build partition of unity.

     23:    Level: intermediate

     25:    Notes:

     27: .seealso: PCBDDC
     28: @*/
     29: PetscErrorCode PCISSetUseStiffnessScaling(PC pc, PetscBool use)
     30: {

     36:   PetscTryMethod(pc,"PCISSetUseStiffnessScaling_C",(PC,PetscBool),(pc,use));
     37:   return(0);
     38: }

     40: static PetscErrorCode PCISSetSubdomainDiagonalScaling_IS(PC pc, Vec scaling_factors)
     41: {
     43:   PC_IS          *pcis = (PC_IS*)pc->data;

     46:   PetscObjectReference((PetscObject)scaling_factors);
     47:   VecDestroy(&pcis->D);
     48:   pcis->D = scaling_factors;
     49:   if (pc->setupcalled) {
     50:     PetscInt sn;

     52:     VecGetSize(pcis->D,&sn);
     53:     if (sn == pcis->n) {
     54:       VecScatterBegin(pcis->N_to_B,pcis->D,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
     55:       VecScatterEnd(pcis->N_to_B,pcis->D,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
     56:       VecDestroy(&pcis->D);
     57:       VecDuplicate(pcis->vec1_B,&pcis->D);
     58:       VecCopy(pcis->vec1_B,pcis->D);
     59:     } else if (sn != pcis->n_B) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Invalid size for scaling vector. Expected %D (or full %D), found %D",pcis->n_B,pcis->n,sn);
     60:   }
     61:   return(0);
     62: }

     64: /*@
     65:  PCISSetSubdomainDiagonalScaling - Set diagonal scaling for PCIS.

     67:    Not collective

     69:    Input Parameters:
     70: +  pc - the preconditioning context
     71: -  scaling_factors - scaling factors for the subdomain

     73:    Level: intermediate

     75:    Notes:
     76:    Intended to use with jumping coefficients cases.

     78: .seealso: PCBDDC
     79: @*/
     80: PetscErrorCode PCISSetSubdomainDiagonalScaling(PC pc, Vec scaling_factors)
     81: {

     87:   PetscTryMethod(pc,"PCISSetSubdomainDiagonalScaling_C",(PC,Vec),(pc,scaling_factors));
     88:   return(0);
     89: }

     91: static PetscErrorCode PCISSetSubdomainScalingFactor_IS(PC pc, PetscScalar scal)
     92: {
     93:   PC_IS *pcis = (PC_IS*)pc->data;

     96:   pcis->scaling_factor = scal;
     97:   if (pcis->D) {

    100:     VecSet(pcis->D,pcis->scaling_factor);
    101:   }
    102:   return(0);
    103: }

    105: /*@
    106:  PCISSetSubdomainScalingFactor - Set scaling factor for PCIS.

    108:    Not collective

    110:    Input Parameters:
    111: +  pc - the preconditioning context
    112: -  scal - scaling factor for the subdomain

    114:    Level: intermediate

    116:    Notes:
    117:    Intended to use with jumping coefficients cases.

    119: .seealso: PCBDDC
    120: @*/
    121: PetscErrorCode PCISSetSubdomainScalingFactor(PC pc, PetscScalar scal)
    122: {

    127:   PetscTryMethod(pc,"PCISSetSubdomainScalingFactor_C",(PC,PetscScalar),(pc,scal));
    128:   return(0);
    129: }

    132: /* -------------------------------------------------------------------------- */
    133: /*
    134:    PCISSetUp -
    135: */
    136: PetscErrorCode  PCISSetUp(PC pc, PetscBool computematrices, PetscBool computesolvers)
    137: {
    138:   PC_IS          *pcis  = (PC_IS*)(pc->data);
    139:   Mat_IS         *matis;
    140:   MatReuse       reuse;
    142:   PetscBool      flg,issbaij;

    145:   PetscObjectTypeCompare((PetscObject)pc->pmat,MATIS,&flg);
    146:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Requires preconditioning matrix of type MATIS");
    147:   matis = (Mat_IS*)pc->pmat->data;
    148:   if (pc->useAmat) {
    149:     PetscObjectTypeCompare((PetscObject)pc->mat,MATIS,&flg);
    150:     if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Requires linear system matrix of type MATIS");
    151:   }

    153:   /* first time creation, get info on substructuring */
    154:   if (!pc->setupcalled) {
    155:     PetscInt    n_I;
    156:     PetscInt    *idx_I_local,*idx_B_local,*idx_I_global,*idx_B_global;
    157:     PetscBT     bt;
    158:     PetscInt    i,j;

    160:     /* get info on mapping */
    161:     PetscObjectReference((PetscObject)pc->pmat->rmap->mapping);
    162:     ISLocalToGlobalMappingDestroy(&pcis->mapping);
    163:     pcis->mapping = pc->pmat->rmap->mapping;
    164:     ISLocalToGlobalMappingGetSize(pcis->mapping,&pcis->n);
    165:     ISLocalToGlobalMappingGetInfo(pcis->mapping,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));

    167:     /* Identifying interior and interface nodes, in local numbering */
    168:     PetscBTCreate(pcis->n,&bt);
    169:     for (i=0;i<pcis->n_neigh;i++)
    170:       for (j=0;j<pcis->n_shared[i];j++) {
    171:         PetscBTSet(bt,pcis->shared[i][j]);
    172:       }

    174:     /* Creating local and global index sets for interior and inteface nodes. */
    175:     PetscMalloc1(pcis->n,&idx_I_local);
    176:     PetscMalloc1(pcis->n,&idx_B_local);
    177:     for (i=0, pcis->n_B=0, n_I=0; i<pcis->n; i++) {
    178:       if (!PetscBTLookup(bt,i)) {
    179:         idx_I_local[n_I] = i;
    180:         n_I++;
    181:       } else {
    182:         idx_B_local[pcis->n_B] = i;
    183:         pcis->n_B++;
    184:       }
    185:     }

    187:     /* Getting the global numbering */
    188:     idx_B_global = idx_I_local + n_I; /* Just avoiding allocating extra memory, since we have vacant space */
    189:     idx_I_global = idx_B_local + pcis->n_B;
    190:     ISLocalToGlobalMappingApply(pcis->mapping,pcis->n_B,idx_B_local,idx_B_global);
    191:     ISLocalToGlobalMappingApply(pcis->mapping,n_I,idx_I_local,idx_I_global);

    193:     /* Creating the index sets */
    194:     ISCreateGeneral(PETSC_COMM_SELF,pcis->n_B,idx_B_local,PETSC_COPY_VALUES, &pcis->is_B_local);
    195:     ISCreateGeneral(PetscObjectComm((PetscObject)pc),pcis->n_B,idx_B_global,PETSC_COPY_VALUES,&pcis->is_B_global);
    196:     ISCreateGeneral(PETSC_COMM_SELF,n_I,idx_I_local,PETSC_COPY_VALUES, &pcis->is_I_local);
    197:     ISCreateGeneral(PetscObjectComm((PetscObject)pc),n_I,idx_I_global,PETSC_COPY_VALUES,&pcis->is_I_global);

    199:     /* Freeing memory */
    200:     PetscFree(idx_B_local);
    201:     PetscFree(idx_I_local);
    202:     PetscBTDestroy(&bt);

    204:     /* Creating work vectors and arrays */
    205:     VecDuplicate(matis->x,&pcis->vec1_N);
    206:     VecDuplicate(pcis->vec1_N,&pcis->vec2_N);
    207:     VecCreate(PETSC_COMM_SELF,&pcis->vec1_D);
    208:     VecSetSizes(pcis->vec1_D,pcis->n-pcis->n_B,PETSC_DECIDE);
    209:     VecSetType(pcis->vec1_D,((PetscObject)pcis->vec1_N)->type_name);
    210:     VecDuplicate(pcis->vec1_D,&pcis->vec2_D);
    211:     VecDuplicate(pcis->vec1_D,&pcis->vec3_D);
    212:     VecDuplicate(pcis->vec1_D,&pcis->vec4_D);
    213:     VecCreate(PETSC_COMM_SELF,&pcis->vec1_B);
    214:     VecSetSizes(pcis->vec1_B,pcis->n_B,PETSC_DECIDE);
    215:     VecSetType(pcis->vec1_B,((PetscObject)pcis->vec1_N)->type_name);
    216:     VecDuplicate(pcis->vec1_B,&pcis->vec2_B);
    217:     VecDuplicate(pcis->vec1_B,&pcis->vec3_B);
    218:     MatCreateVecs(pc->pmat,&pcis->vec1_global,NULL);
    219:     PetscMalloc1(pcis->n,&pcis->work_N);
    220:     /* scaling vector */
    221:     if (!pcis->D) { /* it can happen that the user passed in a scaling vector via PCISSetSubdomainDiagonalScaling */
    222:       VecDuplicate(pcis->vec1_B,&pcis->D);
    223:       VecSet(pcis->D,pcis->scaling_factor);
    224:     }

    226:     /* Creating the scatter contexts */
    227:     VecScatterCreate(pcis->vec1_N,pcis->is_I_local,pcis->vec1_D,(IS)0,&pcis->N_to_D);
    228:     VecScatterCreate(pcis->vec1_global,pcis->is_I_global,pcis->vec1_D,(IS)0,&pcis->global_to_D);
    229:     VecScatterCreate(pcis->vec1_N,pcis->is_B_local,pcis->vec1_B,(IS)0,&pcis->N_to_B);
    230:     VecScatterCreate(pcis->vec1_global,pcis->is_B_global,pcis->vec1_B,(IS)0,&pcis->global_to_B);

    232:     /* map from boundary to local */
    233:     ISLocalToGlobalMappingCreateIS(pcis->is_B_local,&pcis->BtoNmap);
    234:   }

    236:   {
    237:     PetscInt sn;

    239:     VecGetSize(pcis->D,&sn);
    240:     if (sn == pcis->n) {
    241:       VecScatterBegin(pcis->N_to_B,pcis->D,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
    242:       VecScatterEnd(pcis->N_to_B,pcis->D,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
    243:       VecDestroy(&pcis->D);
    244:       VecDuplicate(pcis->vec1_B,&pcis->D);
    245:       VecCopy(pcis->vec1_B,pcis->D);
    246:     } else if (sn != pcis->n_B) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Invalid size for scaling vector. Expected %D (or full %D), found %D",pcis->n_B,pcis->n,sn);
    247:   }

    249:   /*
    250:     Extracting the blocks A_II, A_BI, A_IB and A_BB from A. If the numbering
    251:     is such that interior nodes come first than the interface ones, we have

    253:         [ A_II | A_IB ]
    254:     A = [------+------]
    255:         [ A_BI | A_BB ]
    256:   */
    257:   if (computematrices) {
    258:     PetscBool amat = (PetscBool)(pc->mat != pc->pmat && pc->useAmat);
    259:     PetscInt  bs,ibs;

    261:     reuse = MAT_INITIAL_MATRIX;
    262:     if (pcis->reusesubmatrices && pc->setupcalled) {
    263:       if (pc->flag == SAME_NONZERO_PATTERN) {
    264:         reuse = MAT_REUSE_MATRIX;
    265:       } else {
    266:         reuse = MAT_INITIAL_MATRIX;
    267:       }
    268:     }
    269:     if (reuse == MAT_INITIAL_MATRIX) {
    270:       MatDestroy(&pcis->A_II);
    271:       MatDestroy(&pcis->pA_II);
    272:       MatDestroy(&pcis->A_IB);
    273:       MatDestroy(&pcis->A_BI);
    274:       MatDestroy(&pcis->A_BB);
    275:     }

    277:     ISLocalToGlobalMappingGetBlockSize(pcis->mapping,&ibs);
    278:     MatGetBlockSize(matis->A,&bs);
    279:     MatCreateSubMatrix(matis->A,pcis->is_I_local,pcis->is_I_local,reuse,&pcis->pA_II);
    280:     if (amat) {
    281:       Mat_IS *amatis = (Mat_IS*)pc->mat->data;
    282:       MatCreateSubMatrix(amatis->A,pcis->is_I_local,pcis->is_I_local,reuse,&pcis->A_II);
    283:     } else {
    284:       PetscObjectReference((PetscObject)pcis->pA_II);
    285:       MatDestroy(&pcis->A_II);
    286:       pcis->A_II = pcis->pA_II;
    287:     }
    288:     MatSetBlockSize(pcis->A_II,bs == ibs ? bs : 1);
    289:     MatSetBlockSize(pcis->pA_II,bs == ibs ? bs : 1);
    290:     MatCreateSubMatrix(matis->A,pcis->is_B_local,pcis->is_B_local,reuse,&pcis->A_BB);
    291:     PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);
    292:     if (!issbaij) {
    293:       MatCreateSubMatrix(matis->A,pcis->is_I_local,pcis->is_B_local,reuse,&pcis->A_IB);
    294:       MatCreateSubMatrix(matis->A,pcis->is_B_local,pcis->is_I_local,reuse,&pcis->A_BI);
    295:     } else {
    296:       Mat newmat;

    298:       MatConvert(matis->A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&newmat);
    299:       MatCreateSubMatrix(newmat,pcis->is_I_local,pcis->is_B_local,reuse,&pcis->A_IB);
    300:       MatCreateSubMatrix(newmat,pcis->is_B_local,pcis->is_I_local,reuse,&pcis->A_BI);
    301:       MatDestroy(&newmat);
    302:     }
    303:     MatSetBlockSize(pcis->A_BB,bs == ibs ? bs : 1);
    304:   }

    306:   /* Creating scaling vector D */
    307:   PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_is_use_stiffness_scaling",&pcis->use_stiffness_scaling,NULL);
    308:   if (pcis->use_stiffness_scaling) {
    309:     PetscScalar *a;
    310:     PetscInt    i,n;

    312:     if (pcis->A_BB) {
    313:       MatGetDiagonal(pcis->A_BB,pcis->D);
    314:     } else {
    315:       MatGetDiagonal(matis->A,pcis->vec1_N);
    316:       VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);
    317:       VecScatterEnd(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);
    318:     }
    319:     VecAbs(pcis->D);
    320:     VecGetLocalSize(pcis->D,&n);
    321:     VecGetArray(pcis->D,&a);
    322:     for (i=0;i<n;i++) if (PetscAbsScalar(a[i])<PETSC_SMALL) a[i] = 1.0;
    323:     VecRestoreArray(pcis->D,&a);
    324:   }
    325:   VecSet(pcis->vec1_global,0.0);
    326:   VecScatterBegin(pcis->global_to_B,pcis->D,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
    327:   VecScatterEnd(pcis->global_to_B,pcis->D,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
    328:   VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
    329:   VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
    330:   VecPointwiseDivide(pcis->D,pcis->D,pcis->vec1_B);
    331:   /* See historical note 01, at the bottom of this file. */

    333:   /* Creating the KSP contexts for the local Dirichlet and Neumann problems */
    334:   if (computesolvers) {
    335:     PC pc_ctx;

    337:     pcis->pure_neumann = matis->pure_neumann;
    338:     /* Dirichlet */
    339:     KSPCreate(PETSC_COMM_SELF,&pcis->ksp_D);
    340:     KSPSetErrorIfNotConverged(pcis->ksp_D,pc->erroriffailure);
    341:     PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_D,(PetscObject)pc,1);
    342:     KSPSetOperators(pcis->ksp_D,pcis->A_II,pcis->A_II);
    343:     KSPSetOptionsPrefix(pcis->ksp_D,"is_localD_");
    344:     KSPGetPC(pcis->ksp_D,&pc_ctx);
    345:     PCSetType(pc_ctx,PCLU);
    346:     KSPSetType(pcis->ksp_D,KSPPREONLY);
    347:     KSPSetFromOptions(pcis->ksp_D);
    348:     /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
    349:     KSPSetUp(pcis->ksp_D);
    350:     /* Neumann */
    351:     KSPCreate(PETSC_COMM_SELF,&pcis->ksp_N);
    352:     KSPSetErrorIfNotConverged(pcis->ksp_N,pc->erroriffailure);
    353:     PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_N,(PetscObject)pc,1);
    354:     KSPSetOperators(pcis->ksp_N,matis->A,matis->A);
    355:     KSPSetOptionsPrefix(pcis->ksp_N,"is_localN_");
    356:     KSPGetPC(pcis->ksp_N,&pc_ctx);
    357:     PCSetType(pc_ctx,PCLU);
    358:     KSPSetType(pcis->ksp_N,KSPPREONLY);
    359:     KSPSetFromOptions(pcis->ksp_N);
    360:     {
    361:       PetscBool damp_fixed                    = PETSC_FALSE,
    362:                 remove_nullspace_fixed        = PETSC_FALSE,
    363:                 set_damping_factor_floating   = PETSC_FALSE,
    364:                 not_damp_floating             = PETSC_FALSE,
    365:                 not_remove_nullspace_floating = PETSC_FALSE;
    366:       PetscReal fixed_factor,
    367:                 floating_factor;

    369:       PetscOptionsGetReal(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&fixed_factor,&damp_fixed);
    370:       if (!damp_fixed) fixed_factor = 0.0;
    371:       PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&damp_fixed,NULL);

    373:       PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_remove_nullspace_fixed",&remove_nullspace_fixed,NULL);

    375:       PetscOptionsGetReal(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_set_damping_factor_floating",
    376:                               &floating_factor,&set_damping_factor_floating);
    377:       if (!set_damping_factor_floating) floating_factor = 0.0;
    378:       PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_set_damping_factor_floating",&set_damping_factor_floating,NULL);
    379:       if (!set_damping_factor_floating) floating_factor = 1.e-12;

    381:       PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_not_damp_floating",&not_damp_floating,NULL);

    383:       PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_not_remove_nullspace_floating",&not_remove_nullspace_floating,NULL);

    385:       if (pcis->pure_neumann) {  /* floating subdomain */
    386:         if (!(not_damp_floating)) {
    387:           PCFactorSetShiftType(pc_ctx,MAT_SHIFT_NONZERO);
    388:           PCFactorSetShiftAmount(pc_ctx,floating_factor);
    389:         }
    390:         if (!(not_remove_nullspace_floating)) {
    391:           MatNullSpace nullsp;
    392:           MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_TRUE,0,NULL,&nullsp);
    393:           MatSetNullSpace(matis->A,nullsp);
    394:           MatNullSpaceDestroy(&nullsp);
    395:         }
    396:       } else {  /* fixed subdomain */
    397:         if (damp_fixed) {
    398:           PCFactorSetShiftType(pc_ctx,MAT_SHIFT_NONZERO);
    399:           PCFactorSetShiftAmount(pc_ctx,floating_factor);
    400:         }
    401:         if (remove_nullspace_fixed) {
    402:           MatNullSpace nullsp;
    403:           MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_TRUE,0,NULL,&nullsp);
    404:           MatSetNullSpace(matis->A,nullsp);
    405:           MatNullSpaceDestroy(&nullsp);
    406:         }
    407:       }
    408:     }
    409:     /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
    410:     KSPSetUp(pcis->ksp_N);
    411:   }
    412:   return(0);
    413: }

    415: /* -------------------------------------------------------------------------- */
    416: /*
    417:    PCISDestroy -
    418: */
    419: PetscErrorCode  PCISDestroy(PC pc)
    420: {
    421:   PC_IS          *pcis = (PC_IS*)(pc->data);

    425:   ISDestroy(&pcis->is_B_local);
    426:   ISDestroy(&pcis->is_I_local);
    427:   ISDestroy(&pcis->is_B_global);
    428:   ISDestroy(&pcis->is_I_global);
    429:   MatDestroy(&pcis->A_II);
    430:   MatDestroy(&pcis->pA_II);
    431:   MatDestroy(&pcis->A_IB);
    432:   MatDestroy(&pcis->A_BI);
    433:   MatDestroy(&pcis->A_BB);
    434:   VecDestroy(&pcis->D);
    435:   KSPDestroy(&pcis->ksp_N);
    436:   KSPDestroy(&pcis->ksp_D);
    437:   VecDestroy(&pcis->vec1_N);
    438:   VecDestroy(&pcis->vec2_N);
    439:   VecDestroy(&pcis->vec1_D);
    440:   VecDestroy(&pcis->vec2_D);
    441:   VecDestroy(&pcis->vec3_D);
    442:   VecDestroy(&pcis->vec4_D);
    443:   VecDestroy(&pcis->vec1_B);
    444:   VecDestroy(&pcis->vec2_B);
    445:   VecDestroy(&pcis->vec3_B);
    446:   VecDestroy(&pcis->vec1_global);
    447:   VecScatterDestroy(&pcis->global_to_D);
    448:   VecScatterDestroy(&pcis->N_to_B);
    449:   VecScatterDestroy(&pcis->N_to_D);
    450:   VecScatterDestroy(&pcis->global_to_B);
    451:   PetscFree(pcis->work_N);
    452:   if (pcis->n_neigh > -1) {
    453:     ISLocalToGlobalMappingRestoreInfo(pcis->mapping,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));
    454:   }
    455:   ISLocalToGlobalMappingDestroy(&pcis->mapping);
    456:   ISLocalToGlobalMappingDestroy(&pcis->BtoNmap);
    457:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetUseStiffnessScaling_C",NULL);
    458:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainScalingFactor_C",NULL);
    459:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainDiagonalScaling_C",NULL);
    460:   return(0);
    461: }

    463: /* -------------------------------------------------------------------------- */
    464: /*
    465:    PCISCreate -
    466: */
    467: PetscErrorCode  PCISCreate(PC pc)
    468: {
    469:   PC_IS          *pcis = (PC_IS*)(pc->data);

    473:   pcis->n_neigh          = -1;
    474:   pcis->scaling_factor   = 1.0;
    475:   pcis->reusesubmatrices = PETSC_TRUE;
    476:   /* composing functions */
    477:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetUseStiffnessScaling_C",PCISSetUseStiffnessScaling_IS);
    478:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainScalingFactor_C",PCISSetSubdomainScalingFactor_IS);
    479:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainDiagonalScaling_C",PCISSetSubdomainDiagonalScaling_IS);
    480:   return(0);
    481: }

    483: /* -------------------------------------------------------------------------- */
    484: /*
    485:    PCISApplySchur -

    487:    Input parameters:
    488: .  pc - preconditioner context
    489: .  v - vector to which the Schur complement is to be applied (it is NOT modified inside this function, UNLESS vec2_B is null)

    491:    Output parameters:
    492: .  vec1_B - result of Schur complement applied to chunk
    493: .  vec2_B - garbage (used as work space), or null (and v is used as workspace)
    494: .  vec1_D - garbage (used as work space)
    495: .  vec2_D - garbage (used as work space)

    497: */
    498: PetscErrorCode  PCISApplySchur(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D)
    499: {
    501:   PC_IS          *pcis = (PC_IS*)(pc->data);

    504:   if (!vec2_B) vec2_B = v;

    506:   MatMult(pcis->A_BB,v,vec1_B);
    507:   MatMult(pcis->A_IB,v,vec1_D);
    508:   KSPSolve(pcis->ksp_D,vec1_D,vec2_D);
    509:   KSPCheckSolve(pcis->ksp_D,pc,vec2_D);
    510:   MatMult(pcis->A_BI,vec2_D,vec2_B);
    511:   VecAXPY(vec1_B,-1.0,vec2_B);
    512:   return(0);
    513: }

    515: /* -------------------------------------------------------------------------- */
    516: /*
    517:    PCISScatterArrayNToVecB - Scatters interface node values from a big array (of all local nodes, interior or interface,
    518:    including ghosts) into an interface vector, when in SCATTER_FORWARD mode, or vice-versa, when in SCATTER_REVERSE
    519:    mode.

    521:    Input parameters:
    522: .  pc - preconditioner context
    523: .  array_N - [when in SCATTER_FORWARD mode] Array to be scattered into the vector
    524: .  v_B - [when in SCATTER_REVERSE mode] Vector to be scattered into the array

    526:    Output parameter:
    527: .  array_N - [when in SCATTER_REVERSE mode] Array to receive the scattered vector
    528: .  v_B - [when in SCATTER_FORWARD mode] Vector to receive the scattered array

    530:    Notes:
    531:    The entries in the array that do not correspond to interface nodes remain unaltered.
    532: */
    533: PetscErrorCode  PCISScatterArrayNToVecB(PetscScalar *array_N, Vec v_B, InsertMode imode, ScatterMode smode, PC pc)
    534: {
    535:   PetscInt       i;
    536:   const PetscInt *idex;
    538:   PetscScalar    *array_B;
    539:   PC_IS          *pcis = (PC_IS*)(pc->data);

    542:   VecGetArray(v_B,&array_B);
    543:   ISGetIndices(pcis->is_B_local,&idex);

    545:   if (smode == SCATTER_FORWARD) {
    546:     if (imode == INSERT_VALUES) {
    547:       for (i=0; i<pcis->n_B; i++) array_B[i] = array_N[idex[i]];
    548:     } else {  /* ADD_VALUES */
    549:       for (i=0; i<pcis->n_B; i++) array_B[i] += array_N[idex[i]];
    550:     }
    551:   } else {  /* SCATTER_REVERSE */
    552:     if (imode == INSERT_VALUES) {
    553:       for (i=0; i<pcis->n_B; i++) array_N[idex[i]] = array_B[i];
    554:     } else {  /* ADD_VALUES */
    555:       for (i=0; i<pcis->n_B; i++) array_N[idex[i]] += array_B[i];
    556:     }
    557:   }
    558:   ISRestoreIndices(pcis->is_B_local,&idex);
    559:   VecRestoreArray(v_B,&array_B);
    560:   return(0);
    561: }

    563: /* -------------------------------------------------------------------------- */
    564: /*
    565:    PCISApplyInvSchur - Solves the Neumann problem related to applying the inverse of the Schur complement.
    566:    More precisely, solves the problem:
    567:                                         [ A_II  A_IB ] [ . ]   [ 0 ]
    568:                                         [            ] [   ] = [   ]
    569:                                         [ A_BI  A_BB ] [ x ]   [ b ]

    571:    Input parameters:
    572: .  pc - preconditioner context
    573: .  b - vector of local interface nodes (including ghosts)

    575:    Output parameters:
    576: .  x - vector of local interface nodes (including ghosts); returns the application of the inverse of the Schur
    577:        complement to b
    578: .  vec1_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space)
    579: .  vec2_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space)

    581: */
    582: PetscErrorCode  PCISApplyInvSchur(PC pc, Vec b, Vec x, Vec vec1_N, Vec vec2_N)
    583: {
    585:   PC_IS          *pcis = (PC_IS*)(pc->data);

    588:   /*
    589:     Neumann solvers.
    590:     Applying the inverse of the local Schur complement, i.e, solving a Neumann
    591:     Problem with zero at the interior nodes of the RHS and extracting the interface
    592:     part of the solution. inverse Schur complement is applied to b and the result
    593:     is stored in x.
    594:   */
    595:   /* Setting the RHS vec1_N */
    596:   VecSet(vec1_N,0.0);
    597:   VecScatterBegin(pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);
    598:   VecScatterEnd  (pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);
    599:   /* Checking for consistency of the RHS */
    600:   {
    601:     PetscBool flg = PETSC_FALSE;
    602:     PetscOptionsGetBool(NULL,NULL,"-pc_is_check_consistency",&flg,NULL);
    603:     if (flg) {
    604:       PetscScalar average;
    605:       PetscViewer viewer;
    606:       PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);

    608:       VecSum(vec1_N,&average);
    609:       average = average / ((PetscReal)pcis->n);
    610:       PetscViewerASCIIPushSynchronized(viewer);
    611:       if (pcis->pure_neumann) {
    612:         PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is floating. Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));
    613:       } else {
    614:         PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is fixed.    Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));
    615:       }
    616:       PetscViewerFlush(viewer);
    617:       PetscViewerASCIIPopSynchronized(viewer);
    618:     }
    619:   }
    620:   /* Solving the system for vec2_N */
    621:   KSPSolve(pcis->ksp_N,vec1_N,vec2_N);
    622:   KSPCheckSolve(pcis->ksp_N,pc,vec2_N);
    623:   /* Extracting the local interface vector out of the solution */
    624:   VecScatterBegin(pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);
    625:   VecScatterEnd  (pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);
    626:   return(0);
    627: }