xref: /petsc/src/ksp/pc/impls/is/pcis.c (revision f86f754480c33b36e217a142e1a0364e3fc198eb) 
1 
2 #include <../src/ksp/pc/impls/is/pcis.h> /*I "petscpc.h" I*/
3 
4 static PetscErrorCode PCISSetUseStiffnessScaling_IS(PC pc, PetscBool use)
5 {
6   PC_IS *pcis = (PC_IS*)pc->data;
7 
8   PetscFunctionBegin;
9   pcis->use_stiffness_scaling = use;
10   PetscFunctionReturn(0);
11 }
12 
13 /*@
14  PCISSetUseStiffnessScaling - Tells PCIS to construct partition of unity using
15                               local matrices' diagonal.
16 
17    Not collective
18 
19    Input Parameters:
20 +  pc - the preconditioning context
21 -  use - whether or not pcis use matrix diagonal to build partition of unity.
22 
23    Level: intermediate
24 
25    Notes:
26 
27 .seealso: PCBDDC
28 @*/
29 PetscErrorCode PCISSetUseStiffnessScaling(PC pc, PetscBool use)
30 {
31   PetscErrorCode ierr;
32 
33   PetscFunctionBegin;
34   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
35   ierr = PetscTryMethod(pc,"PCISSetUseStiffnessScaling_C",(PC,PetscBool),(pc,use));CHKERRQ(ierr);
36   PetscFunctionReturn(0);
37 }
38 
39 static PetscErrorCode PCISSetSubdomainDiagonalScaling_IS(PC pc, Vec scaling_factors)
40 {
41   PetscErrorCode ierr;
42   PC_IS          *pcis = (PC_IS*)pc->data;
43 
44   PetscFunctionBegin;
45   ierr    = PetscObjectReference((PetscObject)scaling_factors);CHKERRQ(ierr);
46   ierr    = VecDestroy(&pcis->D);CHKERRQ(ierr);
47   pcis->D = scaling_factors;
48   PetscFunctionReturn(0);
49 }
50 
51 /*@
52  PCISSetSubdomainDiagonalScaling - Set diagonal scaling for PCIS.
53 
54    Not collective
55 
56    Input Parameters:
57 +  pc - the preconditioning context
58 -  scaling_factors - scaling factors for the subdomain
59 
60    Level: intermediate
61 
62    Notes:
63    Intended to use with jumping coefficients cases.
64 
65 .seealso: PCBDDC
66 @*/
67 PetscErrorCode PCISSetSubdomainDiagonalScaling(PC pc, Vec scaling_factors)
68 {
69   PetscErrorCode ierr;
70 
71   PetscFunctionBegin;
72   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
73   ierr = PetscTryMethod(pc,"PCISSetSubdomainDiagonalScaling_C",(PC,Vec),(pc,scaling_factors));CHKERRQ(ierr);
74   PetscFunctionReturn(0);
75 }
76 
77 static PetscErrorCode PCISSetSubdomainScalingFactor_IS(PC pc, PetscScalar scal)
78 {
79   PC_IS *pcis = (PC_IS*)pc->data;
80 
81   PetscFunctionBegin;
82   pcis->scaling_factor = scal;
83   PetscFunctionReturn(0);
84 }
85 
86 /*@
87  PCISSetSubdomainScalingFactor - Set scaling factor for PCIS.
88 
89    Not collective
90 
91    Input Parameters:
92 +  pc - the preconditioning context
93 -  scal - scaling factor for the subdomain
94 
95    Level: intermediate
96 
97    Notes:
98    Intended to use with jumping coefficients cases.
99 
100 .seealso: PCBDDC
101 @*/
102 PetscErrorCode PCISSetSubdomainScalingFactor(PC pc, PetscScalar scal)
103 {
104   PetscErrorCode ierr;
105 
106   PetscFunctionBegin;
107   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
108   ierr = PetscTryMethod(pc,"PCISSetSubdomainScalingFactor_C",(PC,PetscScalar),(pc,scal));CHKERRQ(ierr);
109   PetscFunctionReturn(0);
110 }
111 
112 
113 /* -------------------------------------------------------------------------- */
114 /*
115    PCISSetUp -
116 */
117 PetscErrorCode  PCISSetUp(PC pc, PetscBool computematrices, PetscBool computesolvers)
118 {
119   PC_IS          *pcis  = (PC_IS*)(pc->data);
120   Mat_IS         *matis;
121   MatReuse       reuse;
122   PetscErrorCode ierr;
123   PetscBool      flg,issbaij;
124 
125   PetscFunctionBegin;
126   ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATIS,&flg);CHKERRQ(ierr);
127   if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Requires preconditioning matrix of type MATIS");
128   matis = (Mat_IS*)pc->pmat->data;
129   if (pc->useAmat) {
130     ierr = PetscObjectTypeCompare((PetscObject)pc->mat,MATIS,&flg);CHKERRQ(ierr);
131     if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Requires linear system matrix of type MATIS");
132   }
133 
134   /* first time creation, get info on substructuring */
135   if (!pc->setupcalled) {
136     PetscInt    n_I;
137     PetscInt    *idx_I_local,*idx_B_local,*idx_I_global,*idx_B_global;
138     PetscBT     bt;
139     PetscInt    i,j;
140 
141     /* get info on mapping */
142     ierr = PetscObjectReference((PetscObject)pc->pmat->rmap->mapping);CHKERRQ(ierr);
143     ierr = ISLocalToGlobalMappingDestroy(&pcis->mapping);CHKERRQ(ierr);
144     pcis->mapping = pc->pmat->rmap->mapping;
145     ierr = ISLocalToGlobalMappingGetSize(pcis->mapping,&pcis->n);CHKERRQ(ierr);
146     ierr = ISLocalToGlobalMappingGetInfo(pcis->mapping,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));CHKERRQ(ierr);
147 
148     /* Identifying interior and interface nodes, in local numbering */
149     ierr = PetscBTCreate(pcis->n,&bt);CHKERRQ(ierr);
150     for (i=0;i<pcis->n_neigh;i++)
151       for (j=0;j<pcis->n_shared[i];j++) {
152           ierr = PetscBTSet(bt,pcis->shared[i][j]);CHKERRQ(ierr);
153       }
154 
155     /* Creating local and global index sets for interior and inteface nodes. */
156     ierr = PetscMalloc1(pcis->n,&idx_I_local);CHKERRQ(ierr);
157     ierr = PetscMalloc1(pcis->n,&idx_B_local);CHKERRQ(ierr);
158     for (i=0, pcis->n_B=0, n_I=0; i<pcis->n; i++) {
159       if (!PetscBTLookup(bt,i)) {
160         idx_I_local[n_I] = i;
161         n_I++;
162       } else {
163         idx_B_local[pcis->n_B] = i;
164         pcis->n_B++;
165       }
166     }
167 
168     /* Getting the global numbering */
169     idx_B_global = idx_I_local + n_I; /* Just avoiding allocating extra memory, since we have vacant space */
170     idx_I_global = idx_B_local + pcis->n_B;
171     ierr         = ISLocalToGlobalMappingApply(pcis->mapping,pcis->n_B,idx_B_local,idx_B_global);CHKERRQ(ierr);
172     ierr         = ISLocalToGlobalMappingApply(pcis->mapping,n_I,idx_I_local,idx_I_global);CHKERRQ(ierr);
173 
174     /* Creating the index sets */
175     ierr = ISCreateGeneral(PETSC_COMM_SELF,pcis->n_B,idx_B_local,PETSC_COPY_VALUES, &pcis->is_B_local);CHKERRQ(ierr);
176     ierr = ISCreateGeneral(PetscObjectComm((PetscObject)pc),pcis->n_B,idx_B_global,PETSC_COPY_VALUES,&pcis->is_B_global);CHKERRQ(ierr);
177     ierr = ISCreateGeneral(PETSC_COMM_SELF,n_I,idx_I_local,PETSC_COPY_VALUES, &pcis->is_I_local);CHKERRQ(ierr);
178     ierr = ISCreateGeneral(PetscObjectComm((PetscObject)pc),n_I,idx_I_global,PETSC_COPY_VALUES,&pcis->is_I_global);CHKERRQ(ierr);
179 
180     /* Freeing memory */
181     ierr = PetscFree(idx_B_local);CHKERRQ(ierr);
182     ierr = PetscFree(idx_I_local);CHKERRQ(ierr);
183     ierr = PetscBTDestroy(&bt);CHKERRQ(ierr);
184 
185     /* Creating work vectors and arrays */
186     ierr = VecDuplicate(matis->x,&pcis->vec1_N);CHKERRQ(ierr);
187     ierr = VecDuplicate(pcis->vec1_N,&pcis->vec2_N);CHKERRQ(ierr);
188     ierr = VecCreate(PETSC_COMM_SELF,&pcis->vec1_D);CHKERRQ(ierr);
189     ierr = VecSetSizes(pcis->vec1_D,pcis->n-pcis->n_B,PETSC_DECIDE);CHKERRQ(ierr);
190     ierr = VecSetType(pcis->vec1_D,((PetscObject)pcis->vec1_N)->type_name);CHKERRQ(ierr);
191     ierr = VecDuplicate(pcis->vec1_D,&pcis->vec2_D);CHKERRQ(ierr);
192     ierr = VecDuplicate(pcis->vec1_D,&pcis->vec3_D);CHKERRQ(ierr);
193     ierr = VecDuplicate(pcis->vec1_D,&pcis->vec4_D);CHKERRQ(ierr);
194     ierr = VecCreate(PETSC_COMM_SELF,&pcis->vec1_B);CHKERRQ(ierr);
195     ierr = VecSetSizes(pcis->vec1_B,pcis->n_B,PETSC_DECIDE);CHKERRQ(ierr);
196     ierr = VecSetType(pcis->vec1_B,((PetscObject)pcis->vec1_N)->type_name);CHKERRQ(ierr);
197     ierr = VecDuplicate(pcis->vec1_B,&pcis->vec2_B);CHKERRQ(ierr);
198     ierr = VecDuplicate(pcis->vec1_B,&pcis->vec3_B);CHKERRQ(ierr);
199     ierr = MatCreateVecs(pc->pmat,&pcis->vec1_global,0);CHKERRQ(ierr);
200     ierr = PetscMalloc1(pcis->n,&pcis->work_N);CHKERRQ(ierr);
201     /* scaling vector */
202     if (!pcis->D) { /* it can happen that the user passed in a scaling vector via PCISSetSubdomainDiagonalScaling */
203       ierr = VecDuplicate(pcis->vec1_B,&pcis->D);CHKERRQ(ierr);
204       ierr = VecSet(pcis->D,pcis->scaling_factor);CHKERRQ(ierr);
205     }
206 
207     /* Creating the scatter contexts */
208     ierr = VecScatterCreateWithData(pcis->vec1_N,pcis->is_I_local,pcis->vec1_D,(IS)0,&pcis->N_to_D);CHKERRQ(ierr);
209     ierr = VecScatterCreateWithData(pcis->vec1_global,pcis->is_I_global,pcis->vec1_D,(IS)0,&pcis->global_to_D);CHKERRQ(ierr);
210     ierr = VecScatterCreateWithData(pcis->vec1_N,pcis->is_B_local,pcis->vec1_B,(IS)0,&pcis->N_to_B);CHKERRQ(ierr);
211     ierr = VecScatterCreateWithData(pcis->vec1_global,pcis->is_B_global,pcis->vec1_B,(IS)0,&pcis->global_to_B);CHKERRQ(ierr);
212 
213     /* map from boundary to local */
214     ierr = ISLocalToGlobalMappingCreateIS(pcis->is_B_local,&pcis->BtoNmap);CHKERRQ(ierr);
215   }
216 
217   /*
218     Extracting the blocks A_II, A_BI, A_IB and A_BB from A. If the numbering
219     is such that interior nodes come first than the interface ones, we have
220 
221         [ A_II | A_IB ]
222     A = [------+------]
223         [ A_BI | A_BB ]
224   */
225   if (computematrices) {
226     PetscBool amat = (PetscBool)(pc->mat != pc->pmat && pc->useAmat);
227     PetscInt  bs,ibs;
228 
229     reuse = MAT_INITIAL_MATRIX;
230     if (pcis->reusesubmatrices && pc->setupcalled) {
231       if (pc->flag == SAME_NONZERO_PATTERN) {
232         reuse = MAT_REUSE_MATRIX;
233       } else {
234         reuse = MAT_INITIAL_MATRIX;
235       }
236     }
237     if (reuse == MAT_INITIAL_MATRIX) {
238       ierr = MatDestroy(&pcis->A_II);CHKERRQ(ierr);
239       ierr = MatDestroy(&pcis->pA_II);CHKERRQ(ierr);
240       ierr = MatDestroy(&pcis->A_IB);CHKERRQ(ierr);
241       ierr = MatDestroy(&pcis->A_BI);CHKERRQ(ierr);
242       ierr = MatDestroy(&pcis->A_BB);CHKERRQ(ierr);
243     }
244 
245     ierr = ISLocalToGlobalMappingGetBlockSize(pcis->mapping,&ibs);CHKERRQ(ierr);
246     ierr = MatGetBlockSize(matis->A,&bs);CHKERRQ(ierr);
247     ierr = MatCreateSubMatrix(matis->A,pcis->is_I_local,pcis->is_I_local,reuse,&pcis->pA_II);CHKERRQ(ierr);
248     if (amat) {
249       Mat_IS *amatis = (Mat_IS*)pc->mat->data;
250       ierr = MatCreateSubMatrix(amatis->A,pcis->is_I_local,pcis->is_I_local,reuse,&pcis->A_II);CHKERRQ(ierr);
251     } else {
252       ierr = PetscObjectReference((PetscObject)pcis->pA_II);CHKERRQ(ierr);
253       ierr = MatDestroy(&pcis->A_II);CHKERRQ(ierr);
254       pcis->A_II = pcis->pA_II;
255     }
256     ierr = MatSetBlockSize(pcis->A_II,bs == ibs ? bs : 1);CHKERRQ(ierr);
257     ierr = MatSetBlockSize(pcis->pA_II,bs == ibs ? bs : 1);CHKERRQ(ierr);
258     ierr = MatCreateSubMatrix(matis->A,pcis->is_B_local,pcis->is_B_local,reuse,&pcis->A_BB);CHKERRQ(ierr);
259     ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);CHKERRQ(ierr);
260     if (!issbaij) {
261       ierr = MatCreateSubMatrix(matis->A,pcis->is_I_local,pcis->is_B_local,reuse,&pcis->A_IB);CHKERRQ(ierr);
262       ierr = MatCreateSubMatrix(matis->A,pcis->is_B_local,pcis->is_I_local,reuse,&pcis->A_BI);CHKERRQ(ierr);
263     } else {
264       Mat newmat;
265 
266       ierr = MatConvert(matis->A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&newmat);CHKERRQ(ierr);
267       ierr = MatCreateSubMatrix(newmat,pcis->is_I_local,pcis->is_B_local,reuse,&pcis->A_IB);CHKERRQ(ierr);
268       ierr = MatCreateSubMatrix(newmat,pcis->is_B_local,pcis->is_I_local,reuse,&pcis->A_BI);CHKERRQ(ierr);
269       ierr = MatDestroy(&newmat);CHKERRQ(ierr);
270     }
271     ierr = MatSetBlockSize(pcis->A_BB,bs == ibs ? bs : 1);CHKERRQ(ierr);
272   }
273 
274   /* Creating scaling vector D */
275   ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_is_use_stiffness_scaling",&pcis->use_stiffness_scaling,NULL);CHKERRQ(ierr);
276   if (pcis->use_stiffness_scaling) {
277     PetscScalar *a;
278     PetscInt    i,n;
279 
280     if (pcis->A_BB) {
281       ierr = MatGetDiagonal(pcis->A_BB,pcis->D);CHKERRQ(ierr);
282     } else {
283       ierr = MatGetDiagonal(matis->A,pcis->vec1_N);CHKERRQ(ierr);
284       ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
285       ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
286     }
287     ierr = VecAbs(pcis->D);CHKERRQ(ierr);
288     ierr = VecGetLocalSize(pcis->D,&n);CHKERRQ(ierr);
289     ierr = VecGetArray(pcis->D,&a);CHKERRQ(ierr);
290     for (i=0;i<n;i++) if (PetscAbsScalar(a[i])<PETSC_SMALL) a[i] = 1.0;
291     ierr = VecRestoreArray(pcis->D,&a);CHKERRQ(ierr);
292   }
293   ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr);
294   ierr = VecScatterBegin(pcis->global_to_B,pcis->D,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
295   ierr = VecScatterEnd(pcis->global_to_B,pcis->D,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
296   ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
297   ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
298   ierr = VecPointwiseDivide(pcis->D,pcis->D,pcis->vec1_B);CHKERRQ(ierr);
299 
300   /* See historical note 01, at the bottom of this file. */
301 
302   /* Creating the KSP contexts for the local Dirichlet and Neumann problems */
303   if (computesolvers) {
304     PC pc_ctx;
305 
306     pcis->pure_neumann = matis->pure_neumann;
307     /* Dirichlet */
308     ierr = KSPCreate(PETSC_COMM_SELF,&pcis->ksp_D);CHKERRQ(ierr);
309     ierr = KSPSetErrorIfNotConverged(pcis->ksp_D,pc->erroriffailure);CHKERRQ(ierr);
310     ierr = PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_D,(PetscObject)pc,1);CHKERRQ(ierr);
311     ierr = KSPSetOperators(pcis->ksp_D,pcis->A_II,pcis->A_II);CHKERRQ(ierr);
312     ierr = KSPSetOptionsPrefix(pcis->ksp_D,"is_localD_");CHKERRQ(ierr);
313     ierr = KSPGetPC(pcis->ksp_D,&pc_ctx);CHKERRQ(ierr);
314     ierr = PCSetType(pc_ctx,PCLU);CHKERRQ(ierr);
315     ierr = KSPSetType(pcis->ksp_D,KSPPREONLY);CHKERRQ(ierr);
316     ierr = KSPSetFromOptions(pcis->ksp_D);CHKERRQ(ierr);
317     /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
318     ierr = KSPSetUp(pcis->ksp_D);CHKERRQ(ierr);
319     /* Neumann */
320     ierr = KSPCreate(PETSC_COMM_SELF,&pcis->ksp_N);CHKERRQ(ierr);
321     ierr = KSPSetErrorIfNotConverged(pcis->ksp_N,pc->erroriffailure);CHKERRQ(ierr);
322     ierr = PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_N,(PetscObject)pc,1);CHKERRQ(ierr);
323     ierr = KSPSetOperators(pcis->ksp_N,matis->A,matis->A);CHKERRQ(ierr);
324     ierr = KSPSetOptionsPrefix(pcis->ksp_N,"is_localN_");CHKERRQ(ierr);
325     ierr = KSPGetPC(pcis->ksp_N,&pc_ctx);CHKERRQ(ierr);
326     ierr = PCSetType(pc_ctx,PCLU);CHKERRQ(ierr);
327     ierr = KSPSetType(pcis->ksp_N,KSPPREONLY);CHKERRQ(ierr);
328     ierr = KSPSetFromOptions(pcis->ksp_N);CHKERRQ(ierr);
329     {
330       PetscBool damp_fixed                    = PETSC_FALSE,
331                 remove_nullspace_fixed        = PETSC_FALSE,
332                 set_damping_factor_floating   = PETSC_FALSE,
333                 not_damp_floating             = PETSC_FALSE,
334                 not_remove_nullspace_floating = PETSC_FALSE;
335       PetscReal fixed_factor,
336                 floating_factor;
337 
338       ierr = PetscOptionsGetReal(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&fixed_factor,&damp_fixed);CHKERRQ(ierr);
339       if (!damp_fixed) fixed_factor = 0.0;
340       ierr = PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&damp_fixed,NULL);CHKERRQ(ierr);
341 
342       ierr = PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_remove_nullspace_fixed",&remove_nullspace_fixed,NULL);CHKERRQ(ierr);
343 
344       ierr = PetscOptionsGetReal(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_set_damping_factor_floating",
345                               &floating_factor,&set_damping_factor_floating);CHKERRQ(ierr);
346       if (!set_damping_factor_floating) floating_factor = 0.0;
347       ierr = PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_set_damping_factor_floating",&set_damping_factor_floating,NULL);CHKERRQ(ierr);
348       if (!set_damping_factor_floating) floating_factor = 1.e-12;
349 
350       ierr = PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_not_damp_floating",&not_damp_floating,NULL);CHKERRQ(ierr);
351 
352       ierr = PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_not_remove_nullspace_floating",&not_remove_nullspace_floating,NULL);CHKERRQ(ierr);
353 
354       if (pcis->pure_neumann) {  /* floating subdomain */
355         if (!(not_damp_floating)) {
356           ierr = PCFactorSetShiftType(pc_ctx,MAT_SHIFT_NONZERO);CHKERRQ(ierr);
357           ierr = PCFactorSetShiftAmount(pc_ctx,floating_factor);CHKERRQ(ierr);
358         }
359         if (!(not_remove_nullspace_floating)) {
360           MatNullSpace nullsp;
361           ierr = MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_TRUE,0,NULL,&nullsp);CHKERRQ(ierr);
362           ierr = MatSetNullSpace(matis->A,nullsp);CHKERRQ(ierr);
363           ierr = MatNullSpaceDestroy(&nullsp);CHKERRQ(ierr);
364         }
365       } else {  /* fixed subdomain */
366         if (damp_fixed) {
367           ierr = PCFactorSetShiftType(pc_ctx,MAT_SHIFT_NONZERO);CHKERRQ(ierr);
368           ierr = PCFactorSetShiftAmount(pc_ctx,floating_factor);CHKERRQ(ierr);
369         }
370         if (remove_nullspace_fixed) {
371           MatNullSpace nullsp;
372           ierr = MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_TRUE,0,NULL,&nullsp);CHKERRQ(ierr);
373           ierr = MatSetNullSpace(matis->A,nullsp);CHKERRQ(ierr);
374           ierr = MatNullSpaceDestroy(&nullsp);CHKERRQ(ierr);
375         }
376       }
377     }
378     /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
379     ierr = KSPSetUp(pcis->ksp_N);CHKERRQ(ierr);
380   }
381   PetscFunctionReturn(0);
382 }
383 
384 /* -------------------------------------------------------------------------- */
385 /*
386    PCISDestroy -
387 */
388 PetscErrorCode  PCISDestroy(PC pc)
389 {
390   PC_IS          *pcis = (PC_IS*)(pc->data);
391   PetscErrorCode ierr;
392 
393   PetscFunctionBegin;
394   ierr = ISDestroy(&pcis->is_B_local);CHKERRQ(ierr);
395   ierr = ISDestroy(&pcis->is_I_local);CHKERRQ(ierr);
396   ierr = ISDestroy(&pcis->is_B_global);CHKERRQ(ierr);
397   ierr = ISDestroy(&pcis->is_I_global);CHKERRQ(ierr);
398   ierr = MatDestroy(&pcis->A_II);CHKERRQ(ierr);
399   ierr = MatDestroy(&pcis->pA_II);CHKERRQ(ierr);
400   ierr = MatDestroy(&pcis->A_IB);CHKERRQ(ierr);
401   ierr = MatDestroy(&pcis->A_BI);CHKERRQ(ierr);
402   ierr = MatDestroy(&pcis->A_BB);CHKERRQ(ierr);
403   ierr = VecDestroy(&pcis->D);CHKERRQ(ierr);
404   ierr = KSPDestroy(&pcis->ksp_N);CHKERRQ(ierr);
405   ierr = KSPDestroy(&pcis->ksp_D);CHKERRQ(ierr);
406   ierr = VecDestroy(&pcis->vec1_N);CHKERRQ(ierr);
407   ierr = VecDestroy(&pcis->vec2_N);CHKERRQ(ierr);
408   ierr = VecDestroy(&pcis->vec1_D);CHKERRQ(ierr);
409   ierr = VecDestroy(&pcis->vec2_D);CHKERRQ(ierr);
410   ierr = VecDestroy(&pcis->vec3_D);CHKERRQ(ierr);
411   ierr = VecDestroy(&pcis->vec4_D);CHKERRQ(ierr);
412   ierr = VecDestroy(&pcis->vec1_B);CHKERRQ(ierr);
413   ierr = VecDestroy(&pcis->vec2_B);CHKERRQ(ierr);
414   ierr = VecDestroy(&pcis->vec3_B);CHKERRQ(ierr);
415   ierr = VecDestroy(&pcis->vec1_global);CHKERRQ(ierr);
416   ierr = VecScatterDestroy(&pcis->global_to_D);CHKERRQ(ierr);
417   ierr = VecScatterDestroy(&pcis->N_to_B);CHKERRQ(ierr);
418   ierr = VecScatterDestroy(&pcis->N_to_D);CHKERRQ(ierr);
419   ierr = VecScatterDestroy(&pcis->global_to_B);CHKERRQ(ierr);
420   ierr = PetscFree(pcis->work_N);CHKERRQ(ierr);
421   if (pcis->n_neigh > -1) {
422     ierr = ISLocalToGlobalMappingRestoreInfo(pcis->mapping,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));CHKERRQ(ierr);
423   }
424   ierr = ISLocalToGlobalMappingDestroy(&pcis->mapping);CHKERRQ(ierr);
425   ierr = ISLocalToGlobalMappingDestroy(&pcis->BtoNmap);CHKERRQ(ierr);
426   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCISSetUseStiffnessScaling_C",NULL);CHKERRQ(ierr);
427   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainScalingFactor_C",NULL);CHKERRQ(ierr);
428   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainDiagonalScaling_C",NULL);CHKERRQ(ierr);
429   PetscFunctionReturn(0);
430 }
431 
432 /* -------------------------------------------------------------------------- */
433 /*
434    PCISCreate -
435 */
436 PetscErrorCode  PCISCreate(PC pc)
437 {
438   PC_IS          *pcis = (PC_IS*)(pc->data);
439   PetscErrorCode ierr;
440 
441   PetscFunctionBegin;
442   pcis->n_neigh          = -1;
443   pcis->scaling_factor   = 1.0;
444   pcis->reusesubmatrices = PETSC_TRUE;
445   /* composing functions */
446   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCISSetUseStiffnessScaling_C",PCISSetUseStiffnessScaling_IS);CHKERRQ(ierr);
447   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainScalingFactor_C",PCISSetSubdomainScalingFactor_IS);CHKERRQ(ierr);
448   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainDiagonalScaling_C",PCISSetSubdomainDiagonalScaling_IS);CHKERRQ(ierr);
449   PetscFunctionReturn(0);
450 }
451 
452 /* -------------------------------------------------------------------------- */
453 /*
454    PCISApplySchur -
455 
456    Input parameters:
457 .  pc - preconditioner context
458 .  v - vector to which the Schur complement is to be applied (it is NOT modified inside this function, UNLESS vec2_B is null)
459 
460    Output parameters:
461 .  vec1_B - result of Schur complement applied to chunk
462 .  vec2_B - garbage (used as work space), or null (and v is used as workspace)
463 .  vec1_D - garbage (used as work space)
464 .  vec2_D - garbage (used as work space)
465 
466 */
467 PetscErrorCode  PCISApplySchur(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D)
468 {
469   PetscErrorCode ierr;
470   PC_IS          *pcis = (PC_IS*)(pc->data);
471 
472   PetscFunctionBegin;
473   if (!vec2_B) vec2_B = v;
474 
475   ierr = MatMult(pcis->A_BB,v,vec1_B);CHKERRQ(ierr);
476   ierr = MatMult(pcis->A_IB,v,vec1_D);CHKERRQ(ierr);
477   ierr = KSPSolve(pcis->ksp_D,vec1_D,vec2_D);CHKERRQ(ierr);
478   ierr = KSPCheckSolve(pcis->ksp_D,pc,vec2_D);CHKERRQ(ierr);
479   ierr = MatMult(pcis->A_BI,vec2_D,vec2_B);CHKERRQ(ierr);
480   ierr = VecAXPY(vec1_B,-1.0,vec2_B);CHKERRQ(ierr);
481   PetscFunctionReturn(0);
482 }
483 
484 /* -------------------------------------------------------------------------- */
485 /*
486    PCISScatterArrayNToVecB - Scatters interface node values from a big array (of all local nodes, interior or interface,
487    including ghosts) into an interface vector, when in SCATTER_FORWARD mode, or vice-versa, when in SCATTER_REVERSE
488    mode.
489 
490    Input parameters:
491 .  pc - preconditioner context
492 .  array_N - [when in SCATTER_FORWARD mode] Array to be scattered into the vector
493 .  v_B - [when in SCATTER_REVERSE mode] Vector to be scattered into the array
494 
495    Output parameter:
496 .  array_N - [when in SCATTER_REVERSE mode] Array to receive the scattered vector
497 .  v_B - [when in SCATTER_FORWARD mode] Vector to receive the scattered array
498 
499    Notes:
500    The entries in the array that do not correspond to interface nodes remain unaltered.
501 */
502 PetscErrorCode  PCISScatterArrayNToVecB(PetscScalar *array_N, Vec v_B, InsertMode imode, ScatterMode smode, PC pc)
503 {
504   PetscInt       i;
505   const PetscInt *idex;
506   PetscErrorCode ierr;
507   PetscScalar    *array_B;
508   PC_IS          *pcis = (PC_IS*)(pc->data);
509 
510   PetscFunctionBegin;
511   ierr = VecGetArray(v_B,&array_B);CHKERRQ(ierr);
512   ierr = ISGetIndices(pcis->is_B_local,&idex);CHKERRQ(ierr);
513 
514   if (smode == SCATTER_FORWARD) {
515     if (imode == INSERT_VALUES) {
516       for (i=0; i<pcis->n_B; i++) array_B[i] = array_N[idex[i]];
517     } else {  /* ADD_VALUES */
518       for (i=0; i<pcis->n_B; i++) array_B[i] += array_N[idex[i]];
519     }
520   } else {  /* SCATTER_REVERSE */
521     if (imode == INSERT_VALUES) {
522       for (i=0; i<pcis->n_B; i++) array_N[idex[i]] = array_B[i];
523     } else {  /* ADD_VALUES */
524       for (i=0; i<pcis->n_B; i++) array_N[idex[i]] += array_B[i];
525     }
526   }
527   ierr = ISRestoreIndices(pcis->is_B_local,&idex);CHKERRQ(ierr);
528   ierr = VecRestoreArray(v_B,&array_B);CHKERRQ(ierr);
529   PetscFunctionReturn(0);
530 }
531 
532 /* -------------------------------------------------------------------------- */
533 /*
534    PCISApplyInvSchur - Solves the Neumann problem related to applying the inverse of the Schur complement.
535    More precisely, solves the problem:
536                                         [ A_II  A_IB ] [ . ]   [ 0 ]
537                                         [            ] [   ] = [   ]
538                                         [ A_BI  A_BB ] [ x ]   [ b ]
539 
540    Input parameters:
541 .  pc - preconditioner context
542 .  b - vector of local interface nodes (including ghosts)
543 
544    Output parameters:
545 .  x - vector of local interface nodes (including ghosts); returns the application of the inverse of the Schur
546        complement to b
547 .  vec1_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space)
548 .  vec2_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space)
549 
550 */
551 PetscErrorCode  PCISApplyInvSchur(PC pc, Vec b, Vec x, Vec vec1_N, Vec vec2_N)
552 {
553   PetscErrorCode ierr;
554   PC_IS          *pcis = (PC_IS*)(pc->data);
555 
556   PetscFunctionBegin;
557   /*
558     Neumann solvers.
559     Applying the inverse of the local Schur complement, i.e, solving a Neumann
560     Problem with zero at the interior nodes of the RHS and extracting the interface
561     part of the solution. inverse Schur complement is applied to b and the result
562     is stored in x.
563   */
564   /* Setting the RHS vec1_N */
565   ierr = VecSet(vec1_N,0.0);CHKERRQ(ierr);
566   ierr = VecScatterBegin(pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
567   ierr = VecScatterEnd  (pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
568   /* Checking for consistency of the RHS */
569   {
570     PetscBool flg = PETSC_FALSE;
571     ierr = PetscOptionsGetBool(NULL,NULL,"-pc_is_check_consistency",&flg,NULL);CHKERRQ(ierr);
572     if (flg) {
573       PetscScalar average;
574       PetscViewer viewer;
575       ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);CHKERRQ(ierr);
576 
577       ierr    = VecSum(vec1_N,&average);CHKERRQ(ierr);
578       average = average / ((PetscReal)pcis->n);
579       ierr    = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr);
580       if (pcis->pure_neumann) {
581         ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is floating. Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));CHKERRQ(ierr);
582       } else {
583         ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is fixed.    Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));CHKERRQ(ierr);
584       }
585       ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
586       ierr = PetscViewerASCIIPopSynchronized(viewer);CHKERRQ(ierr);
587     }
588   }
589   /* Solving the system for vec2_N */
590   ierr = KSPSolve(pcis->ksp_N,vec1_N,vec2_N);CHKERRQ(ierr);
591   ierr = KSPCheckSolve(pcis->ksp_N,pc,vec2_N);CHKERRQ(ierr);
592   /* Extracting the local interface vector out of the solution */
593   ierr = VecScatterBegin(pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
594   ierr = VecScatterEnd  (pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
595   PetscFunctionReturn(0);
596 }
597