Actual source code: bddcnullspace.c
1: #include <../src/ksp/pc/impls/bddc/bddc.h>
2: #include <../src/ksp/pc/impls/bddc/bddcprivate.h>
3: #include <../src/mat/impls/dense/seq/dense.h>
5: /* E + small_solve */
6: static PetscErrorCode PCBDDCNullSpaceCorrPreSolve(KSP ksp,Vec y,Vec x, void* ctx)
7: {
8: NullSpaceCorrection_ctx corr_ctx = (NullSpaceCorrection_ctx)ctx;
9: Mat K;
11: PetscLogEventBegin(corr_ctx->evapply,ksp,0,0,0);
12: MatMultTranspose(corr_ctx->basis_mat,y,corr_ctx->sw[0]);
13: if (corr_ctx->symm) {
14: MatMult(corr_ctx->inv_smat,corr_ctx->sw[0],corr_ctx->sw[1]);
15: } else {
16: MatMultTranspose(corr_ctx->inv_smat,corr_ctx->sw[0],corr_ctx->sw[1]);
17: }
18: VecScale(corr_ctx->sw[1],-1.0);
19: MatMult(corr_ctx->basis_mat,corr_ctx->sw[1],corr_ctx->fw[0]);
20: VecScale(corr_ctx->sw[1],-1.0);
21: KSPGetOperators(ksp,&K,NULL);
22: MatMultAdd(K,corr_ctx->fw[0],y,y);
23: PetscLogEventEnd(corr_ctx->evapply,ksp,0,0,0);
24: return 0;
25: }
27: /* E^t + small */
28: static PetscErrorCode PCBDDCNullSpaceCorrPostSolve(KSP ksp,Vec y,Vec x, void* ctx)
29: {
30: NullSpaceCorrection_ctx corr_ctx = (NullSpaceCorrection_ctx)ctx;
31: Mat K;
33: PetscLogEventBegin(corr_ctx->evapply,ksp,0,0,0);
34: KSPGetOperators(ksp,&K,NULL);
35: if (corr_ctx->symm) {
36: MatMult(K,x,corr_ctx->fw[0]);
37: } else {
38: MatMultTranspose(K,x,corr_ctx->fw[0]);
39: }
40: MatMultTranspose(corr_ctx->basis_mat,corr_ctx->fw[0],corr_ctx->sw[0]);
41: VecScale(corr_ctx->sw[0],-1.0);
42: MatMult(corr_ctx->inv_smat,corr_ctx->sw[0],corr_ctx->sw[2]);
43: MatMultAdd(corr_ctx->basis_mat,corr_ctx->sw[2],x,corr_ctx->fw[0]);
44: VecScale(corr_ctx->fw[0],corr_ctx->scale);
45: /* Sum contributions from approximate solver and projected system */
46: MatMultAdd(corr_ctx->basis_mat,corr_ctx->sw[1],corr_ctx->fw[0],x);
47: PetscLogEventEnd(corr_ctx->evapply,ksp,0,0,0);
48: return 0;
49: }
51: static PetscErrorCode PCBDDCNullSpaceCorrDestroy(void * ctx)
52: {
53: NullSpaceCorrection_ctx corr_ctx = (NullSpaceCorrection_ctx)ctx;
55: VecDestroyVecs(3,&corr_ctx->sw);
56: VecDestroyVecs(1,&corr_ctx->fw);
57: MatDestroy(&corr_ctx->basis_mat);
58: MatDestroy(&corr_ctx->inv_smat);
59: PetscFree(corr_ctx);
60: return 0;
61: }
63: PetscErrorCode PCBDDCNullSpaceAssembleCorrection(PC pc, PetscBool isdir, PetscBool needscaling)
64: {
65: PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
66: MatNullSpace NullSpace = NULL;
67: KSP local_ksp;
68: NullSpaceCorrection_ctx shell_ctx;
69: Mat local_mat,local_pmat,dmat,Kbasis_mat;
70: Vec v;
71: PetscContainer c;
72: PetscInt basis_size;
73: IS zerorows;
74: PetscBool iscusp;
76: if (isdir) local_ksp = pcbddc->ksp_D; /* Dirichlet solver */
77: else local_ksp = pcbddc->ksp_R; /* Neumann solver */
78: KSPGetOperators(local_ksp,&local_mat,&local_pmat);
79: MatGetNearNullSpace(local_pmat,&NullSpace);
80: if (!NullSpace) {
81: if (pcbddc->dbg_flag) {
82: PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d doesn't have local (near) nullspace: no need for correction in %s solver \n",PetscGlobalRank,isdir ? "Dirichlet" : "Neumann");
83: }
84: return 0;
85: }
86: PetscObjectQuery((PetscObject)NullSpace,"_PBDDC_Null_dmat",(PetscObject*)&dmat);
88: PetscLogEventBegin(PC_BDDC_ApproxSetUp[pcbddc->current_level],pc,0,0,0);
90: PetscNew(&shell_ctx);
91: shell_ctx->scale = 1.0;
92: PetscObjectReference((PetscObject)dmat);
93: shell_ctx->basis_mat = dmat;
94: MatGetSize(dmat,NULL,&basis_size);
95: shell_ctx->evapply = PC_BDDC_ApproxApply[pcbddc->current_level];
97: MatGetOption(local_mat,MAT_SYMMETRIC,&shell_ctx->symm);
99: /* explicit construct (Phi^T K Phi)^-1 */
100: PetscObjectTypeCompare((PetscObject)local_mat,MATSEQAIJCUSPARSE,&iscusp);
101: if (iscusp) {
102: MatConvert(shell_ctx->basis_mat,MATSEQDENSECUDA,MAT_INPLACE_MATRIX,&shell_ctx->basis_mat);
103: }
104: MatMatMult(local_mat,shell_ctx->basis_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&Kbasis_mat);
105: MatTransposeMatMult(Kbasis_mat,shell_ctx->basis_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&shell_ctx->inv_smat);
106: MatDestroy(&Kbasis_mat);
107: MatBindToCPU(shell_ctx->inv_smat,PETSC_TRUE);
108: MatFindZeroRows(shell_ctx->inv_smat,&zerorows);
109: if (zerorows) { /* linearly dependent basis */
110: const PetscInt *idxs;
111: PetscInt i,nz;
113: ISGetLocalSize(zerorows,&nz);
114: ISGetIndices(zerorows,&idxs);
115: for (i=0;i<nz;i++) {
116: MatSetValue(shell_ctx->inv_smat,idxs[i],idxs[i],1.0,INSERT_VALUES);
117: }
118: ISRestoreIndices(zerorows,&idxs);
119: MatAssemblyBegin(shell_ctx->inv_smat,MAT_FINAL_ASSEMBLY);
120: MatAssemblyEnd(shell_ctx->inv_smat,MAT_FINAL_ASSEMBLY);
121: }
122: MatLUFactor(shell_ctx->inv_smat,NULL,NULL,NULL);
123: MatSeqDenseInvertFactors_Private(shell_ctx->inv_smat);
124: if (zerorows) { /* linearly dependent basis */
125: const PetscInt *idxs;
126: PetscInt i,nz;
128: ISGetLocalSize(zerorows,&nz);
129: ISGetIndices(zerorows,&idxs);
130: for (i=0;i<nz;i++) {
131: MatSetValue(shell_ctx->inv_smat,idxs[i],idxs[i],0.0,INSERT_VALUES);
132: }
133: ISRestoreIndices(zerorows,&idxs);
134: MatAssemblyBegin(shell_ctx->inv_smat,MAT_FINAL_ASSEMBLY);
135: MatAssemblyEnd(shell_ctx->inv_smat,MAT_FINAL_ASSEMBLY);
136: }
137: ISDestroy(&zerorows);
139: /* Create work vectors in shell context */
140: MatCreateVecs(shell_ctx->inv_smat,&v,NULL);
141: KSPCreateVecs(local_ksp,1,&shell_ctx->fw,0,NULL);
142: VecDuplicateVecs(v,3,&shell_ctx->sw);
143: VecDestroy(&v);
145: /* add special pre/post solve to KSP (see [1], eq. 48) */
146: KSPSetPreSolve(local_ksp,PCBDDCNullSpaceCorrPreSolve,shell_ctx);
147: KSPSetPostSolve(local_ksp,PCBDDCNullSpaceCorrPostSolve,shell_ctx);
148: PetscContainerCreate(PetscObjectComm((PetscObject)local_ksp),&c);
149: PetscContainerSetPointer(c,shell_ctx);
150: PetscContainerSetUserDestroy(c,PCBDDCNullSpaceCorrDestroy);
151: PetscObjectCompose((PetscObject)local_ksp,"_PCBDDC_Null_PrePost_ctx",(PetscObject)c);
152: PetscContainerDestroy(&c);
154: /* Create ksp object suitable for extreme eigenvalues' estimation */
155: if (needscaling || pcbddc->dbg_flag) {
156: KSP check_ksp;
157: PC local_pc;
158: Vec work1,work2;
159: const char* prefix;
160: PetscReal test_err,lambda_min,lambda_max;
161: PetscInt k,maxit;
163: VecDuplicate(shell_ctx->fw[0],&work1);
164: VecDuplicate(shell_ctx->fw[0],&work2);
165: KSPCreate(PETSC_COMM_SELF,&check_ksp);
166: if (local_mat->spd) {
167: KSPSetType(check_ksp,KSPCG);
168: }
169: PetscObjectIncrementTabLevel((PetscObject)check_ksp,(PetscObject)local_ksp,0);
170: KSPGetOptionsPrefix(local_ksp,&prefix);
171: KSPSetOptionsPrefix(check_ksp,prefix);
172: KSPAppendOptionsPrefix(check_ksp,"approximate_scale_");
173: KSPSetErrorIfNotConverged(check_ksp,PETSC_FALSE);
174: KSPSetOperators(check_ksp,local_mat,local_pmat);
175: KSPSetComputeSingularValues(check_ksp,PETSC_TRUE);
176: KSPSetPreSolve(check_ksp,PCBDDCNullSpaceCorrPreSolve,shell_ctx);
177: KSPSetPostSolve(check_ksp,PCBDDCNullSpaceCorrPostSolve,shell_ctx);
178: KSPSetTolerances(check_ksp,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,PETSC_DEFAULT);
179: KSPSetFromOptions(check_ksp);
180: /* setup with default maxit, then set maxit to min(10,any_set_from_command_line) (bug in computing eigenvalues when chaning the number of iterations */
181: KSPSetUp(check_ksp);
182: KSPGetPC(local_ksp,&local_pc);
183: KSPSetPC(check_ksp,local_pc);
184: KSPGetTolerances(check_ksp,NULL,NULL,NULL,&maxit);
185: KSPSetTolerances(check_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,PetscMin(10,maxit));
186: VecSetRandom(work2,NULL);
187: MatMult(local_mat,work2,work1);
188: KSPSolve(check_ksp,work1,work1);
189: KSPCheckSolve(check_ksp,pc,work1);
190: VecAXPY(work1,-1.,work2);
191: VecNorm(work1,NORM_INFINITY,&test_err);
192: KSPComputeExtremeSingularValues(check_ksp,&lambda_max,&lambda_min);
193: KSPGetIterationNumber(check_ksp,&k);
194: if (pcbddc->dbg_flag) {
195: if (isdir) {
196: PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Dirichlet adapted solver (no scale) %1.14e (it %D, eigs %1.6e %1.6e)\n",PetscGlobalRank,test_err,k,lambda_min,lambda_max);
197: } else {
198: PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Neumann adapted solver (no scale) %1.14e (it %D, eigs %1.6e %1.6e)\n",PetscGlobalRank,test_err,k,lambda_min,lambda_max);
199: }
200: }
201: if (needscaling) shell_ctx->scale = 1.0/lambda_max;
203: if (needscaling && pcbddc->dbg_flag) { /* test for scaling factor */
204: PC new_pc;
206: VecSetRandom(work2,NULL);
207: MatMult(local_mat,work2,work1);
208: PCCreate(PetscObjectComm((PetscObject)check_ksp),&new_pc);
209: PCSetType(new_pc,PCKSP);
210: PCSetOperators(new_pc,local_mat,local_pmat);
211: PCKSPSetKSP(new_pc,local_ksp);
212: KSPSetPC(check_ksp,new_pc);
213: PCDestroy(&new_pc);
214: KSPSetTolerances(check_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,maxit);
215: KSPSetPreSolve(check_ksp,NULL,NULL);
216: KSPSetPostSolve(check_ksp,NULL,NULL);
217: KSPSolve(check_ksp,work1,work1);
218: KSPCheckSolve(check_ksp,pc,work1);
219: VecAXPY(work1,-1.,work2);
220: VecNorm(work1,NORM_INFINITY,&test_err);
221: KSPComputeExtremeSingularValues(check_ksp,&lambda_max,&lambda_min);
222: KSPGetIterationNumber(check_ksp,&k);
223: if (pcbddc->dbg_flag) {
224: if (isdir) {
225: PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Dirichlet adapted solver (scale %g) %1.14e (it %D, eigs %1.6e %1.6e)\n",PetscGlobalRank,(double)PetscRealPart(shell_ctx->scale),test_err,k,lambda_min,lambda_max);
226: } else {
227: PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Neumann adapted solver (scale %g) %1.14e (it %D, eigs %1.6e %1.6e)\n",PetscGlobalRank,(double)PetscRealPart(shell_ctx->scale),test_err,k,lambda_min,lambda_max);
228: }
229: }
230: }
231: KSPDestroy(&check_ksp);
232: VecDestroy(&work1);
233: VecDestroy(&work2);
234: }
235: PetscLogEventEnd(PC_BDDC_ApproxSetUp[pcbddc->current_level],pc,0,0,0);
237: if (pcbddc->dbg_flag) {
238: Vec work1,work2,work3;
239: PetscReal test_err;
241: /* check nullspace basis is solved exactly */
242: VecDuplicate(shell_ctx->fw[0],&work1);
243: VecDuplicate(shell_ctx->fw[0],&work2);
244: VecDuplicate(shell_ctx->fw[0],&work3);
245: VecSetRandom(shell_ctx->sw[0],NULL);
246: MatMult(shell_ctx->basis_mat,shell_ctx->sw[0],work1);
247: VecCopy(work1,work2);
248: MatMult(local_mat,work1,work3);
249: KSPSolve(local_ksp,work3,work1);
250: VecAXPY(work1,-1.,work2);
251: VecNorm(work1,NORM_INFINITY,&test_err);
252: if (isdir) {
253: PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Dirichlet nullspace correction solver: %1.14e\n",PetscGlobalRank,test_err);
254: } else {
255: PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Neumann nullspace correction solver: %1.14e\n",PetscGlobalRank,test_err);
256: }
257: VecDestroy(&work1);
258: VecDestroy(&work2);
259: VecDestroy(&work3);
260: }
261: return 0;
262: }