Actual source code: pcis.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  2: #include <../src/ksp/pc/impls/is/pcis.h> /*I "petscpc.h" I*/

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

 11:   pcis->use_stiffness_scaling = use;
 12:   return(0);
 13: }

 17: /*@
 18:  PCISSetUseStiffnessScaling - Tells PCIS to construct partition of unity using
 19:                               local matrices' diagonal.

 21:    Not collective

 23:    Input Parameters:
 24: +  pc - the preconditioning context
 25: -  use - whether or not pcis use matrix diagonal to build partition of unity.

 27:    Level: intermediate

 29:    Notes:

 31: .seealso: PCBDDC
 32: @*/
 33: PetscErrorCode PCISSetUseStiffnessScaling(PC pc, PetscBool use)
 34: {

 39:   PetscTryMethod(pc,"PCISSetUseStiffnessScaling_C",(PC,PetscBool),(pc,use));
 40:   return(0);
 41: }

 45: static PetscErrorCode PCISSetSubdomainDiagonalScaling_IS(PC pc, Vec scaling_factors)
 46: {
 48:   PC_IS          *pcis = (PC_IS*)pc->data;

 51:   VecDestroy(&pcis->D);
 52:   PetscObjectReference((PetscObject)scaling_factors);
 53:   pcis->D = scaling_factors;
 54:   return(0);
 55: }

 59: /*@
 60:  PCISSetSubdomainDiagonalScaling - Set diagonal scaling for PCIS.

 62:    Not collective

 64:    Input Parameters:
 65: +  pc - the preconditioning context
 66: -  scaling_factors - scaling factors for the subdomain

 68:    Level: intermediate

 70:    Notes:
 71:    Intended to use with jumping coefficients cases.

 73: .seealso: PCBDDC
 74: @*/
 75: PetscErrorCode PCISSetSubdomainDiagonalScaling(PC pc, Vec scaling_factors)
 76: {

 81:   PetscTryMethod(pc,"PCISSetSubdomainDiagonalScaling_C",(PC,Vec),(pc,scaling_factors));
 82:   return(0);
 83: }

 87: static PetscErrorCode PCISSetSubdomainScalingFactor_IS(PC pc, PetscScalar scal)
 88: {
 89:   PC_IS *pcis = (PC_IS*)pc->data;

 92:   pcis->scaling_factor = scal;
 93:   return(0);
 94: }

 98: /*@
 99:  PCISSetSubdomainScalingFactor - Set scaling factor for PCIS.

101:    Not collective

103:    Input Parameters:
104: +  pc - the preconditioning context
105: -  scal - scaling factor for the subdomain

107:    Level: intermediate

109:    Notes:
110:    Intended to use with jumping coefficients cases.

112: .seealso: PCBDDC
113: @*/
114: PetscErrorCode PCISSetSubdomainScalingFactor(PC pc, PetscScalar scal)
115: {

120:   PetscTryMethod(pc,"PCISSetSubdomainScalingFactor_C",(PC,PetscScalar),(pc,scal));
121:   return(0);
122: }


