xref: /petsc/src/ksp/pc/impls/bddc/bddcnullspace.c (revision 98a51de6e8e383a56523efa208c0ffa3caed8230) 
1 #include "bddc.h"
2 #include "bddcprivate.h"
3 
4 #undef __FUNCT__
5 #define __FUNCT__ "PCBDDCNullSpaceAssembleCoarse"
6 PetscErrorCode PCBDDCNullSpaceAssembleCoarse(PC pc, Mat coarse_mat, MatNullSpace* CoarseNullSpace)
7 {
8   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
9   Mat_IS         *matis = (Mat_IS*)pc->pmat->data;
10   MatNullSpace   tempCoarseNullSpace=NULL;
11   const Vec      *nsp_vecs;
12   Vec            *coarse_nsp_vecs,local_vec,local_primal_vec,wcoarse_vec,wcoarse_rhs;
13   PetscInt       nsp_size,coarse_nsp_size,i;
14   PetscBool      nsp_has_cnst;
15   PetscReal      test_null;
16   PetscErrorCode ierr;
17 
18   PetscFunctionBegin;
19   tempCoarseNullSpace = 0;
20   coarse_nsp_size = 0;
21   coarse_nsp_vecs = 0;
22   ierr = MatNullSpaceGetVecs(pcbddc->NullSpace,&nsp_has_cnst,&nsp_size,&nsp_vecs);CHKERRQ(ierr);
23   if (coarse_mat) {
24     ierr = PetscMalloc((nsp_size+1)*sizeof(Vec),&coarse_nsp_vecs);CHKERRQ(ierr);
25     for (i=0;i<nsp_size+1;i++) {
26       ierr = MatGetVecs(coarse_mat,&coarse_nsp_vecs[i],NULL);CHKERRQ(ierr);
27     }
28     if (pcbddc->dbg_flag) {
29       ierr = MatGetVecs(coarse_mat,&wcoarse_vec,&wcoarse_rhs);CHKERRQ(ierr);
30     }
31   }
32   ierr = MatGetVecs(pcbddc->ConstraintMatrix,&local_vec,&local_primal_vec);CHKERRQ(ierr);
33   if (nsp_has_cnst) {
34     ierr = VecSet(local_vec,1.0);CHKERRQ(ierr);
35     ierr = MatMult(pcbddc->ConstraintMatrix,local_vec,local_primal_vec);CHKERRQ(ierr);
36     ierr = VecScatterBegin(pcbddc->coarse_loc_to_glob,local_primal_vec,pcbddc->coarse_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
37     ierr = VecScatterEnd(pcbddc->coarse_loc_to_glob,local_primal_vec,pcbddc->coarse_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
38     if (coarse_mat) {
39       PetscScalar *array_out;
40       const PetscScalar *array_in;
41       PetscInt lsize;
42       if (pcbddc->dbg_flag) {
43         PetscViewer dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)coarse_mat));
44         ierr = MatMult(coarse_mat,wcoarse_vec,wcoarse_rhs);CHKERRQ(ierr);
45         ierr = VecNorm(wcoarse_rhs,NORM_INFINITY,&test_null);CHKERRQ(ierr);
46         ierr = PetscViewerASCIIPrintf(dbg_viewer,"Constant coarse null space error % 1.14e\n",test_null);CHKERRQ(ierr);
47         ierr = PetscViewerFlush(dbg_viewer);CHKERRQ(ierr);
48       }
49       ierr = VecGetLocalSize(pcbddc->coarse_vec,&lsize);CHKERRQ(ierr);
50       ierr = VecGetArrayRead(pcbddc->coarse_vec,&array_in);CHKERRQ(ierr);
51       ierr = VecGetArray(coarse_nsp_vecs[coarse_nsp_size],&array_out);CHKERRQ(ierr);
52       ierr = PetscMemcpy(array_out,array_in,lsize*sizeof(PetscScalar));CHKERRQ(ierr);
53       ierr = VecRestoreArray(coarse_nsp_vecs[coarse_nsp_size],&array_out);CHKERRQ(ierr);
54       ierr = VecRestoreArrayRead(pcbddc->coarse_vec,&array_in);CHKERRQ(ierr);
55       coarse_nsp_size++;
56     }
57   }
58   for (i=0;i<nsp_size;i++)  {
59     ierr = VecScatterBegin(matis->ctx,nsp_vecs[i],local_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
60     ierr = VecScatterEnd(matis->ctx,nsp_vecs[i],local_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
61     ierr = MatMult(pcbddc->ConstraintMatrix,local_vec,local_primal_vec);CHKERRQ(ierr);
62     ierr = VecScatterBegin(pcbddc->coarse_loc_to_glob,local_primal_vec,pcbddc->coarse_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
63     ierr = VecScatterEnd(pcbddc->coarse_loc_to_glob,local_primal_vec,pcbddc->coarse_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
64     if (coarse_mat) {
65       PetscScalar *array_out;
66       const PetscScalar *array_in;
67       PetscInt lsize;
68       if (pcbddc->dbg_flag) {
69         PetscViewer dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)coarse_mat));
70         ierr = MatMult(coarse_mat,wcoarse_vec,wcoarse_rhs);CHKERRQ(ierr);
71         ierr = VecNorm(wcoarse_rhs,NORM_2,&test_null);CHKERRQ(ierr);
72         ierr = PetscViewerASCIIPrintf(dbg_viewer,"Vec %d coarse null space error % 1.14e\n",i,test_null);CHKERRQ(ierr);
73         ierr = PetscViewerFlush(dbg_viewer);CHKERRQ(ierr);
74       }
75       ierr = VecGetLocalSize(pcbddc->coarse_vec,&lsize);CHKERRQ(ierr);
76       ierr = VecGetArrayRead(pcbddc->coarse_vec,&array_in);CHKERRQ(ierr);
77       ierr = VecGetArray(coarse_nsp_vecs[coarse_nsp_size],&array_out);CHKERRQ(ierr);
78       ierr = PetscMemcpy(array_out,array_in,lsize*sizeof(PetscScalar));CHKERRQ(ierr);
79       ierr = VecRestoreArray(coarse_nsp_vecs[coarse_nsp_size],&array_out);CHKERRQ(ierr);
80       ierr = VecRestoreArrayRead(pcbddc->coarse_vec,&array_in);CHKERRQ(ierr);
81       coarse_nsp_size++;
82     }
83   }
84   if (coarse_nsp_size > 0) {
85     ierr = PCBDDCOrthonormalizeVecs(coarse_nsp_size,coarse_nsp_vecs);CHKERRQ(ierr);
86     ierr = MatNullSpaceCreate(PetscObjectComm((PetscObject)coarse_mat),PETSC_FALSE,coarse_nsp_size,coarse_nsp_vecs,&tempCoarseNullSpace);CHKERRQ(ierr);
87     for (i=0;i<nsp_size+1;i++) {
88       ierr = VecDestroy(&coarse_nsp_vecs[i]);CHKERRQ(ierr);
89     }
90   }
91   if (coarse_mat) {
92     ierr = PetscFree(coarse_nsp_vecs);CHKERRQ(ierr);
93     if (pcbddc->dbg_flag) {
94       ierr = VecDestroy(&wcoarse_vec);CHKERRQ(ierr);
95       ierr = VecDestroy(&wcoarse_rhs);CHKERRQ(ierr);
96     }
97   }
98   ierr = VecDestroy(&local_vec);CHKERRQ(ierr);
99   ierr = VecDestroy(&local_primal_vec);CHKERRQ(ierr);
100   *CoarseNullSpace = tempCoarseNullSpace;
101   PetscFunctionReturn(0);
102 }
103 
104 
105 #undef __FUNCT__
106 #define __FUNCT__ "PCBDDCApplyNullSpaceCorrectionPC"
107 static PetscErrorCode PCBDDCApplyNullSpaceCorrectionPC(PC pc,Vec x,Vec y)
108 {
109   NullSpaceCorrection_ctx pc_ctx;
110   PetscErrorCode          ierr;
111 
112   PetscFunctionBegin;
113   ierr = PCShellGetContext(pc,(void**)&pc_ctx);CHKERRQ(ierr);
114   /* E */
115   ierr = MatMultTranspose(pc_ctx->Lbasis_mat,x,pc_ctx->work_small_2);CHKERRQ(ierr);
116   ierr = MatMultAdd(pc_ctx->Kbasis_mat,pc_ctx->work_small_2,x,pc_ctx->work_full_1);CHKERRQ(ierr);
117   /* P^-1 */
118   ierr = PCApply(pc_ctx->local_pc,pc_ctx->work_full_1,pc_ctx->work_full_2);CHKERRQ(ierr);
119   /* E^T */
120   ierr = MatMultTranspose(pc_ctx->Kbasis_mat,pc_ctx->work_full_2,pc_ctx->work_small_1);CHKERRQ(ierr);
121   ierr = VecScale(pc_ctx->work_small_1,-1.0);CHKERRQ(ierr);
122   ierr = MatMultAdd(pc_ctx->Lbasis_mat,pc_ctx->work_small_1,pc_ctx->work_full_2,pc_ctx->work_full_1);CHKERRQ(ierr);
123   /* Sum contributions */
124   ierr = MatMultAdd(pc_ctx->basis_mat,pc_ctx->work_small_2,pc_ctx->work_full_1,y);CHKERRQ(ierr);
125   PetscFunctionReturn(0);
126 }
127 
128 #undef __FUNCT__
129 #define __FUNCT__ "PCBDDCDestroyNullSpaceCorrectionPC"
130 static PetscErrorCode PCBDDCDestroyNullSpaceCorrectionPC(PC pc)
131 {
132   NullSpaceCorrection_ctx pc_ctx;
133   PetscErrorCode          ierr;
134 
135   PetscFunctionBegin;
136   ierr = PCShellGetContext(pc,(void**)&pc_ctx);CHKERRQ(ierr);
137   ierr = VecDestroy(&pc_ctx->work_small_1);CHKERRQ(ierr);
138   ierr = VecDestroy(&pc_ctx->work_small_2);CHKERRQ(ierr);
139   ierr = VecDestroy(&pc_ctx->work_full_1);CHKERRQ(ierr);
140   ierr = VecDestroy(&pc_ctx->work_full_2);CHKERRQ(ierr);
141   ierr = MatDestroy(&pc_ctx->basis_mat);CHKERRQ(ierr);
142   ierr = MatDestroy(&pc_ctx->Lbasis_mat);CHKERRQ(ierr);
143   ierr = MatDestroy(&pc_ctx->Kbasis_mat);CHKERRQ(ierr);
144   ierr = PCDestroy(&pc_ctx->local_pc);CHKERRQ(ierr);
145   ierr = PetscFree(pc_ctx);CHKERRQ(ierr);
146   PetscFunctionReturn(0);
147 }
148 
149 /*
150 PETSC_EXTERN PetscErrorCode PCBDDCApplyNullSpaceCorrectionPC(PC,Vec,Vec);
151 PETSC_EXTERN PetscErrorCode PCBDDCDestroyNullSpaceCorrectionPC(PC);
152 */
153 
154 #undef __FUNCT__
155 #define __FUNCT__ "PCBDDCNullSpaceAssembleCorrection"
156 PetscErrorCode PCBDDCNullSpaceAssembleCorrection(PC pc,IS local_dofs)
157 {
158   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
159   PC_IS          *pcis = (PC_IS*)pc->data;
160   Mat_IS*        matis = (Mat_IS*)pc->pmat->data;
161   KSP            *local_ksp;
162   PC             newpc;
163   NullSpaceCorrection_ctx  shell_ctx;
164   Mat            local_mat,local_pmat,small_mat,inv_small_mat;
165   MatStructure   local_mat_struct;
166   Vec            work1,work2;
167   const Vec      *nullvecs;
168   VecScatter     scatter_ctx;
169   IS             is_aux;
170   MatFactorInfo  matinfo;
171   PetscScalar    *basis_mat,*Kbasis_mat,*array,*array_mat;
172   PetscScalar    one = 1.0,zero = 0.0, m_one = -1.0;
173   PetscInt       basis_dofs,basis_size,nnsp_size,i,k,n_I,n_R;
174   PetscBool      nnsp_has_cnst;
175   PetscErrorCode ierr;
176 
177   PetscFunctionBegin;
178   /* Infer the local solver */
179   ierr = ISGetSize(local_dofs,&basis_dofs);CHKERRQ(ierr);
180   ierr = VecGetSize(pcis->vec1_D,&n_I);CHKERRQ(ierr);
181   ierr = VecGetSize(pcbddc->vec1_R,&n_R);CHKERRQ(ierr);
182   if (basis_dofs == n_I) {
183     /* Dirichlet solver */
184     local_ksp = &pcbddc->ksp_D;
185   } else if (basis_dofs == n_R) {
186     /* Neumann solver */
187     local_ksp = &pcbddc->ksp_R;
188   } else {
189     SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in %s: unknown local IS size %d. n_I=%d, n_R=%d)\n",__FUNCT__,basis_dofs,n_I,n_R);
190   }
191   ierr = KSPGetOperators(*local_ksp,&local_mat,&local_pmat,&local_mat_struct);CHKERRQ(ierr);
192 
193   /* Get null space vecs */
194   ierr = MatNullSpaceGetVecs(pcbddc->NullSpace,&nnsp_has_cnst,&nnsp_size,&nullvecs);CHKERRQ(ierr);
195   basis_size = nnsp_size;
196   if (nnsp_has_cnst) {
197     basis_size++;
198   }
199 
200   if (basis_dofs) {
201      /* Create shell ctx */
202      ierr = PetscMalloc(sizeof(*shell_ctx),&shell_ctx);CHKERRQ(ierr);
203 
204      /* Create work vectors in shell context */
205      ierr = VecCreate(PETSC_COMM_SELF,&shell_ctx->work_small_1);CHKERRQ(ierr);
206      ierr = VecSetSizes(shell_ctx->work_small_1,basis_size,basis_size);CHKERRQ(ierr);
207      ierr = VecSetType(shell_ctx->work_small_1,VECSEQ);CHKERRQ(ierr);
208      ierr = VecDuplicate(shell_ctx->work_small_1,&shell_ctx->work_small_2);CHKERRQ(ierr);
209      ierr = VecCreate(PETSC_COMM_SELF,&shell_ctx->work_full_1);CHKERRQ(ierr);
210      ierr = VecSetSizes(shell_ctx->work_full_1,basis_dofs,basis_dofs);CHKERRQ(ierr);
211      ierr = VecSetType(shell_ctx->work_full_1,VECSEQ);CHKERRQ(ierr);
212      ierr = VecDuplicate(shell_ctx->work_full_1,&shell_ctx->work_full_2);CHKERRQ(ierr);
213 
214      /* Allocate workspace */
215      ierr = MatCreateSeqDense(PETSC_COMM_SELF,basis_dofs,basis_size,NULL,&shell_ctx->basis_mat );CHKERRQ(ierr);
216      ierr = MatCreateSeqDense(PETSC_COMM_SELF,basis_dofs,basis_size,NULL,&shell_ctx->Kbasis_mat);CHKERRQ(ierr);
217      ierr = MatDenseGetArray(shell_ctx->basis_mat,&basis_mat);CHKERRQ(ierr);
218      ierr = MatDenseGetArray(shell_ctx->Kbasis_mat,&Kbasis_mat);CHKERRQ(ierr);
219 
220      /* Restrict local null space on selected dofs (Dirichlet or Neumann)
221         and compute matrices N and K*N */
222      ierr = VecDuplicate(shell_ctx->work_full_1,&work1);CHKERRQ(ierr);
223      ierr = VecDuplicate(shell_ctx->work_full_1,&work2);CHKERRQ(ierr);
224      ierr = VecScatterCreate(pcis->vec1_N,local_dofs,work1,(IS)0,&scatter_ctx);CHKERRQ(ierr);
225   }
226 
227   for (k=0;k<nnsp_size;k++) {
228     ierr = VecScatterBegin(matis->ctx,nullvecs[k],pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
229     ierr = VecScatterEnd(matis->ctx,nullvecs[k],pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
230     if (basis_dofs) {
231       ierr = VecPlaceArray(work1,(const PetscScalar*)&basis_mat[k*basis_dofs]);CHKERRQ(ierr);
232       ierr = VecScatterBegin(scatter_ctx,pcis->vec1_N,work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
233       ierr = VecScatterEnd(scatter_ctx,pcis->vec1_N,work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
234       ierr = VecPlaceArray(work2,(const PetscScalar*)&Kbasis_mat[k*basis_dofs]);CHKERRQ(ierr);
235       ierr = MatMult(local_mat,work1,work2);CHKERRQ(ierr);
236       ierr = VecResetArray(work1);CHKERRQ(ierr);
237       ierr = VecResetArray(work2);CHKERRQ(ierr);
238     }
239   }
240 
241   if (basis_dofs) {
242     if (nnsp_has_cnst) {
243       ierr = VecPlaceArray(work1,(const PetscScalar*)&basis_mat[k*basis_dofs]);CHKERRQ(ierr);
244       ierr = VecSet(work1,one);CHKERRQ(ierr);
245       ierr = VecPlaceArray(work2,(const PetscScalar*)&Kbasis_mat[k*basis_dofs]);CHKERRQ(ierr);
246       ierr = MatMult(local_mat,work1,work2);CHKERRQ(ierr);
247       ierr = VecResetArray(work1);CHKERRQ(ierr);
248       ierr = VecResetArray(work2);CHKERRQ(ierr);
249     }
250     ierr = VecDestroy(&work1);CHKERRQ(ierr);
251     ierr = VecDestroy(&work2);CHKERRQ(ierr);
252     ierr = VecScatterDestroy(&scatter_ctx);CHKERRQ(ierr);
253     ierr = MatDenseRestoreArray(shell_ctx->basis_mat,&basis_mat);CHKERRQ(ierr);
254     ierr = MatDenseRestoreArray(shell_ctx->Kbasis_mat,&Kbasis_mat);CHKERRQ(ierr);
255 
256     /* Assemble another Mat object in shell context */
257     ierr = MatTransposeMatMult(shell_ctx->basis_mat,shell_ctx->Kbasis_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&small_mat);CHKERRQ(ierr);
258     ierr = MatFactorInfoInitialize(&matinfo);CHKERRQ(ierr);
259     ierr = ISCreateStride(PETSC_COMM_SELF,basis_size,0,1,&is_aux);CHKERRQ(ierr);
260     ierr = MatLUFactor(small_mat,is_aux,is_aux,&matinfo);CHKERRQ(ierr);
261     ierr = ISDestroy(&is_aux);CHKERRQ(ierr);
262     ierr = PetscMalloc(basis_size*basis_size*sizeof(PetscScalar),&array_mat);CHKERRQ(ierr);
263     for (k=0;k<basis_size;k++) {
264       ierr = VecSet(shell_ctx->work_small_1,zero);CHKERRQ(ierr);
265       ierr = VecSetValue(shell_ctx->work_small_1,k,one,INSERT_VALUES);CHKERRQ(ierr);
266       ierr = VecAssemblyBegin(shell_ctx->work_small_1);CHKERRQ(ierr);
267       ierr = VecAssemblyEnd(shell_ctx->work_small_1);CHKERRQ(ierr);
268       ierr = MatSolve(small_mat,shell_ctx->work_small_1,shell_ctx->work_small_2);CHKERRQ(ierr);
269       ierr = VecGetArrayRead(shell_ctx->work_small_2,(const PetscScalar**)&array);CHKERRQ(ierr);
270       for (i=0;i<basis_size;i++) {
271         array_mat[i*basis_size+k]=array[i];
272       }
273       ierr = VecRestoreArrayRead(shell_ctx->work_small_2,(const PetscScalar**)&array);CHKERRQ(ierr);
274     }
275     ierr = MatCreateSeqDense(PETSC_COMM_SELF,basis_size,basis_size,array_mat,&inv_small_mat);CHKERRQ(ierr);
276     ierr = MatMatMult(shell_ctx->basis_mat,inv_small_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&shell_ctx->Lbasis_mat);CHKERRQ(ierr);
277     ierr = PetscFree(array_mat);CHKERRQ(ierr);
278     ierr = MatDestroy(&inv_small_mat);CHKERRQ(ierr);
279     ierr = MatDestroy(&small_mat);CHKERRQ(ierr);
280     ierr = MatScale(shell_ctx->Kbasis_mat,m_one);CHKERRQ(ierr);
281 
282     /* Rebuild local PC */
283     ierr = KSPGetPC(*local_ksp,&shell_ctx->local_pc);CHKERRQ(ierr);
284     ierr = PetscObjectReference((PetscObject)shell_ctx->local_pc);CHKERRQ(ierr);
285     ierr = PCCreate(PETSC_COMM_SELF,&newpc);CHKERRQ(ierr);
286     ierr = PCSetOperators(newpc,local_mat,local_mat,SAME_PRECONDITIONER);CHKERRQ(ierr);
287     ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr);
288     ierr = PCShellSetContext(newpc,shell_ctx);CHKERRQ(ierr);
289     ierr = PCShellSetApply(newpc,PCBDDCApplyNullSpaceCorrectionPC);CHKERRQ(ierr);
290     ierr = PCShellSetDestroy(newpc,PCBDDCDestroyNullSpaceCorrectionPC);CHKERRQ(ierr);
291     ierr = PCSetUp(newpc);CHKERRQ(ierr);
292     ierr = KSPSetPC(*local_ksp,newpc);CHKERRQ(ierr);
293     ierr = PCDestroy(&newpc);CHKERRQ(ierr);
294     ierr = KSPSetUp(*local_ksp);CHKERRQ(ierr);
295   }
296   /* test */
297   if (pcbddc->dbg_flag && basis_dofs) {
298     KSP         check_ksp;
299     PC          check_pc;
300     Mat         test_mat;
301     Vec         work3;
302     PetscReal   test_err,lambda_min,lambda_max;
303     PetscBool   setsym,issym=PETSC_FALSE;
304     PetscInt    tabs;
305 
306     ierr = PetscViewerASCIIGetTab(pcbddc->dbg_viewer,&tabs);CHKERRQ(ierr);
307     ierr = KSPGetPC(*local_ksp,&check_pc);CHKERRQ(ierr);
308     ierr = VecDuplicate(shell_ctx->work_full_1,&work1);CHKERRQ(ierr);
309     ierr = VecDuplicate(shell_ctx->work_full_1,&work2);CHKERRQ(ierr);
310     ierr = VecDuplicate(shell_ctx->work_full_1,&work3);CHKERRQ(ierr);
311     ierr = VecSetRandom(shell_ctx->work_small_1,NULL);CHKERRQ(ierr);
312     ierr = MatMult(shell_ctx->basis_mat,shell_ctx->work_small_1,work1);CHKERRQ(ierr);
313     ierr = VecCopy(work1,work2);CHKERRQ(ierr);
314     ierr = MatMult(local_mat,work1,work3);CHKERRQ(ierr);
315     ierr = PCApply(check_pc,work3,work1);CHKERRQ(ierr);
316     ierr = VecAXPY(work1,m_one,work2);CHKERRQ(ierr);
317     ierr = VecNorm(work1,NORM_INFINITY,&test_err);CHKERRQ(ierr);
318     ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d error for nullspace correction for ",PetscGlobalRank);
319     ierr = PetscViewerASCIIUseTabs(pcbddc->dbg_viewer,PETSC_FALSE);CHKERRQ(ierr);
320     if (basis_dofs == n_I) {
321       ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Dirichlet ");
322     } else {
323       ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Neumann ");
324     }
325     ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"solver is :%1.14e\n",test_err);
326     ierr = PetscViewerASCIISetTab(pcbddc->dbg_viewer,tabs);CHKERRQ(ierr);
327     ierr = PetscViewerASCIIUseTabs(pcbddc->dbg_viewer,PETSC_TRUE);CHKERRQ(ierr);
328 
329     ierr = MatTransposeMatMult(shell_ctx->Lbasis_mat,shell_ctx->Kbasis_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&test_mat);CHKERRQ(ierr);
330     ierr = MatShift(test_mat,one);CHKERRQ(ierr);
331     ierr = MatNorm(test_mat,NORM_INFINITY,&test_err);CHKERRQ(ierr);
332     ierr = MatDestroy(&test_mat);CHKERRQ(ierr);
333     ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d error for nullspace matrices is :%1.14e\n",PetscGlobalRank,test_err);
334 
335     /* Create ksp object suitable for extreme eigenvalues' estimation */
336     ierr = KSPCreate(PETSC_COMM_SELF,&check_ksp);CHKERRQ(ierr);
337     ierr = KSPSetOperators(check_ksp,local_mat,local_mat,SAME_PRECONDITIONER);CHKERRQ(ierr);
338     ierr = KSPSetTolerances(check_ksp,1.e-8,1.e-8,PETSC_DEFAULT,basis_dofs);CHKERRQ(ierr);
339     ierr = KSPSetComputeSingularValues(check_ksp,PETSC_TRUE);CHKERRQ(ierr);
340     ierr = MatIsSymmetricKnown(pc->pmat,&setsym,&issym);CHKERRQ(ierr);
341     if (issym) {
342       ierr = KSPSetType(check_ksp,KSPCG);CHKERRQ(ierr);
343     }
344     ierr = KSPSetPC(check_ksp,check_pc);CHKERRQ(ierr);
345     ierr = KSPSetUp(check_ksp);CHKERRQ(ierr);
346     ierr = VecSetRandom(work1,NULL);CHKERRQ(ierr);
347     ierr = MatMult(local_mat,work1,work2);CHKERRQ(ierr);
348     ierr = KSPSolve(check_ksp,work2,work2);CHKERRQ(ierr);
349     ierr = VecAXPY(work2,m_one,work1);CHKERRQ(ierr);
350     ierr = VecNorm(work2,NORM_INFINITY,&test_err);CHKERRQ(ierr);
351     ierr = KSPComputeExtremeSingularValues(check_ksp,&lambda_max,&lambda_min);CHKERRQ(ierr);
352     ierr = KSPGetIterationNumber(check_ksp,&k);CHKERRQ(ierr);
353     ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d error for adapted KSP %1.14e (it %d, eigs %1.6e %1.6e)\n",PetscGlobalRank,test_err,k,lambda_min,lambda_max);
354     ierr = KSPDestroy(&check_ksp);CHKERRQ(ierr);
355     ierr = VecDestroy(&work1);CHKERRQ(ierr);
356     ierr = VecDestroy(&work2);CHKERRQ(ierr);
357     ierr = VecDestroy(&work3);CHKERRQ(ierr);
358   }
359   /* all processes shoud call this, even the void ones */
360   if (pcbddc->dbg_flag) {
361     ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr);
362   }
363   PetscFunctionReturn(0);
364 }
365 
366 /* uncomment to test nullspace adaptation when change of basis has been requested */
367 /* #define PCBDDC_TESTNULLSPACE 1 */
368 #undef __FUNCT__
369 #define __FUNCT__ "PCBDDCNullSpaceAdaptGlobal"
370 PetscErrorCode PCBDDCNullSpaceAdaptGlobal(PC pc)
371 {
372   PC_IS*         pcis = (PC_IS*)(pc->data);
373   PC_BDDC*       pcbddc = (PC_BDDC*)(pc->data);
374   KSP            inv_change;
375   PC             pc_change;
376   const Vec      *nsp_vecs;
377   Vec            *new_nsp_vecs;
378   PetscInt       i,nsp_size,new_nsp_size,start_new;
379   PetscBool      nsp_has_cnst;
380   MatNullSpace   new_nsp;
381   PetscErrorCode ierr;
382 
383   PetscFunctionBegin;
384   /* create KSP for change of basis */
385   ierr = KSPCreate(PETSC_COMM_SELF,&inv_change);CHKERRQ(ierr);
386   ierr = KSPSetOperators(inv_change,pcbddc->ChangeOfBasisMatrix,pcbddc->ChangeOfBasisMatrix,SAME_PRECONDITIONER);CHKERRQ(ierr);
387   ierr = KSPSetType(inv_change,KSPPREONLY);CHKERRQ(ierr);
388   ierr = KSPGetPC(inv_change,&pc_change);CHKERRQ(ierr);
389   ierr = PCSetType(pc_change,PCLU);CHKERRQ(ierr);
390   ierr = KSPSetUp(inv_change);CHKERRQ(ierr);
391   /* get nullspace and transform it */
392   ierr = MatNullSpaceGetVecs(pcbddc->NullSpace,&nsp_has_cnst,&nsp_size,&nsp_vecs);CHKERRQ(ierr);
393   new_nsp_size = nsp_size;
394   if (nsp_has_cnst) {
395     new_nsp_size++;
396   }
397   ierr = PetscMalloc(new_nsp_size*sizeof(Vec),&new_nsp_vecs);CHKERRQ(ierr);
398   for (i=0;i<new_nsp_size;i++) {
399     ierr = VecDuplicate(pcis->vec1_global,&new_nsp_vecs[i]);CHKERRQ(ierr);
400   }
401   start_new = 0;
402   if (nsp_has_cnst) {
403     start_new = 1;
404     ierr = VecSet(new_nsp_vecs[0],1.0);CHKERRQ(ierr);
405     ierr = VecSet(pcis->vec1_B,1.0);CHKERRQ(ierr);
406     ierr = KSPSolve(inv_change,pcis->vec1_B,pcis->vec1_B);CHKERRQ(ierr);
407     ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,new_nsp_vecs[0],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
408     ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,new_nsp_vecs[0],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
409   }
410   for (i=0;i<nsp_size;i++) {
411     ierr = VecCopy(nsp_vecs[i],new_nsp_vecs[i+start_new]);CHKERRQ(ierr);
412     ierr = VecScatterBegin(pcis->global_to_B,nsp_vecs[i],pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
413     ierr = VecScatterEnd(pcis->global_to_B,nsp_vecs[i],pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
414     ierr = KSPSolve(inv_change,pcis->vec1_B,pcis->vec1_B);
415     ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,new_nsp_vecs[i+start_new],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
416     ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,new_nsp_vecs[i+start_new],INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
417   }
418   ierr = PCBDDCOrthonormalizeVecs(new_nsp_size,new_nsp_vecs);CHKERRQ(ierr);
419 #ifdef PCBDDC_TESTNULLSPACE
420   {
421     PetscBool nsp_t=PETSC_FALSE;
422     PetscReal test_norm;
423     Mat       temp_mat;
424     Mat_IS*   matis = (Mat_IS*)pc->pmat->data;
425     ierr = PetscPrintf(PetscObjectComm((PetscObject)(pc->pmat)),"Testing BDDC nullspace (mat nullspace)\n",nsp_t);
426     ierr = MatNullSpaceTest(pcbddc->NullSpace,pc->pmat,&nsp_t);CHKERRQ(ierr);
427     ierr = PetscPrintf(PetscObjectComm((PetscObject)(pc->pmat)),"ORI ORI: %d\n",nsp_t);
428     temp_mat = matis->A;
429     matis->A = pcbddc->local_mat;
430     pcbddc->local_mat = temp_mat;
431     ierr = MatNullSpaceTest(pcbddc->NullSpace,pc->pmat,&nsp_t);CHKERRQ(ierr);
432     ierr = PetscPrintf(PetscObjectComm((PetscObject)(pc->pmat)),"NEW ORI: %d\n",nsp_t);
433     for (i=0;i<new_nsp_size;i++) {
434       ierr = MatMult(pc->pmat,new_nsp_vecs[i],pcis->vec1_global);CHKERRQ(ierr);
435       ierr = VecNorm(pcis->vec1_global,NORM_2,&test_norm);CHKERRQ(ierr);
436       if (test_norm > 1.e-12) {
437         ierr = PetscPrintf(PetscObjectComm((PetscObject)(pc->pmat)),"------------ERROR VEC %d------------------\n",i);
438         ierr = VecView(pcis->vec1_global,PETSC_VIEWER_STDOUT_WORLD);
439         ierr = PetscPrintf(PetscObjectComm((PetscObject)(pc->pmat)),"------------------------------------------\n");
440       }
441     }
442   }
443 #endif
444 
445   ierr = KSPDestroy(&inv_change);CHKERRQ(ierr);
446   ierr = MatNullSpaceCreate(PetscObjectComm((PetscObject)pc),PETSC_FALSE,new_nsp_size,new_nsp_vecs,&new_nsp);CHKERRQ(ierr);
447   ierr = PCBDDCSetNullSpace(pc,new_nsp);CHKERRQ(ierr);
448   ierr = MatNullSpaceDestroy(&new_nsp);CHKERRQ(ierr);
449 #ifdef PCBDDC_TESTNULLSPACE
450   {
451     PetscBool nsp_t=PETSC_FALSE;
452     Mat       temp_mat;
453     Mat_IS*   matis = (Mat_IS*)pc->pmat->data;
454     ierr = MatNullSpaceTest(pcbddc->NullSpace,pc->pmat,&nsp_t);CHKERRQ(ierr);
455     ierr = PetscPrintf(PetscObjectComm((PetscObject)(pc->pmat)),"NEW NEW: %d\n",nsp_t);
456     temp_mat = matis->A;
457     matis->A = pcbddc->local_mat;
458     pcbddc->local_mat = temp_mat;
459     ierr = MatNullSpaceTest(pcbddc->NullSpace,pc->pmat,&nsp_t);CHKERRQ(ierr);
460     ierr = PetscPrintf(PetscObjectComm((PetscObject)(pc->pmat)),"ORI NEW: %d\n",nsp_t);
461   }
462 #endif
463 
464   for (i=0;i<new_nsp_size;i++) {
465     ierr = VecDestroy(&new_nsp_vecs[i]);CHKERRQ(ierr);
466   }
467   ierr = PetscFree(new_nsp_vecs);CHKERRQ(ierr);
468   PetscFunctionReturn(0);
469 }
470