Actual source code: pcis.c

petsc-3.8.4 2018-03-24
Report Typos and Errors

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

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

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

 45:   PetscObjectReference((PetscObject)scaling_factors);
 46:   VecDestroy(&pcis->D);
 47:   pcis->D = scaling_factors;
 48:   return(0);
 49: }

 51: /*@
 52:  PCISSetSubdomainDiagonalScaling - Set diagonal scaling for PCIS.

 54:    Not collective

 56:    Input Parameters:
 57: +  pc - the preconditioning context
 58: -  scaling_factors - scaling factors for the subdomain

 60:    Level: intermediate

 62:    Notes:
 63:    Intended to use with jumping coefficients cases.

 65: .seealso: PCBDDC
 66: @*/
 67: PetscErrorCode PCISSetSubdomainDiagonalScaling(PC pc, Vec scaling_factors)
 68: {

 73:   PetscTryMethod(pc,"PCISSetSubdomainDiagonalScaling_C",(PC,Vec),(pc,scaling_factors));
 74:   return(0);
 75: }

 77: static PetscErrorCode PCISSetSubdomainScalingFactor_IS(PC pc, PetscScalar scal)
 78: {
 79:   PC_IS *pcis = (PC_IS*)pc->data;

 82:   pcis->scaling_factor = scal;
 83:   return(0);
 84: }

 86: /*@
 87:  PCISSetSubdomainScalingFactor - Set scaling factor for PCIS.

 89:    Not collective

 91:    Input Parameters:
 92: +  pc - the preconditioning context
 93: -  scal - scaling factor for the subdomain

 95:    Level: intermediate

 97:    Notes:
 98:    Intended to use with jumping coefficients cases.

100: .seealso: PCBDDC
101: @*/
102: PetscErrorCode PCISSetSubdomainScalingFactor(PC pc, PetscScalar scal)
103: {

108:   PetscTryMethod(pc,"PCISSetSubdomainScalingFactor_C",(PC,PetscScalar),(pc,scal));
109:   return(0);
110: }


