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