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, void *ctx) { 7 NullSpaceCorrection_ctx corr_ctx = (NullSpaceCorrection_ctx)ctx; 8 Mat K; 9 10 PetscFunctionBegin; 11 PetscCall(PetscLogEventBegin(corr_ctx->evapply, ksp, 0, 0, 0)); 12 PetscCall(MatMultTranspose(corr_ctx->basis_mat, y, corr_ctx->sw[0])); 13 if (corr_ctx->symm) { 14 PetscCall(MatMult(corr_ctx->inv_smat, corr_ctx->sw[0], corr_ctx->sw[1])); 15 } else { 16 PetscCall(MatMultTranspose(corr_ctx->inv_smat, corr_ctx->sw[0], corr_ctx->sw[1])); 17 } 18 PetscCall(VecScale(corr_ctx->sw[1], -1.0)); 19 PetscCall(MatMult(corr_ctx->basis_mat, corr_ctx->sw[1], corr_ctx->fw[0])); 20 PetscCall(VecScale(corr_ctx->sw[1], -1.0)); 21 PetscCall(KSPGetOperators(ksp, &K, NULL)); 22 PetscCall(MatMultAdd(K, corr_ctx->fw[0], y, y)); 23 PetscCall(PetscLogEventEnd(corr_ctx->evapply, ksp, 0, 0, 0)); 24 PetscFunctionReturn(0); 25 } 26 27 /* E^t + small */ 28 static PetscErrorCode PCBDDCNullSpaceCorrPostSolve(KSP ksp, Vec y, Vec x, void *ctx) { 29 NullSpaceCorrection_ctx corr_ctx = (NullSpaceCorrection_ctx)ctx; 30 Mat K; 31 32 PetscFunctionBegin; 33 PetscCall(PetscLogEventBegin(corr_ctx->evapply, ksp, 0, 0, 0)); 34 PetscCall(KSPGetOperators(ksp, &K, NULL)); 35 if (corr_ctx->symm) { 36 PetscCall(MatMult(K, x, corr_ctx->fw[0])); 37 } else { 38 PetscCall(MatMultTranspose(K, x, corr_ctx->fw[0])); 39 } 40 PetscCall(MatMultTranspose(corr_ctx->basis_mat, corr_ctx->fw[0], corr_ctx->sw[0])); 41 PetscCall(VecScale(corr_ctx->sw[0], -1.0)); 42 PetscCall(MatMult(corr_ctx->inv_smat, corr_ctx->sw[0], corr_ctx->sw[2])); 43 PetscCall(MatMultAdd(corr_ctx->basis_mat, corr_ctx->sw[2], x, corr_ctx->fw[0])); 44 PetscCall(VecScale(corr_ctx->fw[0], corr_ctx->scale)); 45 /* Sum contributions from approximate solver and projected system */ 46 PetscCall(MatMultAdd(corr_ctx->basis_mat, corr_ctx->sw[1], corr_ctx->fw[0], x)); 47 PetscCall(PetscLogEventEnd(corr_ctx->evapply, ksp, 0, 0, 0)); 48 PetscFunctionReturn(0); 49 } 50 51 static PetscErrorCode PCBDDCNullSpaceCorrDestroy(void *ctx) { 52 NullSpaceCorrection_ctx corr_ctx = (NullSpaceCorrection_ctx)ctx; 53 54 PetscFunctionBegin; 55 PetscCall(VecDestroyVecs(3, &corr_ctx->sw)); 56 PetscCall(VecDestroyVecs(1, &corr_ctx->fw)); 57 PetscCall(MatDestroy(&corr_ctx->basis_mat)); 58 PetscCall(MatDestroy(&corr_ctx->inv_smat)); 59 PetscCall(PetscFree(corr_ctx)); 60 PetscFunctionReturn(0); 61 } 62 63 PetscErrorCode PCBDDCNullSpaceAssembleCorrection(PC pc, PetscBool isdir, PetscBool needscaling) { 64 PC_BDDC *pcbddc = (PC_BDDC *)pc->data; 65 MatNullSpace NullSpace = NULL; 66 KSP local_ksp; 67 NullSpaceCorrection_ctx shell_ctx; 68 Mat local_mat, local_pmat, dmat, Kbasis_mat; 69 Vec v; 70 PetscContainer c; 71 PetscInt basis_size; 72 IS zerorows; 73 PetscBool iscusp; 74 75 PetscFunctionBegin; 76 if (isdir) local_ksp = pcbddc->ksp_D; /* Dirichlet solver */ 77 else local_ksp = pcbddc->ksp_R; /* Neumann solver */ 78 PetscCall(KSPGetOperators(local_ksp, &local_mat, &local_pmat)); 79 PetscCall(MatGetNearNullSpace(local_pmat, &NullSpace)); 80 if (!NullSpace) { 81 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")); 82 PetscFunctionReturn(0); 83 } 84 PetscCall(PetscObjectQuery((PetscObject)NullSpace, "_PBDDC_Null_dmat", (PetscObject *)&dmat)); 85 PetscCheck(dmat, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing dense matrix"); 86 PetscCall(PetscLogEventBegin(PC_BDDC_ApproxSetUp[pcbddc->current_level], pc, 0, 0, 0)); 87 88 PetscCall(PetscNew(&shell_ctx)); 89 shell_ctx->scale = 1.0; 90 PetscCall(PetscObjectReference((PetscObject)dmat)); 91 shell_ctx->basis_mat = dmat; 92 PetscCall(MatGetSize(dmat, NULL, &basis_size)); 93 shell_ctx->evapply = PC_BDDC_ApproxApply[pcbddc->current_level]; 94 95 PetscCall(MatIsSymmetric(local_mat, 0.0, &shell_ctx->symm)); 96 97 /* explicit construct (Phi^T K Phi)^-1 */ 98 PetscCall(PetscObjectTypeCompare((PetscObject)local_mat, MATSEQAIJCUSPARSE, &iscusp)); 99 if (iscusp) PetscCall(MatConvert(shell_ctx->basis_mat, MATSEQDENSECUDA, MAT_INPLACE_MATRIX, &shell_ctx->basis_mat)); 100 PetscCall(MatMatMult(local_mat, shell_ctx->basis_mat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Kbasis_mat)); 101 PetscCall(MatTransposeMatMult(Kbasis_mat, shell_ctx->basis_mat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &shell_ctx->inv_smat)); 102 PetscCall(MatDestroy(&Kbasis_mat)); 103 PetscCall(MatBindToCPU(shell_ctx->inv_smat, PETSC_TRUE)); 104 PetscCall(MatFindZeroRows(shell_ctx->inv_smat, &zerorows)); 105 if (zerorows) { /* linearly dependent basis */ 106 const PetscInt *idxs; 107 PetscInt i, nz; 108 109 PetscCall(ISGetLocalSize(zerorows, &nz)); 110 PetscCall(ISGetIndices(zerorows, &idxs)); 111 for (i = 0; i < nz; i++) PetscCall(MatSetValue(shell_ctx->inv_smat, idxs[i], idxs[i], 1.0, INSERT_VALUES)); 112 PetscCall(ISRestoreIndices(zerorows, &idxs)); 113 PetscCall(MatAssemblyBegin(shell_ctx->inv_smat, MAT_FINAL_ASSEMBLY)); 114 PetscCall(MatAssemblyEnd(shell_ctx->inv_smat, MAT_FINAL_ASSEMBLY)); 115 } 116 PetscCall(MatLUFactor(shell_ctx->inv_smat, NULL, NULL, NULL)); 117 PetscCall(MatSeqDenseInvertFactors_Private(shell_ctx->inv_smat)); 118 if (zerorows) { /* linearly dependent basis */ 119 const PetscInt *idxs; 120 PetscInt i, nz; 121 122 PetscCall(ISGetLocalSize(zerorows, &nz)); 123 PetscCall(ISGetIndices(zerorows, &idxs)); 124 for (i = 0; i < nz; i++) PetscCall(MatSetValue(shell_ctx->inv_smat, idxs[i], idxs[i], 0.0, INSERT_VALUES)); 125 PetscCall(ISRestoreIndices(zerorows, &idxs)); 126 PetscCall(MatAssemblyBegin(shell_ctx->inv_smat, MAT_FINAL_ASSEMBLY)); 127 PetscCall(MatAssemblyEnd(shell_ctx->inv_smat, MAT_FINAL_ASSEMBLY)); 128 } 129 PetscCall(ISDestroy(&zerorows)); 130 131 /* Create work vectors in shell context */ 132 PetscCall(MatCreateVecs(shell_ctx->inv_smat, &v, NULL)); 133 PetscCall(KSPCreateVecs(local_ksp, 1, &shell_ctx->fw, 0, NULL)); 134 PetscCall(VecDuplicateVecs(v, 3, &shell_ctx->sw)); 135 PetscCall(VecDestroy(&v)); 136 137 /* add special pre/post solve to KSP (see [1], eq. 48) */ 138 PetscCall(KSPSetPreSolve(local_ksp, PCBDDCNullSpaceCorrPreSolve, shell_ctx)); 139 PetscCall(KSPSetPostSolve(local_ksp, PCBDDCNullSpaceCorrPostSolve, shell_ctx)); 140 PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)local_ksp), &c)); 141 PetscCall(PetscContainerSetPointer(c, shell_ctx)); 142 PetscCall(PetscContainerSetUserDestroy(c, PCBDDCNullSpaceCorrDestroy)); 143 PetscCall(PetscObjectCompose((PetscObject)local_ksp, "_PCBDDC_Null_PrePost_ctx", (PetscObject)c)); 144 PetscCall(PetscContainerDestroy(&c)); 145 146 /* Create ksp object suitable for extreme eigenvalues' estimation */ 147 if (needscaling || pcbddc->dbg_flag) { 148 KSP check_ksp; 149 PC local_pc; 150 Vec work1, work2; 151 const char *prefix; 152 PetscReal test_err, lambda_min, lambda_max; 153 PetscInt k, maxit; 154 PetscBool isspd, isset; 155 156 PetscCall(VecDuplicate(shell_ctx->fw[0], &work1)); 157 PetscCall(VecDuplicate(shell_ctx->fw[0], &work2)); 158 PetscCall(KSPCreate(PETSC_COMM_SELF, &check_ksp)); 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_DEFAULT, PETSC_DEFAULT)); 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_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, 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_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, 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(0); 254 } 255