113: /* -------------------------------------------------------------------------- */
114: /*
115:    PCISSetUp -
116: */
117: PetscErrorCode  PCISSetUp(PC pc, PetscBool computesolvers)
118: {
119:   PC_IS          *pcis  = (PC_IS*)(pc->data);
120:   Mat_IS         *matis;
121:   MatReuse       reuse;
123:   PetscBool      flg,issbaij;
124:   Vec            counter;

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

131:   /* first time creation, get info on substructuring */
132:   if (!pc->setupcalled) {
133:     PetscInt    n_I;
134:     PetscInt    *idx_I_local,*idx_B_local,*idx_I_global,*idx_B_global;
135:     PetscBT     bt;
136:     PetscInt    i,j;

138:     /* get info on mapping */
139:     PetscObjectReference((PetscObject)pc->pmat->rmap->mapping);
140:     ISLocalToGlobalMappingDestroy(&pcis->mapping);
141:     pcis->mapping = pc->pmat->rmap->mapping;
142:     ISLocalToGlobalMappingGetSize(pcis->mapping,&pcis->n);
143:     ISLocalToGlobalMappingGetInfo(pcis->mapping,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));

145:     /* Identifying interior and interface nodes, in local numbering */
146:     PetscBTCreate(pcis->n,&bt);
147:     for (i=0;i<pcis->n_neigh;i++)
148:       for (j=0;j<pcis->n_shared[i];j++) {
149:           PetscBTSet(bt,pcis->shared[i][j]);
150:       }

152:     /* Creating local and global index sets for interior and inteface nodes. */
153:     PetscMalloc1(pcis->n,&idx_I_local);
154:     PetscMalloc1(pcis->n,&idx_B_local);
155:     for (i=0, pcis->n_B=0, n_I=0; i<pcis->n; i++) {
156:       if (!PetscBTLookup(bt,i)) {
157:         idx_I_local[n_I] = i;
158:         n_I++;
159:       } else {
160:         idx_B_local[pcis->n_B] = i;
161:         pcis->n_B++;
162:       }
163:     }

165:     /* Getting the global numbering */
166:     idx_B_global = idx_I_local + n_I; /* Just avoiding allocating extra memory, since we have vacant space */
167:     idx_I_global = idx_B_local + pcis->n_B;
168:     ISLocalToGlobalMappingApply(pcis->mapping,pcis->n_B,idx_B_local,idx_B_global);
169:     ISLocalToGlobalMappingApply(pcis->mapping,n_I,idx_I_local,idx_I_global);

171:     /* Creating the index sets */
172:     ISCreateGeneral(PETSC_COMM_SELF,pcis->n_B,idx_B_local,PETSC_COPY_VALUES, &pcis->is_B_local);
173:     ISCreateGeneral(PETSC_COMM_SELF,pcis->n_B,idx_B_global,PETSC_COPY_VALUES,&pcis->is_B_global);
174:     ISCreateGeneral(PETSC_COMM_SELF,n_I,idx_I_local,PETSC_COPY_VALUES, &pcis->is_I_local);
175:     ISCreateGeneral(PETSC_COMM_SELF,n_I,idx_I_global,PETSC_COPY_VALUES,&pcis->is_I_global);

177:     /* Freeing memory */
178:     PetscFree(idx_B_local);
179:     PetscFree(idx_I_local);
180:     PetscBTDestroy(&bt);

182:     /* Creating work vectors and arrays */
183:     VecDuplicate(matis->x,&pcis->vec1_N);
184:     VecDuplicate(pcis->vec1_N,&pcis->vec2_N);
185:     VecCreateSeq(PETSC_COMM_SELF,pcis->n-pcis->n_B,&pcis->vec1_D);
186:     VecDuplicate(pcis->vec1_D,&pcis->vec2_D);
187:     VecDuplicate(pcis->vec1_D,&pcis->vec3_D);
188:     VecDuplicate(pcis->vec1_D,&pcis->vec4_D);
189:     VecCreateSeq(PETSC_COMM_SELF,pcis->n_B,&pcis->vec1_B);
190:     VecDuplicate(pcis->vec1_B,&pcis->vec2_B);
191:     VecDuplicate(pcis->vec1_B,&pcis->vec3_B);
192:     MatCreateVecs(pc->pmat,&pcis->vec1_global,0);
193:     PetscMalloc1(pcis->n,&pcis->work_N);
194:     /* scaling vector */
195:     if (!pcis->D) { /* it can happen that the user passed in a scaling vector via PCISSetSubdomainDiagonalScaling */
196:       VecDuplicate(pcis->vec1_B,&pcis->D);
197:       VecSet(pcis->D,pcis->scaling_factor);
198:     }

200:     /* Creating the scatter contexts */
201:     VecScatterCreate(pcis->vec1_N,pcis->is_I_local,pcis->vec1_D,(IS)0,&pcis->N_to_D);
202:     VecScatterCreate(pcis->vec1_global,pcis->is_I_global,pcis->vec1_D,(IS)0,&pcis->global_to_D);
203:     VecScatterCreate(pcis->vec1_N,pcis->is_B_local,pcis->vec1_B,(IS)0,&pcis->N_to_B);
204:     VecScatterCreate(pcis->vec1_global,pcis->is_B_global,pcis->vec1_B,(IS)0,&pcis->global_to_B);

206:     /* map from boundary to local */
207:     ISLocalToGlobalMappingCreateIS(pcis->is_B_local,&pcis->BtoNmap);
208:   }

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

214:         [ A_II | A_IB ]
215:     A = [------+------]
216:         [ A_BI | A_BB ]
217:   */
218:   reuse = MAT_INITIAL_MATRIX;
219:   if (pcis->reusesubmatrices && pc->setupcalled) {
220:     if (pc->flag == SAME_NONZERO_PATTERN) {
221:       reuse = MAT_REUSE_MATRIX;
222:     } else {
223:       reuse = MAT_INITIAL_MATRIX;
224:     }
225:   }
226:   if (reuse == MAT_INITIAL_MATRIX) {
227:     MatDestroy(&pcis->A_II);
228:     MatDestroy(&pcis->A_IB);
229:     MatDestroy(&pcis->A_BI);
230:     MatDestroy(&pcis->A_BB);
231:   }