125: /* -------------------------------------------------------------------------- */
126: /*
127:    PCISSetUp -
128: */
131: PetscErrorCode  PCISSetUp(PC pc, PetscBool computesolvers)
132: {
133:   PC_IS          *pcis  = (PC_IS*)(pc->data);
134:   Mat_IS         *matis;
135:   MatReuse       reuse;
137:   PetscBool      flg,issbaij;
138:   Vec            counter;

141:   PetscObjectTypeCompare((PetscObject)pc->pmat,MATIS,&flg);
142:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Preconditioner type of Neumann Neumman requires matrix of type MATIS");
143:   matis = (Mat_IS*)pc->pmat->data;

145:   /* first time creation, get info on substructuring */
146:   if (!pc->setupcalled) {
147:     PetscInt    n_I;
148:     PetscInt    *idx_I_local,*idx_B_local,*idx_I_global,*idx_B_global;
149:     PetscInt    *array;
150:     PetscInt    i,j;

152:     /* get info on mapping */
153:     PetscObjectReference((PetscObject)pc->pmat->rmap->mapping);
154:     ISLocalToGlobalMappingDestroy(&pcis->mapping);
155:     pcis->mapping = pc->pmat->rmap->mapping;
156:     ISLocalToGlobalMappingGetSize(pcis->mapping,&pcis->n);
157:     ISLocalToGlobalMappingGetInfo(pcis->mapping,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));

159:     /* Identifying interior and interface nodes, in local numbering */
160:     PetscMalloc1(pcis->n,&array);
161:     PetscMemzero(array,pcis->n*sizeof(PetscInt));
162:     for (i=0;i<pcis->n_neigh;i++)
163:       for (j=0;j<pcis->n_shared[i];j++)
164:           array[pcis->shared[i][j]] += 1;

166:     /* Creating local and global index sets for interior and inteface nodes. */
167:     PetscMalloc1(pcis->n,&idx_I_local);
168:     PetscMalloc1(pcis->n,&idx_B_local);
169:     for (i=0, pcis->n_B=0, n_I=0; i<pcis->n; i++) {
170:       if (!array[i]) {
171:         idx_I_local[n_I] = i;
172:         n_I++;
173:       } else {
174:         idx_B_local[pcis->n_B] = i;
175:         pcis->n_B++;
176:       }
177:     }
178:     /* Getting the global numbering */
179:     idx_B_global = idx_I_local + n_I; /* Just avoiding allocating extra memory, since we have vacant space */
180:     idx_I_global = idx_B_local + pcis->n_B;
181:     ISLocalToGlobalMappingApply(pcis->mapping,pcis->n_B,idx_B_local,idx_B_global);
182:     ISLocalToGlobalMappingApply(pcis->mapping,n_I,idx_I_local,idx_I_global);

184:     /* Creating the index sets */
185:     ISCreateGeneral(PETSC_COMM_SELF,pcis->n_B,idx_B_local,PETSC_COPY_VALUES, &pcis->is_B_local);
186:     ISCreateGeneral(PETSC_COMM_SELF,pcis->n_B,idx_B_global,PETSC_COPY_VALUES,&pcis->is_B_global);
187:     ISCreateGeneral(PETSC_COMM_SELF,n_I,idx_I_local,PETSC_COPY_VALUES, &pcis->is_I_local);
188:     ISCreateGeneral(PETSC_COMM_SELF,n_I,idx_I_global,PETSC_COPY_VALUES,&pcis->is_I_global);

190:     /* Freeing memory */
191:     PetscFree(idx_B_local);
192:     PetscFree(idx_I_local);
193:     PetscFree(array);

195:     /* Creating work vectors and arrays */
196:     VecDuplicate(matis->x,&pcis->vec1_N);
197:     VecDuplicate(pcis->vec1_N,&pcis->vec2_N);
198:     VecCreateSeq(PETSC_COMM_SELF,pcis->n-pcis->n_B,&pcis->vec1_D);
199:     VecDuplicate(pcis->vec1_D,&pcis->vec2_D);
200:     VecDuplicate(pcis->vec1_D,&pcis->vec3_D);
201:     VecDuplicate(pcis->vec1_D,&pcis->vec4_D);
202:     VecCreateSeq(PETSC_COMM_SELF,pcis->n_B,&pcis->vec1_B);
203:     VecDuplicate(pcis->vec1_B,&pcis->vec2_B);
204:     VecDuplicate(pcis->vec1_B,&pcis->vec3_B);
205:     MatCreateVecs(pc->pmat,&pcis->vec1_global,0);
206:     PetscMalloc1(pcis->n,&pcis->work_N);
207:     /* scaling vector */
208:     VecDuplicate(pcis->vec1_B,&pcis->D);

210:     /* Creating the scatter contexts */
211:     VecScatterCreate(pcis->vec1_global,pcis->is_I_global,pcis->vec1_D,(IS)0,&pcis->global_to_D);
212:     VecScatterCreate(pcis->vec1_N,pcis->is_B_local,pcis->vec1_B,(IS)0,&pcis->N_to_B);
213:     VecScatterCreate(pcis->vec1_global,pcis->is_B_global,pcis->vec1_B,(IS)0,&pcis->global_to_B);

215:     /* map from boundary to local */
216:     ISLocalToGlobalMappingCreateIS(pcis->is_B_local,&pcis->BtoNmap);
217:   }

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

