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