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