223:         [ A_II | A_IB ]
224:     A = [------+------]
225:         [ A_BI | A_BB ]
226:   */
227:   reuse = MAT_INITIAL_MATRIX;
228:   if (pcis->reusesubmatrices && pc->setupcalled) {
229:     if (pc->flag == SAME_NONZERO_PATTERN) {
230:       reuse = MAT_REUSE_MATRIX;
231:     } else {
232:       reuse = MAT_INITIAL_MATRIX;
233:     }
234:   }
235:   if (reuse == MAT_INITIAL_MATRIX) {
236:     MatDestroy(&pcis->A_II);
237:     MatDestroy(&pcis->A_IB);
238:     MatDestroy(&pcis->A_BI);
239:     MatDestroy(&pcis->A_BB);
240:   }

242:   MatGetSubMatrix(matis->A,pcis->is_I_local,pcis->is_I_local,reuse,&pcis->A_II);
243:   MatGetSubMatrix(matis->A,pcis->is_B_local,pcis->is_B_local,reuse,&pcis->A_BB);
244:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);
245:   if (!issbaij) {
246:     MatGetSubMatrix(matis->A,pcis->is_I_local,pcis->is_B_local,reuse,&pcis->A_IB);
247:     MatGetSubMatrix(matis->A,pcis->is_B_local,pcis->is_I_local,reuse,&pcis->A_BI);
248:   } else {
249:     Mat newmat;
250:     MatConvert(matis->A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&newmat);
251:     MatGetSubMatrix(newmat,pcis->is_I_local,pcis->is_B_local,reuse,&pcis->A_IB);
252:     MatGetSubMatrix(newmat,pcis->is_B_local,pcis->is_I_local,reuse,&pcis->A_BI);
253:     MatDestroy(&newmat);
254:   }

