1 #include <petsc/private/pcbddcimpl.h> 2 #include <petsc/private/pcbddcprivateimpl.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, PetscCtx 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(PETSC_SUCCESS); 26 } 27 28 /* E^t + small */ 29 static PetscErrorCode PCBDDCNullSpaceCorrPostSolve(KSP ksp, Vec y, Vec x, PetscCtx 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(PETSC_SUCCESS); 51 } 52 53 static PetscErrorCode PCBDDCNullSpaceCorrDestroy(PetscCtxRt 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(PETSC_SUCCESS); 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 PetscInt basis_size; 75 IS zerorows; 76 PetscBool iscusp; 77 78 PetscFunctionBegin; 79 if (isdir) local_ksp = pcbddc->ksp_D; /* Dirichlet solver */ 80 else local_ksp = pcbddc->ksp_R; /* Neumann solver */ 81 PetscCall(KSPGetOperators(local_ksp, &local_mat, &local_pmat)); 82 PetscCall(MatGetNearNullSpace(local_pmat, &NullSpace)); 83 if (!NullSpace) { 84 if (pcbddc->dbg_flag) 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")); 85 PetscFunctionReturn(PETSC_SUCCESS); 86 } 87 PetscCall(PetscObjectQuery((PetscObject)NullSpace, "_PBDDC_Null_dmat", (PetscObject *)&dmat)); 88 PetscCheck(dmat, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing dense matrix"); 89 PetscCall(PetscLogEventBegin(PC_BDDC_ApproxSetUp[pcbddc->current_level], pc, 0, 0, 0)); 90 91 PetscCall(PetscNew(&shell_ctx)); 92 shell_ctx->scale = 1.0; 93 PetscCall(PetscObjectReference((PetscObject)dmat)); 94 shell_ctx->basis_mat = dmat; 95 PetscCall(MatGetSize(dmat, NULL, &basis_size)); 96 shell_ctx->evapply = PC_BDDC_ApproxApply[pcbddc->current_level]; 97 98 PetscCall(MatIsSymmetric(local_mat, 0.0, &shell_ctx->symm)); 99 100 /* explicit construct (Phi^T K Phi)^-1 */ 101 PetscCall(PetscObjectTypeCompare((PetscObject)local_mat, MATSEQAIJCUSPARSE, &iscusp)); 102 if (iscusp) PetscCall(MatConvert(shell_ctx->basis_mat, MATSEQDENSECUDA, MAT_INPLACE_MATRIX, &shell_ctx->basis_mat)); 103 PetscCall(MatMatMult(local_mat, shell_ctx->basis_mat, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Kbasis_mat)); 104 PetscCall(MatTransposeMatMult(Kbasis_mat, shell_ctx->basis_mat, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &shell_ctx->inv_smat)); 105 PetscCall(MatDestroy(&Kbasis_mat)); 106 PetscCall(MatBindToCPU(shell_ctx->inv_smat, PETSC_TRUE)); 107 PetscCall(MatFindZeroRows(shell_ctx->inv_smat, &zerorows)); 108 if (zerorows) { /* linearly dependent basis */ 109 const PetscInt *idxs; 110 PetscInt i, nz; 111 112 PetscCall(ISGetLocalSize(zerorows, &nz)); 113 PetscCall(ISGetIndices(zerorows, &idxs)); 114 for (i = 0; i < nz; i++) PetscCall(MatSetValue(shell_ctx->inv_smat, idxs[i], idxs[i], 1.0, INSERT_VALUES)); 115 PetscCall(ISRestoreIndices(zerorows, &idxs)); 116 PetscCall(MatAssemblyBegin(shell_ctx->inv_smat, MAT_FINAL_ASSEMBLY)); 117 PetscCall(MatAssemblyEnd(shell_ctx->inv_smat, MAT_FINAL_ASSEMBLY)); 118 } 119 PetscCall(MatLUFactor(shell_ctx->inv_smat, NULL, NULL, NULL)); 120 PetscCall(MatSeqDenseInvertFactors_Private(shell_ctx->inv_smat)); 121 if (zerorows) { /* linearly dependent basis */ 122 const PetscInt *idxs; 123 PetscInt i, nz; 124 125 PetscCall(ISGetLocalSize(zerorows, &nz)); 126 PetscCall(ISGetIndices(zerorows, &idxs)); 127 for (i = 0; i < nz; i++) PetscCall(MatSetValue(shell_ctx->inv_smat, idxs[i], idxs[i], 0.0, INSERT_VALUES)); 128 PetscCall(ISRestoreIndices(zerorows, &idxs)); 129 PetscCall(MatAssemblyBegin(shell_ctx->inv_smat, MAT_FINAL_ASSEMBLY)); 130 PetscCall(MatAssemblyEnd(shell_ctx->inv_smat, MAT_FINAL_ASSEMBLY)); 131 } 132 PetscCall(ISDestroy(&zerorows)); 133 134 /* Create work vectors in shell context */ 135 PetscCall(MatCreateVecs(shell_ctx->inv_smat, &v, NULL)); 136 PetscCall(KSPCreateVecs(local_ksp, 1, &shell_ctx->fw, 0, NULL)); 137 PetscCall(VecDuplicateVecs(v, 3, &shell_ctx->sw)); 138 PetscCall(VecDestroy(&v)); 139 140 /* add special pre/post solve to KSP (see [1], eq. 48) */ 141 PetscCall(KSPSetPreSolve(local_ksp, PCBDDCNullSpaceCorrPreSolve, shell_ctx)); 142 PetscCall(KSPSetPostSolve(local_ksp, PCBDDCNullSpaceCorrPostSolve, shell_ctx)); 143 PetscCall(PetscObjectContainerCompose((PetscObject)local_ksp, "_PCBDDC_Null_PrePost_ctx", shell_ctx, PCBDDCNullSpaceCorrDestroy)); 144 145 /* Create ksp object suitable for extreme eigenvalues' estimation */ 146 if (needscaling || pcbddc->dbg_flag) { 147 KSP check_ksp; 148 PC local_pc; 149 Vec work1, work2; 150 const char *prefix; 151 PetscReal test_err, lambda_min, lambda_max; 152 PetscInt k, maxit; 153 PetscBool isspd, isset; 154 155 PetscCall(VecDuplicate(shell_ctx->fw[0], &work1)); 156 PetscCall(VecDuplicate(shell_ctx->fw[0], &work2)); 157 PetscCall(KSPCreate(PETSC_COMM_SELF, &check_ksp)); 158 PetscCall(KSPSetNestLevel(check_ksp, pc->kspnestlevel)); 159 PetscCall(MatIsSPDKnown(local_mat, &isset, &isspd)); 160 if (isset && isspd) PetscCall(KSPSetType(check_ksp, KSPCG)); 161 PetscCall(PetscObjectIncrementTabLevel((PetscObject)check_ksp, (PetscObject)local_ksp, 0)); 162 PetscCall(KSPGetOptionsPrefix(local_ksp, &prefix)); 163 PetscCall(KSPSetOptionsPrefix(check_ksp, prefix)); 164 PetscCall(KSPAppendOptionsPrefix(check_ksp, "approximate_scale_")); 165 PetscCall(KSPSetErrorIfNotConverged(check_ksp, PETSC_FALSE)); 166 PetscCall(KSPSetOperators(check_ksp, local_mat, local_pmat)); 167 PetscCall(KSPSetComputeSingularValues(check_ksp, PETSC_TRUE)); 168 PetscCall(KSPSetPreSolve(check_ksp, PCBDDCNullSpaceCorrPreSolve, shell_ctx)); 169 PetscCall(KSPSetPostSolve(check_ksp, PCBDDCNullSpaceCorrPostSolve, shell_ctx)); 170 PetscCall(KSPSetTolerances(check_ksp, PETSC_SMALL, PETSC_SMALL, PETSC_CURRENT, PETSC_CURRENT)); 171 PetscCall(KSPSetFromOptions(check_ksp)); 172 /* setup with default maxit, then set maxit to min(10,any_set_from_command_line) (bug in computing eigenvalues when changing the number of iterations */ 173 PetscCall(KSPSetUp(check_ksp)); 174 PetscCall(KSPGetPC(local_ksp, &local_pc)); 175 PetscCall(KSPSetPC(check_ksp, local_pc)); 176 PetscCall(KSPGetTolerances(check_ksp, NULL, NULL, NULL, &maxit)); 177 PetscCall(KSPSetTolerances(check_ksp, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, PetscMin(10, maxit))); 178 PetscCall(VecSetRandom(work2, NULL)); 179 PetscCall(MatMult(local_mat, work2, work1)); 180 PetscCall(KSPSolve(check_ksp, work1, work1)); 181 PetscCall(KSPCheckSolve(check_ksp, pc, work1)); 182 PetscCall(VecAXPY(work1, -1., work2)); 183 PetscCall(VecNorm(work1, NORM_INFINITY, &test_err)); 184 PetscCall(KSPComputeExtremeSingularValues(check_ksp, &lambda_max, &lambda_min)); 185 PetscCall(KSPGetIterationNumber(check_ksp, &k)); 186 if (pcbddc->dbg_flag) { 187 if (isdir) { 188 PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d infinity error for Dirichlet adapted solver (no scale) %1.14e (it %" PetscInt_FMT ", eigs %1.6e %1.6e)\n", PetscGlobalRank, (double)test_err, k, (double)lambda_min, (double)lambda_max)); 189 } else { 190 PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d infinity error for Neumann adapted solver (no scale) %1.14e (it %" PetscInt_FMT ", eigs %1.6e %1.6e)\n", PetscGlobalRank, (double)test_err, k, (double)lambda_min, (double)lambda_max)); 191 } 192 } 193 if (needscaling) shell_ctx->scale = 1.0 / lambda_max; 194 195 if (needscaling && pcbddc->dbg_flag) { /* test for scaling factor */ 196 PC new_pc; 197 198 PetscCall(VecSetRandom(work2, NULL)); 199 PetscCall(MatMult(local_mat, work2, work1)); 200 PetscCall(PCCreate(PetscObjectComm((PetscObject)check_ksp), &new_pc)); 201 PetscCall(PCSetType(new_pc, PCKSP)); 202 PetscCall(PCSetOperators(new_pc, local_mat, local_pmat)); 203 PetscCall(PCKSPSetKSP(new_pc, local_ksp)); 204 PetscCall(KSPSetPC(check_ksp, new_pc)); 205 PetscCall(PCDestroy(&new_pc)); 206 PetscCall(KSPSetTolerances(check_ksp, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, maxit)); 207 PetscCall(KSPSetPreSolve(check_ksp, NULL, NULL)); 208 PetscCall(KSPSetPostSolve(check_ksp, NULL, NULL)); 209 PetscCall(KSPSolve(check_ksp, work1, work1)); 210 PetscCall(KSPCheckSolve(check_ksp, pc, work1)); 211 PetscCall(VecAXPY(work1, -1., work2)); 212 PetscCall(VecNorm(work1, NORM_INFINITY, &test_err)); 213 PetscCall(KSPComputeExtremeSingularValues(check_ksp, &lambda_max, &lambda_min)); 214 PetscCall(KSPGetIterationNumber(check_ksp, &k)); 215 if (pcbddc->dbg_flag) { 216 if (isdir) { 217 PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d infinity error for Dirichlet adapted solver (scale %g) %1.14e (it %" PetscInt_FMT ", eigs %1.6e %1.6e)\n", PetscGlobalRank, (double)PetscRealPart(shell_ctx->scale), (double)test_err, k, (double)lambda_min, (double)lambda_max)); 218 } else { 219 PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d infinity error for Neumann adapted solver (scale %g) %1.14e (it %" PetscInt_FMT ", eigs %1.6e %1.6e)\n", PetscGlobalRank, (double)PetscRealPart(shell_ctx->scale), (double)test_err, k, (double)lambda_min, (double)lambda_max)); 220 } 221 } 222 } 223 PetscCall(KSPDestroy(&check_ksp)); 224 PetscCall(VecDestroy(&work1)); 225 PetscCall(VecDestroy(&work2)); 226 } 227 PetscCall(PetscLogEventEnd(PC_BDDC_ApproxSetUp[pcbddc->current_level], pc, 0, 0, 0)); 228 229 if (pcbddc->dbg_flag) { 230 Vec work1, work2, work3; 231 PetscReal test_err; 232 233 /* check nullspace basis is solved exactly */ 234 PetscCall(VecDuplicate(shell_ctx->fw[0], &work1)); 235 PetscCall(VecDuplicate(shell_ctx->fw[0], &work2)); 236 PetscCall(VecDuplicate(shell_ctx->fw[0], &work3)); 237 PetscCall(VecSetRandom(shell_ctx->sw[0], NULL)); 238 PetscCall(MatMult(shell_ctx->basis_mat, shell_ctx->sw[0], work1)); 239 PetscCall(VecCopy(work1, work2)); 240 PetscCall(MatMult(local_mat, work1, work3)); 241 PetscCall(KSPSolve(local_ksp, work3, work1)); 242 PetscCall(VecAXPY(work1, -1., work2)); 243 PetscCall(VecNorm(work1, NORM_INFINITY, &test_err)); 244 if (isdir) { 245 PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d infinity error for Dirichlet nullspace correction solver: %1.14e\n", PetscGlobalRank, (double)test_err)); 246 } else { 247 PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d infinity error for Neumann nullspace correction solver: %1.14e\n", PetscGlobalRank, (double)test_err)); 248 } 249 PetscCall(VecDestroy(&work1)); 250 PetscCall(VecDestroy(&work2)); 251 PetscCall(VecDestroy(&work3)); 252 } 253 PetscFunctionReturn(PETSC_SUCCESS); 254 } 255