233:   MatCreateSubMatrix(matis->A,pcis->is_I_local,pcis->is_I_local,reuse,&pcis->A_II);
234:   MatCreateSubMatrix(matis->A,pcis->is_B_local,pcis->is_B_local,reuse,&pcis->A_BB);
235:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);
236:   if (!issbaij) {
237:     MatCreateSubMatrix(matis->A,pcis->is_I_local,pcis->is_B_local,reuse,&pcis->A_IB);
238:     MatCreateSubMatrix(matis->A,pcis->is_B_local,pcis->is_I_local,reuse,&pcis->A_BI);
239:   } else {
240:     Mat newmat;
241:     MatConvert(matis->A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&newmat);
242:     MatCreateSubMatrix(newmat,pcis->is_I_local,pcis->is_B_local,reuse,&pcis->A_IB);
243:     MatCreateSubMatrix(newmat,pcis->is_B_local,pcis->is_I_local,reuse,&pcis->A_BI);
244:     MatDestroy(&newmat);
245:   }

247:   /* Creating scaling vector D */
248:   PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_is_use_stiffness_scaling",&pcis->use_stiffness_scaling,NULL);
249:   if (pcis->use_stiffness_scaling) {
250:     MatGetDiagonal(pcis->A_BB,pcis->D);
251:   }
252:   MatCreateVecs(pc->pmat,&counter,0); /* temporary auxiliar vector */
253:   VecSet(counter,0.0);
254:   VecScatterBegin(pcis->global_to_B,pcis->D,counter,ADD_VALUES,SCATTER_REVERSE);
255:   VecScatterEnd(pcis->global_to_B,pcis->D,counter,ADD_VALUES,SCATTER_REVERSE);
256:   VecScatterBegin(pcis->global_to_B,counter,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
257:   VecScatterEnd(pcis->global_to_B,counter,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
258:   VecPointwiseDivide(pcis->D,pcis->D,pcis->vec1_B);
259:   VecDestroy(&counter);

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

263:   /* Creating the KSP contexts for the local Dirichlet and Neumann problems */
264:   if (computesolvers) {
265:     PC pc_ctx;

267:     pcis->pure_neumann = matis->pure_neumann;
268:     /* Dirichlet */
269:     KSPCreate(PETSC_COMM_SELF,&pcis->ksp_D);
270:     KSPSetErrorIfNotConverged(pcis->ksp_D,pc->erroriffailure);
271:     PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_D,(PetscObject)pc,1);
272:     KSPSetOperators(pcis->ksp_D,pcis->A_II,pcis->A_II);
273:     KSPSetOptionsPrefix(pcis->ksp_D,"is_localD_");
274:     KSPGetPC(pcis->ksp_D,&pc_ctx);
275:     PCSetType(pc_ctx,PCLU);
276:     KSPSetType(pcis->ksp_D,KSPPREONLY);
277:     KSPSetFromOptions(pcis->ksp_D);
278:     /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
279:     KSPSetUp(pcis->ksp_D);
280:     /* Neumann */
281:     KSPCreate(PETSC_COMM_SELF,&pcis->ksp_N);
282:     KSPSetErrorIfNotConverged(pcis->ksp_N,pc->erroriffailure);
283:     PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_N,(PetscObject)pc,1);
284:     KSPSetOperators(pcis->ksp_N,matis->A,matis->A);
285:     KSPSetOptionsPrefix(pcis->ksp_N,"is_localN_");
286:     KSPGetPC(pcis->ksp_N,&pc_ctx);
287:     PCSetType(pc_ctx,PCLU);
288:     KSPSetType(pcis->ksp_N,KSPPREONLY);
289:     KSPSetFromOptions(pcis->ksp_N);
290:     {
291:       PetscBool damp_fixed                    = PETSC_FALSE,
292:                 remove_nullspace_fixed        = PETSC_FALSE,
293:                 set_damping_factor_floating   = PETSC_FALSE,
294:                 not_damp_floating             = PETSC_FALSE,
295:                 not_remove_nullspace_floating = PETSC_FALSE;
296:       PetscReal fixed_factor,
297:                 floating_factor;

299:       PetscOptionsGetReal(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&fixed_factor,&damp_fixed);
300:       if (!damp_fixed) fixed_factor = 0.0;
301:       PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&damp_fixed,NULL);

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

305:       PetscOptionsGetReal(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_set_damping_factor_floating",
306:                               &floating_factor,&set_damping_factor_floating);
307:       if (!set_damping_factor_floating) floating_factor = 0.0;
308:       PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_set_damping_factor_floating",&set_damping_factor_floating,NULL);
309:       if (!set_damping_factor_floating) floating_factor = 1.e-12;

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

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

315:       if (pcis->pure_neumann) {  /* floating subdomain */
316:         if (!(not_damp_floating)) {
317:           PCFactorSetShiftType(pc_ctx,MAT_SHIFT_NONZERO);
318:           PCFactorSetShiftAmount(pc_ctx,floating_factor);
319:         }
320:         if (!(not_remove_nullspace_floating)) {
321:           MatNullSpace nullsp;
322:           MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_TRUE,0,NULL,&nullsp);
323:           MatSetNullSpace(matis->A,nullsp);
324:           MatNullSpaceDestroy(&nullsp);
325:         }
326:       } else {  /* fixed subdomain */
327:         if (damp_fixed) {
328:           PCFactorSetShiftType(pc_ctx,MAT_SHIFT_NONZERO);
329:           PCFactorSetShiftAmount(pc_ctx,floating_factor);
330:         }
331:         if (remove_nullspace_fixed) {
332:           MatNullSpace nullsp;
333:           MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_TRUE,0,NULL,&nullsp);
334:           MatSetNullSpace(matis->A,nullsp);
335:           MatNullSpaceDestroy(&nullsp);
336:         }
337:       }
338:     }
339:     /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
340:     KSPSetUp(pcis->ksp_N);
341:   }

