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

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

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

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

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

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

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

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

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

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

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

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

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

235:   {
236:     PetscInt sn;

238:     VecGetSize(pcis->D,&sn);
239:     if (sn == pcis->n) {
240:       VecScatterBegin(pcis->N_to_B,pcis->D,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
241:       VecScatterEnd(pcis->N_to_B,pcis->D,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
242:       VecDestroy(&pcis->D);
243:       VecDuplicate(pcis->vec1_B,&pcis->D);
244:       VecCopy(pcis->vec1_B,pcis->D);
245:     } 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);
246:   }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

482: /* -------------------------------------------------------------------------- */
483: /*
484:    PCISApplySchur -

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

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

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

503:   if (!vec2_B) vec2_B = v;

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

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

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

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

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

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

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

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

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

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

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

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

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