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 CHKERRQ(PetscLogEventBegin(corr_ctx->evapply,ksp,0,0,0)); 13 CHKERRQ(MatMultTranspose(corr_ctx->basis_mat,y,corr_ctx->sw[0])); 14 if (corr_ctx->symm) { 15 CHKERRQ(MatMult(corr_ctx->inv_smat,corr_ctx->sw[0],corr_ctx->sw[1])); 16 } else { 17 CHKERRQ(MatMultTranspose(corr_ctx->inv_smat,corr_ctx->sw[0],corr_ctx->sw[1])); 18 } 19 CHKERRQ(VecScale(corr_ctx->sw[1],-1.0)); 20 CHKERRQ(MatMult(corr_ctx->basis_mat,corr_ctx->sw[1],corr_ctx->fw[0])); 21 CHKERRQ(VecScale(corr_ctx->sw[1],-1.0)); 22 CHKERRQ(KSPGetOperators(ksp,&K,NULL)); 23 CHKERRQ(MatMultAdd(K,corr_ctx->fw[0],y,y)); 24 CHKERRQ(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 CHKERRQ(PetscLogEventBegin(corr_ctx->evapply,ksp,0,0,0)); 36 CHKERRQ(KSPGetOperators(ksp,&K,NULL)); 37 if (corr_ctx->symm) { 38 CHKERRQ(MatMult(K,x,corr_ctx->fw[0])); 39 } else { 40 CHKERRQ(MatMultTranspose(K,x,corr_ctx->fw[0])); 41 } 42 CHKERRQ(MatMultTranspose(corr_ctx->basis_mat,corr_ctx->fw[0],corr_ctx->sw[0])); 43 CHKERRQ(VecScale(corr_ctx->sw[0],-1.0)); 44 CHKERRQ(MatMult(corr_ctx->inv_smat,corr_ctx->sw[0],corr_ctx->sw[2])); 45 CHKERRQ(MatMultAdd(corr_ctx->basis_mat,corr_ctx->sw[2],x,corr_ctx->fw[0])); 46 CHKERRQ(VecScale(corr_ctx->fw[0],corr_ctx->scale)); 47 /* Sum contributions from approximate solver and projected system */ 48 CHKERRQ(MatMultAdd(corr_ctx->basis_mat,corr_ctx->sw[1],corr_ctx->fw[0],x)); 49 CHKERRQ(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 CHKERRQ(VecDestroyVecs(3,&corr_ctx->sw)); 59 CHKERRQ(VecDestroyVecs(1,&corr_ctx->fw)); 60 CHKERRQ(MatDestroy(&corr_ctx->basis_mat)); 61 CHKERRQ(MatDestroy(&corr_ctx->inv_smat)); 62 CHKERRQ(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 CHKERRQ(KSPGetOperators(local_ksp,&local_mat,&local_pmat)); 83 CHKERRQ(MatGetNearNullSpace(local_pmat,&NullSpace)); 84 if (!NullSpace) { 85 if (pcbddc->dbg_flag) { 86 CHKERRQ(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 CHKERRQ(PetscObjectQuery((PetscObject)NullSpace,"_PBDDC_Null_dmat",(PetscObject*)&dmat)); 91 PetscCheck(dmat,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing dense matrix"); 92 CHKERRQ(PetscLogEventBegin(PC_BDDC_ApproxSetUp[pcbddc->current_level],pc,0,0,0)); 93 94 CHKERRQ(PetscNew(&shell_ctx)); 95 shell_ctx->scale = 1.0; 96 CHKERRQ(PetscObjectReference((PetscObject)dmat)); 97 shell_ctx->basis_mat = dmat; 98 CHKERRQ(MatGetSize(dmat,NULL,&basis_size)); 99 shell_ctx->evapply = PC_BDDC_ApproxApply[pcbddc->current_level]; 100 101 CHKERRQ(MatGetOption(local_mat,MAT_SYMMETRIC,&shell_ctx->symm)); 102 103 /* explicit construct (Phi^T K Phi)^-1 */ 104 CHKERRQ(PetscObjectTypeCompare((PetscObject)local_mat,MATSEQAIJCUSPARSE,&iscusp)); 105 if (iscusp) { 106 CHKERRQ(MatConvert(shell_ctx->basis_mat,MATSEQDENSECUDA,MAT_INPLACE_MATRIX,&shell_ctx->basis_mat)); 107 } 108 CHKERRQ(MatMatMult(local_mat,shell_ctx->basis_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&Kbasis_mat)); 109 CHKERRQ(MatTransposeMatMult(Kbasis_mat,shell_ctx->basis_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&shell_ctx->inv_smat)); 110 CHKERRQ(MatDestroy(&Kbasis_mat)); 111 CHKERRQ(MatBindToCPU(shell_ctx->inv_smat,PETSC_TRUE)); 112 CHKERRQ(MatFindZeroRows(shell_ctx->inv_smat,&zerorows)); 113 if (zerorows) { /* linearly dependent basis */ 114 const PetscInt *idxs; 115 PetscInt i,nz; 116 117 CHKERRQ(ISGetLocalSize(zerorows,&nz)); 118 CHKERRQ(ISGetIndices(zerorows,&idxs)); 119 for (i=0;i<nz;i++) { 120 CHKERRQ(MatSetValue(shell_ctx->inv_smat,idxs[i],idxs[i],1.0,INSERT_VALUES)); 121 } 122 CHKERRQ(ISRestoreIndices(zerorows,&idxs)); 123 CHKERRQ(MatAssemblyBegin(shell_ctx->inv_smat,MAT_FINAL_ASSEMBLY)); 124 CHKERRQ(MatAssemblyEnd(shell_ctx->inv_smat,MAT_FINAL_ASSEMBLY)); 125 } 126 CHKERRQ(MatLUFactor(shell_ctx->inv_smat,NULL,NULL,NULL)); 127 CHKERRQ(MatSeqDenseInvertFactors_Private(shell_ctx->inv_smat)); 128 if (zerorows) { /* linearly dependent basis */ 129 const PetscInt *idxs; 130 PetscInt i,nz; 131 132 CHKERRQ(ISGetLocalSize(zerorows,&nz)); 133 CHKERRQ(ISGetIndices(zerorows,&idxs)); 134 for (i=0;i<nz;i++) { 135 CHKERRQ(MatSetValue(shell_ctx->inv_smat,idxs[i],idxs[i],0.0,INSERT_VALUES)); 136 } 137 CHKERRQ(ISRestoreIndices(zerorows,&idxs)); 138 CHKERRQ(MatAssemblyBegin(shell_ctx->inv_smat,MAT_FINAL_ASSEMBLY)); 139 CHKERRQ(MatAssemblyEnd(shell_ctx->inv_smat,MAT_FINAL_ASSEMBLY)); 140 } 141 CHKERRQ(ISDestroy(&zerorows)); 142 143 /* Create work vectors in shell context */ 144 CHKERRQ(MatCreateVecs(shell_ctx->inv_smat,&v,NULL)); 145 CHKERRQ(KSPCreateVecs(local_ksp,1,&shell_ctx->fw,0,NULL)); 146 CHKERRQ(VecDuplicateVecs(v,3,&shell_ctx->sw)); 147 CHKERRQ(VecDestroy(&v)); 148 149 /* add special pre/post solve to KSP (see [1], eq. 48) */ 150 CHKERRQ(KSPSetPreSolve(local_ksp,PCBDDCNullSpaceCorrPreSolve,shell_ctx)); 151 CHKERRQ(KSPSetPostSolve(local_ksp,PCBDDCNullSpaceCorrPostSolve,shell_ctx)); 152 CHKERRQ(PetscContainerCreate(PetscObjectComm((PetscObject)local_ksp),&c)); 153 CHKERRQ(PetscContainerSetPointer(c,shell_ctx)); 154 CHKERRQ(PetscContainerSetUserDestroy(c,PCBDDCNullSpaceCorrDestroy)); 155 CHKERRQ(PetscObjectCompose((PetscObject)local_ksp,"_PCBDDC_Null_PrePost_ctx",(PetscObject)c)); 156 CHKERRQ(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 CHKERRQ(VecDuplicate(shell_ctx->fw[0],&work1)); 168 CHKERRQ(VecDuplicate(shell_ctx->fw[0],&work2)); 169 CHKERRQ(KSPCreate(PETSC_COMM_SELF,&check_ksp)); 170 if (local_mat->spd) { 171 CHKERRQ(KSPSetType(check_ksp,KSPCG)); 172 } 173 CHKERRQ(PetscObjectIncrementTabLevel((PetscObject)check_ksp,(PetscObject)local_ksp,0)); 174 CHKERRQ(KSPGetOptionsPrefix(local_ksp,&prefix)); 175 CHKERRQ(KSPSetOptionsPrefix(check_ksp,prefix)); 176 CHKERRQ(KSPAppendOptionsPrefix(check_ksp,"approximate_scale_")); 177 CHKERRQ(KSPSetErrorIfNotConverged(check_ksp,PETSC_FALSE)); 178 CHKERRQ(KSPSetOperators(check_ksp,local_mat,local_pmat)); 179 CHKERRQ(KSPSetComputeSingularValues(check_ksp,PETSC_TRUE)); 180 CHKERRQ(KSPSetPreSolve(check_ksp,PCBDDCNullSpaceCorrPreSolve,shell_ctx)); 181 CHKERRQ(KSPSetPostSolve(check_ksp,PCBDDCNullSpaceCorrPostSolve,shell_ctx)); 182 CHKERRQ(KSPSetTolerances(check_ksp,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,PETSC_DEFAULT)); 183 CHKERRQ(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 CHKERRQ(KSPSetUp(check_ksp)); 186 CHKERRQ(KSPGetPC(local_ksp,&local_pc)); 187 CHKERRQ(KSPSetPC(check_ksp,local_pc)); 188 CHKERRQ(KSPGetTolerances(check_ksp,NULL,NULL,NULL,&maxit)); 189 CHKERRQ(KSPSetTolerances(check_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,PetscMin(10,maxit))); 190 CHKERRQ(VecSetRandom(work2,NULL)); 191 CHKERRQ(MatMult(local_mat,work2,work1)); 192 CHKERRQ(KSPSolve(check_ksp,work1,work1)); 193 CHKERRQ(KSPCheckSolve(check_ksp,pc,work1)); 194 CHKERRQ(VecAXPY(work1,-1.,work2)); 195 CHKERRQ(VecNorm(work1,NORM_INFINITY,&test_err)); 196 CHKERRQ(KSPComputeExtremeSingularValues(check_ksp,&lambda_max,&lambda_min)); 197 CHKERRQ(KSPGetIterationNumber(check_ksp,&k)); 198 if (pcbddc->dbg_flag) { 199 if (isdir) { 200 CHKERRQ(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 CHKERRQ(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 CHKERRQ(VecSetRandom(work2,NULL)); 211 CHKERRQ(MatMult(local_mat,work2,work1)); 212 CHKERRQ(PCCreate(PetscObjectComm((PetscObject)check_ksp),&new_pc)); 213 CHKERRQ(PCSetType(new_pc,PCKSP)); 214 CHKERRQ(PCSetOperators(new_pc,local_mat,local_pmat)); 215 CHKERRQ(PCKSPSetKSP(new_pc,local_ksp)); 216 CHKERRQ(KSPSetPC(check_ksp,new_pc)); 217 CHKERRQ(PCDestroy(&new_pc)); 218 CHKERRQ(KSPSetTolerances(check_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,maxit)); 219 CHKERRQ(KSPSetPreSolve(check_ksp,NULL,NULL)); 220 CHKERRQ(KSPSetPostSolve(check_ksp,NULL,NULL)); 221 CHKERRQ(KSPSolve(check_ksp,work1,work1)); 222 CHKERRQ(KSPCheckSolve(check_ksp,pc,work1)); 223 CHKERRQ(VecAXPY(work1,-1.,work2)); 224 CHKERRQ(VecNorm(work1,NORM_INFINITY,&test_err)); 225 CHKERRQ(KSPComputeExtremeSingularValues(check_ksp,&lambda_max,&lambda_min)); 226 CHKERRQ(KSPGetIterationNumber(check_ksp,&k)); 227 if (pcbddc->dbg_flag) { 228 if (isdir) { 229 CHKERRQ(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 CHKERRQ(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 CHKERRQ(KSPDestroy(&check_ksp)); 236 CHKERRQ(VecDestroy(&work1)); 237 CHKERRQ(VecDestroy(&work2)); 238 } 239 CHKERRQ(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 CHKERRQ(VecDuplicate(shell_ctx->fw[0],&work1)); 247 CHKERRQ(VecDuplicate(shell_ctx->fw[0],&work2)); 248 CHKERRQ(VecDuplicate(shell_ctx->fw[0],&work3)); 249 CHKERRQ(VecSetRandom(shell_ctx->sw[0],NULL)); 250 CHKERRQ(MatMult(shell_ctx->basis_mat,shell_ctx->sw[0],work1)); 251 CHKERRQ(VecCopy(work1,work2)); 252 CHKERRQ(MatMult(local_mat,work1,work3)); 253 CHKERRQ(KSPSolve(local_ksp,work3,work1)); 254 CHKERRQ(VecAXPY(work1,-1.,work2)); 255 CHKERRQ(VecNorm(work1,NORM_INFINITY,&test_err)); 256 if (isdir) { 257 CHKERRQ(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Dirichlet nullspace correction solver: %1.14e\n",PetscGlobalRank,test_err)); 258 } else { 259 CHKERRQ(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Neumann nullspace correction solver: %1.14e\n",PetscGlobalRank,test_err)); 260 } 261 CHKERRQ(VecDestroy(&work1)); 262 CHKERRQ(VecDestroy(&work2)); 263 CHKERRQ(VecDestroy(&work3)); 264 } 265 PetscFunctionReturn(0); 266 } 267