343:   return(0);
344: }

346: /* -------------------------------------------------------------------------- */
347: /*
348:    PCISDestroy -
349: */
350: PetscErrorCode  PCISDestroy(PC pc)
351: {
352:   PC_IS          *pcis = (PC_IS*)(pc->data);

356:   ISDestroy(&pcis->is_B_local);
357:   ISDestroy(&pcis->is_I_local);
358:   ISDestroy(&pcis->is_B_global);
359:   ISDestroy(&pcis->is_I_global);
360:   MatDestroy(&pcis->A_II);
361:   MatDestroy(&pcis->A_IB);
362:   MatDestroy(&pcis->A_BI);
363:   MatDestroy(&pcis->A_BB);
364:   VecDestroy(&pcis->D);
365:   KSPDestroy(&pcis->ksp_N);
366:   KSPDestroy(&pcis->ksp_D);
367:   VecDestroy(&pcis->vec1_N);
368:   VecDestroy(&pcis->vec2_N);
369:   VecDestroy(&pcis->vec1_D);
370:   VecDestroy(&pcis->vec2_D);
371:   VecDestroy(&pcis->vec3_D);
372:   VecDestroy(&pcis->vec4_D);
373:   VecDestroy(&pcis->vec1_B);
374:   VecDestroy(&pcis->vec2_B);
375:   VecDestroy(&pcis->vec3_B);
376:   VecDestroy(&pcis->vec1_global);
377:   VecScatterDestroy(&pcis->global_to_D);
378:   VecScatterDestroy(&pcis->N_to_B);
379:   VecScatterDestroy(&pcis->N_to_D);
380:   VecScatterDestroy(&pcis->global_to_B);
381:   PetscFree(pcis->work_N);
382:   if (pcis->n_neigh > -1) {
383:     ISLocalToGlobalMappingRestoreInfo(pcis->mapping,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));
384:   }
385:   ISLocalToGlobalMappingDestroy(&pcis->mapping);
386:   ISLocalToGlobalMappingDestroy(&pcis->BtoNmap);
387:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetUseStiffnessScaling_C",NULL);
388:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainScalingFactor_C",NULL);
389:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainDiagonalScaling_C",NULL);
390:   return(0);
391: }

