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