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