256:   /* Creating scaling "matrix" D */
257:   PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_is_use_stiffness_scaling",&pcis->use_stiffness_scaling,NULL);
258:   if (!pcis->use_stiffness_scaling) {
259:     VecSet(pcis->D,pcis->scaling_factor);
260:   } else {
261:     MatGetDiagonal(matis->A,pcis->vec1_N);
262:     VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);
263:     VecScatterEnd(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);
264:   }
265:   VecCopy(pcis->D,pcis->vec1_B);
266:   MatCreateVecs(pc->pmat,&counter,0); /* temporary auxiliar vector */
267:   VecSet(counter,0.0);
268:   VecScatterBegin(pcis->global_to_B,pcis->vec1_B,counter,ADD_VALUES,SCATTER_REVERSE);
269:   VecScatterEnd(pcis->global_to_B,pcis->vec1_B,counter,ADD_VALUES,SCATTER_REVERSE);
270:   VecScatterBegin(pcis->global_to_B,counter,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
271:   VecScatterEnd(pcis->global_to_B,counter,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
272:   VecPointwiseDivide(pcis->D,pcis->D,pcis->vec1_B);
273:   VecDestroy(&counter);

275:   /* See historical note 01, at the bottom of this file. */

277:   /* Creating the KSP contexts for the local Dirichlet and Neumann problems */
278:   if (computesolvers) {
279:     PC pc_ctx;

281:     pcis->pure_neumann = matis->pure_neumann;
282:     /* Dirichlet */
283:     KSPCreate(PETSC_COMM_SELF,&pcis->ksp_D);
284:     KSPSetErrorIfNotConverged(pcis->ksp_D,pc->erroriffailure);
285:     PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_D,(PetscObject)pc,1);
286:     KSPSetOperators(pcis->ksp_D,pcis->A_II,pcis->A_II);
287:     KSPSetOptionsPrefix(pcis->ksp_D,"is_localD_");
288:     KSPGetPC(pcis->ksp_D,&pc_ctx);
289:     PCSetType(pc_ctx,PCLU);
290:     KSPSetType(pcis->ksp_D,KSPPREONLY);
291:     KSPSetFromOptions(pcis->ksp_D);
292:     /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
293:     KSPSetUp(pcis->ksp_D);
294:     /* Neumann */
295:     KSPCreate(PETSC_COMM_SELF,&pcis->ksp_N);
296:     KSPSetErrorIfNotConverged(pcis->ksp_N,pc->erroriffailure);
297:     PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_N,(PetscObject)pc,1);
298:     KSPSetOperators(pcis->ksp_N,matis->A,matis->A);
299:     KSPSetOptionsPrefix(pcis->ksp_N,"is_localN_");
300:     KSPGetPC(pcis->ksp_N,&pc_ctx);
301:     PCSetType(pc_ctx,PCLU);
302:     KSPSetType(pcis->ksp_N,KSPPREONLY);
303:     KSPSetFromOptions(pcis->ksp_N);
304:     {
305:       PetscBool damp_fixed                    = PETSC_FALSE,
306:                 remove_nullspace_fixed        = PETSC_FALSE,
307:                 set_damping_factor_floating   = PETSC_FALSE,
308:                 not_damp_floating             = PETSC_FALSE,
309:                 not_remove_nullspace_floating = PETSC_FALSE;
310:       PetscReal fixed_factor,
311:                 floating_factor;

313:       PetscOptionsGetReal(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&fixed_factor,&damp_fixed);
314:       if (!damp_fixed) fixed_factor = 0.0;
315:       PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&damp_fixed,NULL);

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

319:       PetscOptionsGetReal(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_set_damping_factor_floating",
320:                               &floating_factor,&set_damping_factor_floating);
321:       if (!set_damping_factor_floating) floating_factor = 0.0;
322:       PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_set_damping_factor_floating",&set_damping_factor_floating,NULL);
323:       if (!set_damping_factor_floating) floating_factor = 1.e-12;

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

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

329:       if (pcis->pure_neumann) {  /* floating subdomain */
330:         if (!(not_damp_floating)) {
331:           PCFactorSetShiftType(pc_ctx,MAT_SHIFT_NONZERO);
332:           PCFactorSetShiftAmount(pc_ctx,floating_factor);
333:         }
334:         if (!(not_remove_nullspace_floating)) {
335:           MatNullSpace nullsp;
336:           MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_TRUE,0,NULL,&nullsp);
337:           MatSetNullSpace(matis->A,nullsp);
338:           MatNullSpaceDestroy(&nullsp);
339:         }
340:       } else {  /* fixed subdomain */
341:         if (damp_fixed) {
342:           PCFactorSetShiftType(pc_ctx,MAT_SHIFT_NONZERO);
343:           PCFactorSetShiftAmount(pc_ctx,floating_factor);
344:         }
345:         if (remove_nullspace_fixed) {
346:           MatNullSpace nullsp;
347:           MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_TRUE,0,NULL,&nullsp);
348:           MatSetNullSpace(matis->A,nullsp);
349:           MatNullSpaceDestroy(&nullsp);
350:         }
351:       }
352:     }
353:     /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
354:     KSPSetUp(pcis->ksp_N);
355:   }

357:   return(0);
358: }

360: /* -------------------------------------------------------------------------- */
361: /*
362:    PCISDestroy -
363: */
366: PetscErrorCode  PCISDestroy(PC pc)
367: {
368:   PC_IS          *pcis = (PC_IS*)(pc->data);

372:   ISDestroy(&pcis->is_B_local);
373:   ISDestroy(&pcis->is_I_local);
374:   ISDestroy(&pcis->is_B_global);
375:   ISDestroy(&pcis->is_I_global);
376:   MatDestroy(&pcis->A_II);
377:   MatDestroy(&pcis->A_IB);
378:   MatDestroy(&pcis->A_BI);
379:   MatDestroy(&pcis->A_BB);
380:   VecDestroy(&pcis->D);
381:   KSPDestroy(&pcis->ksp_N);
382:   KSPDestroy(&pcis->ksp_D);
383:   VecDestroy(&pcis->vec1_N);
384:   VecDestroy(&pcis->vec2_N);
385:   VecDestroy(&pcis->vec1_D);
386:   VecDestroy(&pcis->vec2_D);
387:   VecDestroy(&pcis->vec3_D);
388:   VecDestroy(&pcis->vec4_D);
389:   VecDestroy(&pcis->vec1_B);
390:   VecDestroy(&pcis->vec2_B);
391:   VecDestroy(&pcis->vec3_B);
392:   VecDestroy(&pcis->vec1_global);
393:   VecScatterDestroy(&pcis->global_to_D);
394:   VecScatterDestroy(&pcis->N_to_B);
395:   VecScatterDestroy(&pcis->global_to_B);
396:   PetscFree(pcis->work_N);
397:   if (pcis->n_neigh > -1) {
398:     ISLocalToGlobalMappingRestoreInfo(pcis->mapping,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));
399:   }
400:   ISLocalToGlobalMappingDestroy(&pcis->mapping);
401:   ISLocalToGlobalMappingDestroy(&pcis->BtoNmap);
402:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetUseStiffnessScaling_C",NULL);
403:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainScalingFactor_C",NULL);
404:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainDiagonalScaling_C",NULL);
405:   return(0);
406: }