393: /* -------------------------------------------------------------------------- */
394: /*
395:    PCISCreate -
396: */
397: PetscErrorCode  PCISCreate(PC pc)
398: {
399:   PC_IS          *pcis = (PC_IS*)(pc->data);

403:   pcis->is_B_local       = 0;
404:   pcis->is_I_local       = 0;
405:   pcis->is_B_global      = 0;
406:   pcis->is_I_global      = 0;
407:   pcis->A_II             = 0;
408:   pcis->A_IB             = 0;
409:   pcis->A_BI             = 0;
410:   pcis->A_BB             = 0;
411:   pcis->D                = 0;
412:   pcis->ksp_N            = 0;
413:   pcis->ksp_D            = 0;
414:   pcis->vec1_N           = 0;
415:   pcis->vec2_N           = 0;
416:   pcis->vec1_D           = 0;
417:   pcis->vec2_D           = 0;
418:   pcis->vec3_D           = 0;
419:   pcis->vec1_B           = 0;
420:   pcis->vec2_B           = 0;
421:   pcis->vec3_B           = 0;
422:   pcis->vec1_global      = 0;
423:   pcis->work_N           = 0;
424:   pcis->global_to_D      = 0;
425:   pcis->N_to_B           = 0;
426:   pcis->N_to_D           = 0;
427:   pcis->global_to_B      = 0;
428:   pcis->mapping          = 0;
429:   pcis->BtoNmap          = 0;
430:   pcis->n_neigh          = -1;
431:   pcis->scaling_factor   = 1.0;
432:   pcis->reusesubmatrices = PETSC_TRUE;
433:   /* composing functions */
434:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetUseStiffnessScaling_C",PCISSetUseStiffnessScaling_IS);
435:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainScalingFactor_C",PCISSetSubdomainScalingFactor_IS);
436:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainDiagonalScaling_C",PCISSetSubdomainDiagonalScaling_IS);
437:   return(0);
438: }

440: /* -------------------------------------------------------------------------- */
441: /*
442:    PCISApplySchur -

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

448:    Output parameters:
449: .  vec1_B - result of Schur complement applied to chunk
450: .  vec2_B - garbage (used as work space), or null (and v is used as workspace)
451: .  vec1_D - garbage (used as work space)
452: .  vec2_D - garbage (used as work space)

454: */
455: PetscErrorCode  PCISApplySchur(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D)
456: {
458:   PC_IS          *pcis = (PC_IS*)(pc->data);

461:   if (!vec2_B) vec2_B = v;

463:   MatMult(pcis->A_BB,v,vec1_B);
464:   MatMult(pcis->A_IB,v,vec1_D);
465:   KSPSolve(pcis->ksp_D,vec1_D,vec2_D);
466:   MatMult(pcis->A_BI,vec2_D,vec2_B);
467:   VecAXPY(vec1_B,-1.0,vec2_B);
468:   return(0);
469: }

