Actual source code: pcis.c
petsc-3.4.5 2014-06-29
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",¬_damp_floating,NULL);
306: PetscOptionsGetBool(((PetscObject)pc_ctx)->prefix,"-pc_is_not_remove_nullspace_floating",¬_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: }