408: /* -------------------------------------------------------------------------- */
409: /*
410:    PCISCreate -
411: */
414: PetscErrorCode  PCISCreate(PC pc)
415: {
416:   PC_IS          *pcis = (PC_IS*)(pc->data);

420:   pcis->is_B_local       = 0;
421:   pcis->is_I_local       = 0;
422:   pcis->is_B_global      = 0;
423:   pcis->is_I_global      = 0;
424:   pcis->A_II             = 0;
425:   pcis->A_IB             = 0;
426:   pcis->A_BI             = 0;
427:   pcis->A_BB             = 0;
428:   pcis->D                = 0;
429:   pcis->ksp_N            = 0;
430:   pcis->ksp_D            = 0;
431:   pcis->vec1_N           = 0;
432:   pcis->vec2_N           = 0;
433:   pcis->vec1_D           = 0;
434:   pcis->vec2_D           = 0;
435:   pcis->vec3_D           = 0;
436:   pcis->vec1_B           = 0;
437:   pcis->vec2_B           = 0;
438:   pcis->vec3_B           = 0;
439:   pcis->vec1_global      = 0;
440:   pcis->work_N           = 0;
441:   pcis->global_to_D      = 0;
442:   pcis->N_to_B           = 0;
443:   pcis->global_to_B      = 0;
444:   pcis->mapping          = 0;
445:   pcis->BtoNmap          = 0;
446:   pcis->n_neigh          = -1;
447:   pcis->scaling_factor   = 1.0;
448:   pcis->reusesubmatrices = PETSC_TRUE;
449:   /* composing functions */
450:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetUseStiffnessScaling_C",PCISSetUseStiffnessScaling_IS);
451:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainScalingFactor_C",PCISSetSubdomainScalingFactor_IS);
452:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainDiagonalScaling_C",PCISSetSubdomainDiagonalScaling_IS);
453:   return(0);
454: }

456: /* -------------------------------------------------------------------------- */
457: /*
458:    PCISApplySchur -

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

464:    Output parameters:
465: .  vec1_B - result of Schur complement applied to chunk
466: .  vec2_B - garbage (used as work space), or null (and v is used as workspace)
467: .  vec1_D - garbage (used as work space)
468: .  vec2_D - garbage (used as work space)

470: */
473: PetscErrorCode  PCISApplySchur(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D)
474: {
476:   PC_IS          *pcis = (PC_IS*)(pc->data);

479:   if (!vec2_B) vec2_B = v;

481:   MatMult(pcis->A_BB,v,vec1_B);
482:   MatMult(pcis->A_IB,v,vec1_D);
483:   KSPSolve(pcis->ksp_D,vec1_D,vec2_D);
484:   MatMult(pcis->A_BI,vec2_D,vec2_B);
485:   VecAXPY(vec1_B,-1.0,vec2_B);
486:   return(0);
487: }