471: /* -------------------------------------------------------------------------- */
472: /*
473:    PCISScatterArrayNToVecB - Scatters interface node values from a big array (of all local nodes, interior or interface,
474:    including ghosts) into an interface vector, when in SCATTER_FORWARD mode, or vice-versa, when in SCATTER_REVERSE
475:    mode.

477:    Input parameters:
478: .  pc - preconditioner context
479: .  array_N - [when in SCATTER_FORWARD mode] Array to be scattered into the vector
480: .  v_B - [when in SCATTER_REVERSE mode] Vector to be scattered into the array

482:    Output parameter:
483: .  array_N - [when in SCATTER_REVERSE mode] Array to receive the scattered vector
484: .  v_B - [when in SCATTER_FORWARD mode] Vector to receive the scattered array

486:    Notes:
487:    The entries in the array that do not correspond to interface nodes remain unaltered.
488: */
489: PetscErrorCode  PCISScatterArrayNToVecB(PetscScalar *array_N, Vec v_B, InsertMode imode, ScatterMode smode, PC pc)
490: {
491:   PetscInt       i;
492:   const PetscInt *idex;
494:   PetscScalar    *array_B;
495:   PC_IS          *pcis = (PC_IS*)(pc->data);

498:   VecGetArray(v_B,&array_B);
499:   ISGetIndices(pcis->is_B_local,&idex);

501:   if (smode == SCATTER_FORWARD) {
502:     if (imode == INSERT_VALUES) {
503:       for (i=0; i<pcis->n_B; i++) array_B[i] = array_N[idex[i]];
504:     } else {  /* ADD_VALUES */
505:       for (i=0; i<pcis->n_B; i++) array_B[i] += array_N[idex[i]];
506:     }
507:   } else {  /* SCATTER_REVERSE */
508:     if (imode == INSERT_VALUES) {
509:       for (i=0; i<pcis->n_B; i++) array_N[idex[i]] = array_B[i];
510:     } else {  /* ADD_VALUES */
511:       for (i=0; i<pcis->n_B; i++) array_N[idex[i]] += array_B[i];
512:     }
513:   }
514:   ISRestoreIndices(pcis->is_B_local,&idex);
515:   VecRestoreArray(v_B,&array_B);
516:   return(0);
517: }

519: /* -------------------------------------------------------------------------- */
520: /*
521:    PCISApplyInvSchur - Solves the Neumann problem related to applying the inverse of the Schur complement.
522:    More precisely, solves the problem:
523:                                         [ A_II  A_IB ] [ . ]   [ 0 ]
524:                                         [            ] [   ] = [   ]
525:                                         [ A_BI  A_BB ] [ x ]   [ b ]

527:    Input parameters:
528: .  pc - preconditioner context
529: .  b - vector of local interface nodes (including ghosts)

531:    Output parameters:
532: .  x - vector of local interface nodes (including ghosts); returns the application of the inverse of the Schur
533:        complement to b
534: .  vec1_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space)
535: .  vec2_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space)

537: */
538: PetscErrorCode  PCISApplyInvSchur(PC pc, Vec b, Vec x, Vec vec1_N, Vec vec2_N)
539: {
541:   PC_IS          *pcis = (PC_IS*)(pc->data);

544:   /*
545:     Neumann solvers.
546:     Applying the inverse of the local Schur complement, i.e, solving a Neumann
547:     Problem with zero at the interior nodes of the RHS and extracting the interface
548:     part of the solution. inverse Schur complement is applied to b and the result
549:     is stored in x.
550:   */
551:   /* Setting the RHS vec1_N */
552:   VecSet(vec1_N,0.0);
553:   VecScatterBegin(pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);
554:   VecScatterEnd  (pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);
555:   /* Checking for consistency of the RHS */
556:   {
557:     PetscBool flg = PETSC_FALSE;
558:     PetscOptionsGetBool(NULL,NULL,"-pc_is_check_consistency",&flg,NULL);
559:     if (flg) {
560:       PetscScalar average;
561:       PetscViewer viewer;
562:       PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);

564:       VecSum(vec1_N,&average);
565:       average = average / ((PetscReal)pcis->n);
566:       PetscViewerASCIIPushSynchronized(viewer);
567:       if (pcis->pure_neumann) {
568:         PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is floating. Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));
569:       } else {
570:         PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is fixed.    Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));
571:       }
572:       PetscViewerFlush(viewer);
573:       PetscViewerASCIIPopSynchronized(viewer);
574:     }
575:   }
576:   /* Solving the system for vec2_N */
577:   KSPSolve(pcis->ksp_N,vec1_N,vec2_N);
578:   /* Extracting the local interface vector out of the solution */
579:   VecScatterBegin(pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);
580:   VecScatterEnd  (pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);
581:   return(0);
582: }