Actual source code: pcis.c

petsc-3.4.4 2014-03-13
 2:  #include ../src/ksp/pc/impls/is/pcis.h

  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)
132: {
133:   PC_IS          *pcis  = (PC_IS*)(pc->data);
134:   Mat_IS         *matis = (Mat_IS*)pc->mat->data;
135:   PetscInt       i;
137:   PetscBool      flg;
138:   Vec            counter;

141:   PetscObjectTypeCompare((PetscObject)pc->mat,MATIS,&flg);
142:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Preconditioner type of Neumann Neumman requires matrix of type MATIS");

144:   pcis->pure_neumann = matis->pure_neumann;

146:   /*
147:     Creating the local vector vec1_N, containing the inverse of the number
148:     of subdomains to which each local node (either owned or ghost)
149:     pertains. To accomplish that, we scatter local vectors of 1's to
150:     a global vector (adding the values); scatter the result back to
151:     local vectors and finally invert the result.
152:   */
153:   VecDuplicate(matis->x,&pcis->vec1_N);
154:   MatGetVecs(pc->pmat,&counter,0); /* temporary auxiliar vector */
155:   VecSet(counter,0.0);
156:   VecSet(pcis->vec1_N,1.0);
157:   VecScatterBegin(matis->ctx,pcis->vec1_N,counter,ADD_VALUES,SCATTER_REVERSE);
158:   VecScatterEnd  (matis->ctx,pcis->vec1_N,counter,ADD_VALUES,SCATTER_REVERSE);
159:   VecScatterBegin(matis->ctx,counter,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);
160:   VecScatterEnd  (matis->ctx,counter,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);
161:   /*
162:     Creating local and global index sets for interior and
163:     inteface nodes. Notice that interior nodes have D[i]==1.0.
164:   */
165:   {
166:     PetscInt    n_I;
167:     PetscInt    *idx_I_local,*idx_B_local,*idx_I_global,*idx_B_global;
168:     PetscScalar *array;
169:     /* Identifying interior and interface nodes, in local numbering */
170:     VecGetSize(pcis->vec1_N,&pcis->n);
171:     VecGetArray(pcis->vec1_N,&array);
172:     PetscMalloc(pcis->n*sizeof(PetscInt),&idx_I_local);
173:     PetscMalloc(pcis->n*sizeof(PetscInt),&idx_B_local);
174:     for (i=0, pcis->n_B=0, n_I=0; i<pcis->n; i++) {
175:       if (array[i] == 1.0) {
176:         idx_I_local[n_I] = i;
177:         n_I++;
178:       } else {
179:         idx_B_local[pcis->n_B] = i;
180:         pcis->n_B++;
181:       }
182:     }
183:     /* Getting the global numbering */
184:     idx_B_global = idx_I_local + n_I; /* Just avoiding allocating extra memory, since we have vacant space */
185:     idx_I_global = idx_B_local + pcis->n_B;
186:     ISLocalToGlobalMappingApply(matis->mapping,pcis->n_B,idx_B_local,idx_B_global);
187:     ISLocalToGlobalMappingApply(matis->mapping,n_I,      idx_I_local,idx_I_global);

189:     /* Creating the index sets. */
190:     ISCreateGeneral(MPI_COMM_SELF,pcis->n_B,idx_B_local,PETSC_COPY_VALUES, &pcis->is_B_local);
191:     ISCreateGeneral(MPI_COMM_SELF,pcis->n_B,idx_B_global,PETSC_COPY_VALUES,&pcis->is_B_global);
192:     ISCreateGeneral(MPI_COMM_SELF,n_I,idx_I_local,PETSC_COPY_VALUES, &pcis->is_I_local);
193:     ISCreateGeneral(MPI_COMM_SELF,n_I,idx_I_global,PETSC_COPY_VALUES,&pcis->is_I_global);

195:     /* Freeing memory and restoring arrays */
196:     PetscFree(idx_B_local);
197:     PetscFree(idx_I_local);
198:     VecRestoreArray(pcis->vec1_N,&array);
199:   }

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

205:     [           |      ]
206:     [    A_II   | A_IB ]
207:     A = [           |      ]
208:     [-----------+------]
209:     [    A_BI   | A_BB ]
210:   */

212:   MatGetSubMatrix(matis->A,pcis->is_I_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&pcis->A_II);
213:   MatGetSubMatrix(matis->A,pcis->is_I_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_IB);
214:   MatGetSubMatrix(matis->A,pcis->is_B_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&pcis->A_BI);
215:   MatGetSubMatrix(matis->A,pcis->is_B_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_BB);

217:   /*
218:     Creating work vectors and arrays
219:   */
220:   /* pcis->vec1_N has already been created */
221:   VecDuplicate(pcis->vec1_N,&pcis->vec2_N);
222:   VecCreateSeq(PETSC_COMM_SELF,pcis->n-pcis->n_B,&pcis->vec1_D);
223:   VecDuplicate(pcis->vec1_D,&pcis->vec2_D);
224:   VecDuplicate(pcis->vec1_D,&pcis->vec3_D);
225:   VecCreateSeq(PETSC_COMM_SELF,pcis->n_B,&pcis->vec1_B);
226:   VecDuplicate(pcis->vec1_B,&pcis->vec2_B);
227:   VecDuplicate(pcis->vec1_B,&pcis->vec3_B);
228:   MatGetVecs(pc->pmat,&pcis->vec1_global,0);
229:   PetscMalloc((pcis->n)*sizeof(PetscScalar),&pcis->work_N);

231:   /* Creating the scatter contexts */
232:   VecScatterCreate(pcis->vec1_global,pcis->is_I_global,pcis->vec1_D,(IS)0,&pcis->global_to_D);
233:   VecScatterCreate(pcis->vec1_N,pcis->is_B_local,pcis->vec1_B,(IS)0,&pcis->N_to_B);
234:   VecScatterCreate(pcis->vec1_global,pcis->is_B_global,pcis->vec1_B,(IS)0,&pcis->global_to_B);

236:   /* Creating scaling "matrix" D */
237:   PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_is_use_stiffness_scaling",&pcis->use_stiffness_scaling,NULL);
238:   if (!pcis->D) {
239:     VecDuplicate(pcis->vec1_B,&pcis->D);
240:     if (!pcis->use_stiffness_scaling) {
241:       VecSet(pcis->D,pcis->scaling_factor);
242:     } else {
243:       MatGetDiagonal(matis->A,pcis->vec1_N);
244:       VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);
245:       VecScatterEnd  (pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);
246:     }
247:   }
248:   VecCopy(pcis->D,pcis->vec1_B);
249:   VecSet(counter,0.0);
250:   VecScatterBegin(pcis->global_to_B,pcis->vec1_B,counter,ADD_VALUES,SCATTER_REVERSE);
251:   VecScatterEnd  (pcis->global_to_B,pcis->vec1_B,counter,ADD_VALUES,SCATTER_REVERSE);
252:   VecScatterBegin(pcis->global_to_B,counter,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
253:   VecScatterEnd  (pcis->global_to_B,counter,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
254:   VecPointwiseDivide(pcis->D,pcis->D,pcis->vec1_B);

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

258:   /*
259:     Creating the KSP contexts for the local Dirichlet and Neumann problems.
260:   */
261:   {
262:     PC pc_ctx;
263:     /* Dirichlet */
264:     KSPCreate(PETSC_COMM_SELF,&pcis->ksp_D);
265:     PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_D,(PetscObject)pc,1);
266:     KSPSetOperators(pcis->ksp_D,pcis->A_II,pcis->A_II,SAME_PRECONDITIONER);
267:     KSPSetOptionsPrefix(pcis->ksp_D,"is_localD_");
268:     KSPGetPC(pcis->ksp_D,&pc_ctx);
269:     PCSetType(pc_ctx,PCLU);
270:     KSPSetType(pcis->ksp_D,KSPPREONLY);
271:     KSPSetFromOptions(pcis->ksp_D);
272:     /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
273:     KSPSetUp(pcis->ksp_D);
274:     /* Neumann */
275:     KSPCreate(PETSC_COMM_SELF,&pcis->ksp_N);
276:     PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_N,(PetscObject)pc,1);
277:     KSPSetOperators(pcis->ksp_N,matis->A,matis->A,SAME_PRECONDITIONER);
278:     KSPSetOptionsPrefix(pcis->ksp_N,"is_localN_");
279:     KSPGetPC(pcis->ksp_N,&pc_ctx);
280:     PCSetType(pc_ctx,PCLU);
281:     KSPSetType(pcis->ksp_N,KSPPREONLY);
282:     KSPSetFromOptions(pcis->ksp_N);
283:     {
284:       PetscBool damp_fixed                    = PETSC_FALSE,
285:                 remove_nullspace_fixed        = PETSC_FALSE,
286:                 set_damping_factor_floating   = PETSC_FALSE,
287:                 not_damp_floating             = PETSC_FALSE,
288:                 not_remove_nullspace_floating = PETSC_FALSE;
289:       PetscReal fixed_factor,
290:                 floating_factor;

292:       PetscOptionsGetReal(((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&fixed_factor,&damp_fixed);
293:       if (!damp_fixed) fixed_factor = 0.0;
294:       PetscOptionsGetBool(((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&damp_fixed,NULL);

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

298:       PetscOptionsGetReal(((PetscObject)pc_ctx)->prefix,"-pc_is_set_damping_factor_floating",
299:                               &floating_factor,&set_damping_factor_floating);
300:       if (!set_damping_factor_floating) floating_factor = 0.0;
301:       PetscOptionsGetBool(((PetscObject)pc_ctx)->prefix,"-pc_is_set_damping_factor_floating",&set_damping_factor_floating,NULL);
302:       if (!set_damping_factor_floating) floating_factor = 1.e-12;

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

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

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

336:   ISLocalToGlobalMappingGetInfo(((Mat_IS*)(pc->mat->data))->mapping,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));

338:   pcis->ISLocalToGlobalMappingGetInfoWasCalled = PETSC_TRUE;

340:   VecDestroy(&counter);
341:   return(0);
342: }

344: /* -------------------------------------------------------------------------- */
345: /*
346:    PCISDestroy -
347: */
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->vec1_B);
373:   VecDestroy(&pcis->vec2_B);
374:   VecDestroy(&pcis->vec3_B);
375:   VecDestroy(&pcis->vec1_global);
376:   VecScatterDestroy(&pcis->global_to_D);
377:   VecScatterDestroy(&pcis->N_to_B);
378:   VecScatterDestroy(&pcis->global_to_B);
379:   PetscFree(pcis->work_N);
380:   if (pcis->ISLocalToGlobalMappingGetInfoWasCalled) {
381:     ISLocalToGlobalMappingRestoreInfo((ISLocalToGlobalMapping)0,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));
382:   }
383:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetUseStiffnessScaling_C",NULL);
384:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainScalingFactor_C",NULL);
385:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainDiagonalScaling_C",NULL);
386:   return(0);
387: }

389: /* -------------------------------------------------------------------------- */
390: /*
391:    PCISCreate -
392: */
395: PetscErrorCode  PCISCreate(PC pc)
396: {
397:   PC_IS          *pcis = (PC_IS*)(pc->data);

401:   pcis->is_B_local  = 0;
402:   pcis->is_I_local  = 0;
403:   pcis->is_B_global = 0;
404:   pcis->is_I_global = 0;
405:   pcis->A_II        = 0;
406:   pcis->A_IB        = 0;
407:   pcis->A_BI        = 0;
408:   pcis->A_BB        = 0;
409:   pcis->D           = 0;
410:   pcis->ksp_N       = 0;
411:   pcis->ksp_D      = 0;
412:   pcis->vec1_N      = 0;
413:   pcis->vec2_N      = 0;
414:   pcis->vec1_D      = 0;
415:   pcis->vec2_D      = 0;
416:   pcis->vec3_D      = 0;
417:   pcis->vec1_B      = 0;
418:   pcis->vec2_B      = 0;
419:   pcis->vec3_B      = 0;
420:   pcis->vec1_global = 0;
421:   pcis->work_N      = 0;
422:   pcis->global_to_D = 0;
423:   pcis->N_to_B      = 0;
424:   pcis->global_to_B = 0;

426:   pcis->ISLocalToGlobalMappingGetInfoWasCalled = PETSC_FALSE;

428:   pcis->scaling_factor = 1.0;
429:   /* composing functions */
430:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetUseStiffnessScaling_C",PCISSetUseStiffnessScaling_IS);
431:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainScalingFactor_C",PCISSetSubdomainScalingFactor_IS);
432:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainDiagonalScaling_C",PCISSetSubdomainDiagonalScaling_IS);
433:   return(0);
434: }

436: /* -------------------------------------------------------------------------- */
437: /*
438:    PCISApplySchur -

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

444:    Output parameters:
445: .  vec1_B - result of Schur complement applied to chunk
446: .  vec2_B - garbage (used as work space), or null (and v is used as workspace)
447: .  vec1_D - garbage (used as work space)
448: .  vec2_D - garbage (used as work space)

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

459:   if (!vec2_B) vec2_B = v;

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

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

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

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

484:    Notes:
485:    The entries in the array that do not correspond to interface nodes remain unaltered.
486: */
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: */
540: PetscErrorCode  PCISApplyInvSchur(PC pc, Vec b, Vec x, Vec vec1_N, Vec vec2_N)
541: {
543:   PC_IS          *pcis = (PC_IS*)(pc->data);

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

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