489: /* -------------------------------------------------------------------------- */
490: /*
491:    PCISScatterArrayNToVecB - Scatters interface node values from a big array (of all local nodes, interior or interface,
492:    including ghosts) into an interface vector, when in SCATTER_FORWARD mode, or vice-versa, when in SCATTER_REVERSE
493:    mode.

495:    Input parameters:
496: .  pc - preconditioner context
497: .  array_N - [when in SCATTER_FORWARD mode] Array to be scattered into the vector
498: .  v_B - [when in SCATTER_REVERSE mode] Vector to be scattered into the array

500:    Output parameter:
501: .  array_N - [when in SCATTER_REVERSE mode] Array to receive the scattered vector
502: .  v_B - [when in SCATTER_FORWARD mode] Vector to receive the scattered array

504:    Notes:
505:    The entries in the array that do not correspond to interface nodes remain unaltered.
506: */
509: PetscErrorCode  PCISScatterArrayNToVecB(PetscScalar *array_N, Vec v_B, InsertMode imode, ScatterMode smode, PC pc)
510: {
511:   PetscInt       i;
512:   const PetscInt *idex;
514:   PetscScalar    *array_B;
515:   PC_IS          *pcis = (PC_IS*)(pc->data);

518:   VecGetArray(v_B,&array_B);
519:   ISGetIndices(pcis->is_B_local,&idex);

521:   if (smode == SCATTER_FORWARD) {
522:     if (imode == INSERT_VALUES) {
523:       for (i=0; i<pcis->n_B; i++) array_B[i] = array_N[idex[i]];
524:     } else {  /* ADD_VALUES */
525:       for (i=0; i<pcis->n_B; i++) array_B[i] += array_N[idex[i]];
526:     }
527:   } else {  /* SCATTER_REVERSE */
528:     if (imode == INSERT_VALUES) {
529:       for (i=0; i<pcis->n_B; i++) array_N[idex[i]] = array_B[i];
530:     } else {  /* ADD_VALUES */
531:       for (i=0; i<pcis->n_B; i++) array_N[idex[i]] += array_B[i];
532:     }
533:   }
534:   ISRestoreIndices(pcis->is_B_local,&idex);
535:   VecRestoreArray(v_B,&array_B);
536:   return(0);
537: }

539: /* -------------------------------------------------------------------------- */
540: /*
541:    PCISApplyInvSchur - Solves the Neumann problem related to applying the inverse of the Schur complement.
542:    More precisely, solves the problem:
543:                                         [ A_II  A_IB ] [ . ]   [ 0 ]
544:                                         [            ] [   ] = [   ]
545:                                         [ A_BI  A_BB ] [ x ]   [ b ]

547:    Input parameters:
548: .  pc - preconditioner context
549: .  b - vector of local interface nodes (including ghosts)

551:    Output parameters:
552: .  x - vector of local interface nodes (including ghosts); returns the application of the inverse of the Schur
553:        complement to b
554: .  vec1_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space)
555: .  vec2_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space)

557: */
560: PetscErrorCode  PCISApplyInvSchur(PC pc, Vec b, Vec x, Vec vec1_N, Vec vec2_N)
561: {
563:   PC_IS          *pcis = (PC_IS*)(pc->data);

566:   /*
567:     Neumann solvers.
568:     Applying the inverse of the local Schur complement, i.e, solving a Neumann
569:     Problem with zero at the interior nodes of the RHS and extracting the interface
570:     part of the solution. inverse Schur complement is applied to b and the result
571:     is stored in x.
572:   */
573:   /* Setting the RHS vec1_N */
574:   VecSet(vec1_N,0.0);
575:   VecScatterBegin(pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);
576:   VecScatterEnd  (pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);
577:   /* Checking for consistency of the RHS */
578:   {
579:     PetscBool flg = PETSC_FALSE;
580:     PetscOptionsGetBool(NULL,NULL,"-pc_is_check_consistency",&flg,NULL);
581:     if (flg) {
582:       PetscScalar average;
583:       PetscViewer viewer;
584:       PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);

586:       VecSum(vec1_N,&average);
587:       average = average / ((PetscReal)pcis->n);
588:       PetscViewerASCIIPushSynchronized(viewer);
589:       if (pcis->pure_neumann) {
590:         PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is floating. Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));
591:       } else {
592:         PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is fixed.    Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));
593:       }
594:       PetscViewerFlush(viewer);
595:       PetscViewerASCIIPopSynchronized(viewer);
596:     }
597:   }
598:   /* Solving the system for vec2_N */
599:   KSPSolve(pcis->ksp_N,vec1_N,vec2_N);
600:   /* Extracting the local interface vector out of the solution */
601:   VecScatterBegin(pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);
602:   VecScatterEnd  (pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);
603:   return(0);
604: }