1 #include <../src/ksp/pc/impls/bddc/bddc.h> 2 #include <../src/ksp/pc/impls/bddc/bddcprivate.h> 3 4 static PetscErrorCode PCBDDCViewNullSpaceCorrectionPC(PC pc,PetscViewer view) 5 { 6 PetscErrorCode ierr; 7 PetscBool isascii; 8 NullSpaceCorrection_ctx pc_ctx; 9 10 PetscFunctionBegin; 11 ierr = PCShellGetContext(pc,(void**)&pc_ctx);CHKERRQ(ierr); 12 ierr = PetscObjectTypeCompare((PetscObject)view,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 13 if (isascii) { 14 ierr = PetscViewerASCIIPrintf(view,"inner preconditioner:\n");CHKERRQ(ierr); 15 ierr = PetscViewerASCIIPushTab(view);CHKERRQ(ierr); 16 ierr = PCView(pc_ctx->local_pc,view);CHKERRQ(ierr); 17 ierr = PetscViewerASCIIPopTab(view);CHKERRQ(ierr); 18 19 ierr = PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 20 21 ierr = PetscViewerASCIIPrintf(view,"Lbasis:\n");CHKERRQ(ierr); 22 ierr = PetscViewerASCIIPushTab(view);CHKERRQ(ierr); 23 ierr = MatView(pc_ctx->Lbasis_mat,view);CHKERRQ(ierr); 24 ierr = PetscViewerASCIIPopTab(view);CHKERRQ(ierr); 25 26 ierr = PetscViewerASCIIPrintf(view,"Kbasis:\n");CHKERRQ(ierr); 27 ierr = PetscViewerASCIIPushTab(view);CHKERRQ(ierr); 28 ierr = MatView(pc_ctx->Kbasis_mat,view);CHKERRQ(ierr); 29 ierr = PetscViewerASCIIPopTab(view);CHKERRQ(ierr); 30 31 ierr = PetscViewerPopFormat(view);CHKERRQ(ierr); 32 } 33 PetscFunctionReturn(0); 34 } 35 36 static PetscErrorCode PCBDDCApplyNullSpaceCorrectionPC(PC pc,Vec x,Vec y) 37 { 38 NullSpaceCorrection_ctx pc_ctx; 39 PetscErrorCode ierr; 40 41 PetscFunctionBegin; 42 ierr = PCShellGetContext(pc,(void**)&pc_ctx);CHKERRQ(ierr); 43 /* E */ 44 ierr = MatMultTranspose(pc_ctx->Lbasis_mat,x,pc_ctx->work_small_2);CHKERRQ(ierr); 45 ierr = MatMultAdd(pc_ctx->Kbasis_mat,pc_ctx->work_small_2,x,pc_ctx->work_full_1);CHKERRQ(ierr); 46 /* P^-1 */ 47 ierr = PCApply(pc_ctx->local_pc,pc_ctx->work_full_1,pc_ctx->work_full_2);CHKERRQ(ierr); 48 /* E^T */ 49 ierr = MatMultTranspose(pc_ctx->Kbasis_mat,pc_ctx->work_full_2,pc_ctx->work_small_1);CHKERRQ(ierr); 50 ierr = VecScale(pc_ctx->work_small_1,-1.0);CHKERRQ(ierr); 51 ierr = MatMultAdd(pc_ctx->Lbasis_mat,pc_ctx->work_small_1,pc_ctx->work_full_2,pc_ctx->work_full_1);CHKERRQ(ierr); 52 /* Sum contributions */ 53 ierr = MatMultAdd(pc_ctx->basis_mat,pc_ctx->work_small_2,pc_ctx->work_full_1,y);CHKERRQ(ierr); 54 if (pc_ctx->apply_scaling) { 55 ierr = VecScale(y,pc_ctx->scale);CHKERRQ(ierr); 56 } 57 PetscFunctionReturn(0); 58 } 59 60 static PetscErrorCode PCBDDCDestroyNullSpaceCorrectionPC(PC pc) 61 { 62 NullSpaceCorrection_ctx pc_ctx; 63 PetscErrorCode ierr; 64 65 PetscFunctionBegin; 66 ierr = PCShellGetContext(pc,(void**)&pc_ctx);CHKERRQ(ierr); 67 ierr = VecDestroy(&pc_ctx->work_small_1);CHKERRQ(ierr); 68 ierr = VecDestroy(&pc_ctx->work_small_2);CHKERRQ(ierr); 69 ierr = VecDestroy(&pc_ctx->work_full_1);CHKERRQ(ierr); 70 ierr = VecDestroy(&pc_ctx->work_full_2);CHKERRQ(ierr); 71 ierr = MatDestroy(&pc_ctx->basis_mat);CHKERRQ(ierr); 72 ierr = MatDestroy(&pc_ctx->Lbasis_mat);CHKERRQ(ierr); 73 ierr = MatDestroy(&pc_ctx->Kbasis_mat);CHKERRQ(ierr); 74 ierr = PCDestroy(&pc_ctx->local_pc);CHKERRQ(ierr); 75 ierr = PetscFree(pc_ctx);CHKERRQ(ierr); 76 PetscFunctionReturn(0); 77 } 78 79 PetscErrorCode PCBDDCNullSpaceAssembleCorrection(PC pc, PetscBool isdir, PetscBool needscaling) 80 { 81 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 82 PC_IS *pcis = (PC_IS*)pc->data; 83 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 84 MatNullSpace NullSpace = NULL; 85 KSP local_ksp; 86 PC newpc; 87 NullSpaceCorrection_ctx shell_ctx; 88 Mat local_mat,local_pmat,small_mat,inv_small_mat; 89 Vec *nullvecs; 90 VecScatter scatter_ctx; 91 IS is_aux,local_dofs; 92 MatFactorInfo matinfo; 93 PetscScalar *basis_mat,*Kbasis_mat,*array,*array_mat; 94 PetscScalar one = 1.0,zero = 0.0, m_one = -1.0; 95 PetscInt basis_dofs,basis_size,nnsp_size,i,k; 96 PetscBool nnsp_has_cnst; 97 PetscReal test_err,lambda_min,lambda_max; 98 PetscErrorCode ierr; 99 100 PetscFunctionBegin; 101 ierr = MatGetNullSpace(matis->A,&NullSpace);CHKERRQ(ierr); 102 if (!NullSpace) { 103 ierr = MatGetNearNullSpace(matis->A,&NullSpace);CHKERRQ(ierr); 104 } 105 if (!NullSpace) { 106 if (pcbddc->dbg_flag) { 107 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d doesn't have local (near) nullspace: no need for correction in %s solver \n",PetscGlobalRank,isdir ? "Dirichlet" : "Neumann");CHKERRQ(ierr); 108 } 109 PetscFunctionReturn(0); 110 } 111 112 /* Infer the local solver */ 113 if (isdir) { 114 /* Dirichlet solver */ 115 local_ksp = pcbddc->ksp_D; 116 local_dofs = pcis->is_I_local; 117 } else { 118 /* Neumann solver */ 119 local_ksp = pcbddc->ksp_R; 120 local_dofs = pcbddc->is_R_local; 121 } 122 ierr = ISGetSize(local_dofs,&basis_dofs);CHKERRQ(ierr); 123 ierr = KSPGetOperators(local_ksp,&local_mat,&local_pmat);CHKERRQ(ierr); 124 125 /* Get null space vecs */ 126 ierr = MatNullSpaceGetVecs(NullSpace,&nnsp_has_cnst,&nnsp_size,(const Vec**)&nullvecs);CHKERRQ(ierr); 127 basis_size = nnsp_size; 128 if (nnsp_has_cnst) basis_size++; 129 130 /* Create shell ctx */ 131 ierr = PetscNew(&shell_ctx);CHKERRQ(ierr); 132 shell_ctx->apply_scaling = needscaling; 133 shell_ctx->scale = 1.0; 134 135 /* Create work vectors in shell context */ 136 ierr = VecCreate(PETSC_COMM_SELF,&shell_ctx->work_small_1);CHKERRQ(ierr); 137 ierr = VecSetSizes(shell_ctx->work_small_1,basis_size,basis_size);CHKERRQ(ierr); 138 ierr = VecSetType(shell_ctx->work_small_1,VECSEQ);CHKERRQ(ierr); 139 ierr = VecDuplicate(shell_ctx->work_small_1,&shell_ctx->work_small_2);CHKERRQ(ierr); 140 ierr = VecCreate(PETSC_COMM_SELF,&shell_ctx->work_full_1);CHKERRQ(ierr); 141 ierr = VecSetSizes(shell_ctx->work_full_1,basis_dofs,basis_dofs);CHKERRQ(ierr); 142 ierr = VecSetType(shell_ctx->work_full_1,VECSEQ);CHKERRQ(ierr); 143 ierr = VecDuplicate(shell_ctx->work_full_1,&shell_ctx->work_full_2);CHKERRQ(ierr); 144 145 /* Allocate workspace */ 146 ierr = MatCreateSeqDense(PETSC_COMM_SELF,basis_dofs,basis_size,NULL,&shell_ctx->basis_mat);CHKERRQ(ierr); 147 ierr = MatCreateSeqDense(PETSC_COMM_SELF,basis_dofs,basis_size,NULL,&shell_ctx->Kbasis_mat);CHKERRQ(ierr); 148 ierr = MatDenseGetArray(shell_ctx->basis_mat,&basis_mat);CHKERRQ(ierr); 149 ierr = MatDenseGetArray(shell_ctx->Kbasis_mat,&Kbasis_mat);CHKERRQ(ierr); 150 151 /* Restrict local null space on selected dofs 152 and compute matrices N and K*N */ 153 ierr = VecScatterCreate(pcis->vec1_N,local_dofs,shell_ctx->work_full_1,(IS)0,&scatter_ctx);CHKERRQ(ierr); 154 for (k=0;k<nnsp_size;k++) { 155 ierr = VecPlaceArray(shell_ctx->work_full_1,(const PetscScalar*)&basis_mat[k*basis_dofs]);CHKERRQ(ierr); 156 ierr = VecScatterBegin(scatter_ctx,nullvecs[k],shell_ctx->work_full_1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 157 ierr = VecScatterEnd(scatter_ctx,nullvecs[k],shell_ctx->work_full_1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 158 ierr = VecResetArray(shell_ctx->work_full_1);CHKERRQ(ierr); 159 } 160 if (nnsp_has_cnst) { 161 ierr = VecPlaceArray(shell_ctx->work_full_1,(const PetscScalar*)&basis_mat[k*basis_dofs]);CHKERRQ(ierr); 162 ierr = VecSet(shell_ctx->work_full_1,one);CHKERRQ(ierr); 163 ierr = VecResetArray(shell_ctx->work_full_1);CHKERRQ(ierr); 164 } 165 166 ierr = PetscMalloc1(basis_size,&nullvecs);CHKERRQ(ierr); 167 for (k=0;k<basis_size;k++) { 168 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,basis_dofs,basis_mat + k*basis_dofs,&nullvecs[k]);CHKERRQ(ierr); 169 } 170 ierr = PCBDDCOrthonormalizeVecs(basis_size,nullvecs);CHKERRQ(ierr); 171 ierr = MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_FALSE,basis_size,nullvecs,&NullSpace);CHKERRQ(ierr); 172 ierr = MatSetNearNullSpace(local_mat,NullSpace);CHKERRQ(ierr); 173 ierr = MatNullSpaceDestroy(&NullSpace);CHKERRQ(ierr); 174 for (k=0;k<basis_size;k++) { 175 ierr = VecDestroy(&nullvecs[k]);CHKERRQ(ierr); 176 } 177 ierr = PetscFree(nullvecs);CHKERRQ(ierr); 178 179 for (k=0;k<basis_size;k++) { 180 ierr = VecPlaceArray(shell_ctx->work_full_1,(const PetscScalar*)&basis_mat[k*basis_dofs]);CHKERRQ(ierr); 181 ierr = VecPlaceArray(shell_ctx->work_full_2,(const PetscScalar*)&Kbasis_mat[k*basis_dofs]);CHKERRQ(ierr); 182 ierr = MatMult(local_mat,shell_ctx->work_full_1,shell_ctx->work_full_2);CHKERRQ(ierr); 183 ierr = VecResetArray(shell_ctx->work_full_1);CHKERRQ(ierr); 184 ierr = VecResetArray(shell_ctx->work_full_2);CHKERRQ(ierr); 185 } 186 187 ierr = VecScatterDestroy(&scatter_ctx);CHKERRQ(ierr); 188 ierr = MatDenseRestoreArray(shell_ctx->basis_mat,&basis_mat);CHKERRQ(ierr); 189 ierr = MatDenseRestoreArray(shell_ctx->Kbasis_mat,&Kbasis_mat);CHKERRQ(ierr); 190 /* Assemble another Mat object in shell context */ 191 ierr = MatTransposeMatMult(shell_ctx->basis_mat,shell_ctx->Kbasis_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&small_mat);CHKERRQ(ierr); 192 ierr = MatFactorInfoInitialize(&matinfo);CHKERRQ(ierr); 193 ierr = ISCreateStride(PETSC_COMM_SELF,basis_size,0,1,&is_aux);CHKERRQ(ierr); 194 ierr = MatLUFactor(small_mat,is_aux,is_aux,&matinfo);CHKERRQ(ierr); 195 ierr = ISDestroy(&is_aux);CHKERRQ(ierr); 196 ierr = PetscMalloc1(basis_size*basis_size,&array_mat);CHKERRQ(ierr); 197 for (k=0;k<basis_size;k++) { 198 ierr = VecSet(shell_ctx->work_small_1,zero);CHKERRQ(ierr); 199 ierr = VecSetValue(shell_ctx->work_small_1,k,one,INSERT_VALUES);CHKERRQ(ierr); 200 ierr = VecAssemblyBegin(shell_ctx->work_small_1);CHKERRQ(ierr); 201 ierr = VecAssemblyEnd(shell_ctx->work_small_1);CHKERRQ(ierr); 202 ierr = MatSolve(small_mat,shell_ctx->work_small_1,shell_ctx->work_small_2);CHKERRQ(ierr); 203 ierr = VecGetArrayRead(shell_ctx->work_small_2,(const PetscScalar**)&array);CHKERRQ(ierr); 204 for (i=0;i<basis_size;i++) { 205 array_mat[i*basis_size+k]=array[i]; 206 } 207 ierr = VecRestoreArrayRead(shell_ctx->work_small_2,(const PetscScalar**)&array);CHKERRQ(ierr); 208 } 209 ierr = MatCreateSeqDense(PETSC_COMM_SELF,basis_size,basis_size,array_mat,&inv_small_mat);CHKERRQ(ierr); 210 ierr = MatMatMult(shell_ctx->basis_mat,inv_small_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&shell_ctx->Lbasis_mat);CHKERRQ(ierr); 211 ierr = PetscFree(array_mat);CHKERRQ(ierr); 212 ierr = MatDestroy(&inv_small_mat);CHKERRQ(ierr); 213 ierr = MatDestroy(&small_mat);CHKERRQ(ierr); 214 ierr = MatScale(shell_ctx->Kbasis_mat,m_one);CHKERRQ(ierr); 215 216 /* Rebuild local PC */ 217 ierr = KSPGetPC(local_ksp,&shell_ctx->local_pc);CHKERRQ(ierr); 218 ierr = PetscObjectReference((PetscObject)shell_ctx->local_pc);CHKERRQ(ierr); 219 ierr = PCCreate(PETSC_COMM_SELF,&newpc);CHKERRQ(ierr); 220 ierr = PCSetOperators(newpc,local_mat,local_mat);CHKERRQ(ierr); 221 ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr); 222 ierr = PCShellSetContext(newpc,shell_ctx);CHKERRQ(ierr); 223 ierr = PCShellSetApply(newpc,PCBDDCApplyNullSpaceCorrectionPC);CHKERRQ(ierr); 224 ierr = PCShellSetView(newpc,PCBDDCViewNullSpaceCorrectionPC);CHKERRQ(ierr); 225 ierr = PCShellSetDestroy(newpc,PCBDDCDestroyNullSpaceCorrectionPC);CHKERRQ(ierr); 226 ierr = PCSetUp(newpc);CHKERRQ(ierr); 227 ierr = KSPSetPC(local_ksp,newpc);CHKERRQ(ierr); 228 ierr = KSPSetUp(local_ksp);CHKERRQ(ierr); 229 230 /* Create ksp object suitable for extreme eigenvalues' estimation */ 231 if (needscaling) { 232 KSP check_ksp; 233 Vec *workv; 234 const char* prefix; 235 236 ierr = KSPCreate(PETSC_COMM_SELF,&check_ksp);CHKERRQ(ierr); 237 ierr = KSPGetOptionsPrefix(local_ksp,&prefix);CHKERRQ(ierr); 238 ierr = KSPSetOptionsPrefix(check_ksp,prefix);CHKERRQ(ierr); 239 ierr = KSPAppendOptionsPrefix(check_ksp,"approxscale_");CHKERRQ(ierr); 240 ierr = KSPSetErrorIfNotConverged(check_ksp,PETSC_FALSE);CHKERRQ(ierr); 241 ierr = KSPSetOperators(check_ksp,local_mat,local_mat);CHKERRQ(ierr); 242 ierr = KSPSetTolerances(check_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,10);CHKERRQ(ierr); 243 ierr = KSPSetComputeSingularValues(check_ksp,PETSC_TRUE);CHKERRQ(ierr); 244 ierr = VecDuplicateVecs(shell_ctx->work_full_1,2,&workv);CHKERRQ(ierr); 245 ierr = KSPSetFromOptions(check_ksp);CHKERRQ(ierr); 246 ierr = KSPSetPC(check_ksp,newpc);CHKERRQ(ierr); 247 ierr = KSPSetUp(check_ksp);CHKERRQ(ierr); 248 ierr = VecSetRandom(workv[1],NULL);CHKERRQ(ierr); 249 ierr = MatMult(local_mat,workv[1],workv[0]);CHKERRQ(ierr); 250 ierr = KSPSolve(check_ksp,workv[0],workv[0]);CHKERRQ(ierr); 251 ierr = VecAXPY(workv[0],-1.,workv[1]);CHKERRQ(ierr); 252 ierr = VecNorm(workv[0],NORM_INFINITY,&test_err);CHKERRQ(ierr); 253 ierr = KSPComputeExtremeSingularValues(check_ksp,&lambda_max,&lambda_min);CHKERRQ(ierr); 254 ierr = KSPGetIterationNumber(check_ksp,&k);CHKERRQ(ierr); 255 if (pcbddc->dbg_flag) { 256 if (isdir) { 257 ierr = 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);CHKERRQ(ierr); 258 } else { 259 ierr = 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);CHKERRQ(ierr); 260 } 261 } 262 shell_ctx->scale = 1./lambda_max; 263 ierr = KSPDestroy(&check_ksp);CHKERRQ(ierr); 264 ierr = VecDestroyVecs(2,&workv);CHKERRQ(ierr); 265 } 266 ierr = PCDestroy(&newpc);CHKERRQ(ierr); 267 PetscFunctionReturn(0); 268 } 269 270 PetscErrorCode PCBDDCNullSpaceCheckCorrection(PC pc, PetscBool isdir) 271 { 272 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 273 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 274 KSP check_ksp,local_ksp; 275 MatNullSpace NullSpace = NULL; 276 NullSpaceCorrection_ctx shell_ctx; 277 PC check_pc; 278 Mat test_mat,local_mat; 279 PetscReal test_err,lambda_min,lambda_max; 280 Vec work1,work2,work3; 281 PetscInt k; 282 const char *prefix; 283 PetscErrorCode ierr; 284 285 PetscFunctionBegin; 286 ierr = MatGetNullSpace(matis->A,&NullSpace);CHKERRQ(ierr); 287 if (!NullSpace) { 288 ierr = MatGetNearNullSpace(matis->A,&NullSpace);CHKERRQ(ierr); 289 } 290 if (!NullSpace) { 291 PetscFunctionReturn(0); 292 } 293 if (!pcbddc->dbg_flag) PetscFunctionReturn(0); 294 if (isdir) local_ksp = pcbddc->ksp_D; 295 else local_ksp = pcbddc->ksp_R; 296 ierr = KSPGetOperators(local_ksp,&local_mat,NULL);CHKERRQ(ierr); 297 ierr = KSPGetPC(local_ksp,&check_pc);CHKERRQ(ierr); 298 ierr = PCShellGetContext(check_pc,(void**)&shell_ctx);CHKERRQ(ierr); 299 ierr = VecDuplicate(shell_ctx->work_full_1,&work1);CHKERRQ(ierr); 300 ierr = VecDuplicate(shell_ctx->work_full_1,&work2);CHKERRQ(ierr); 301 ierr = VecDuplicate(shell_ctx->work_full_1,&work3);CHKERRQ(ierr); 302 303 ierr = VecSetRandom(shell_ctx->work_small_1,NULL);CHKERRQ(ierr); 304 ierr = MatMult(shell_ctx->basis_mat,shell_ctx->work_small_1,work1);CHKERRQ(ierr); 305 ierr = VecCopy(work1,work2);CHKERRQ(ierr); 306 ierr = MatMult(local_mat,work1,work3);CHKERRQ(ierr); 307 ierr = PCApply(check_pc,work3,work1);CHKERRQ(ierr); 308 ierr = VecAXPY(work1,-1.,work2);CHKERRQ(ierr); 309 ierr = VecNorm(work1,NORM_INFINITY,&test_err);CHKERRQ(ierr); 310 if (isdir) { 311 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Dirichlet nullspace correction solver: %1.14e\n",PetscGlobalRank,test_err);CHKERRQ(ierr); 312 } else { 313 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Neumann nullspace correction solver: %1.14e\n",PetscGlobalRank,test_err);CHKERRQ(ierr); 314 } 315 316 ierr = MatTransposeMatMult(shell_ctx->Lbasis_mat,shell_ctx->Kbasis_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&test_mat);CHKERRQ(ierr); 317 ierr = MatShift(test_mat,1.);CHKERRQ(ierr); 318 ierr = MatNorm(test_mat,NORM_INFINITY,&test_err);CHKERRQ(ierr); 319 ierr = MatDestroy(&test_mat);CHKERRQ(ierr); 320 if (isdir) { 321 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Dirichlet nullspace matrices: %1.14e\n",PetscGlobalRank,test_err);CHKERRQ(ierr); 322 } else { 323 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Neumann nullspace matrices: %1.14e\n",PetscGlobalRank,test_err);CHKERRQ(ierr); 324 } 325 326 /* Create ksp object suitable for extreme eigenvalues' estimation */ 327 ierr = KSPCreate(PETSC_COMM_SELF,&check_ksp);CHKERRQ(ierr); 328 ierr = KSPGetOptionsPrefix(local_ksp,&prefix);CHKERRQ(ierr); 329 ierr = KSPSetOptionsPrefix(check_ksp,prefix);CHKERRQ(ierr); 330 ierr = KSPAppendOptionsPrefix(check_ksp,"approxcheck_");CHKERRQ(ierr); 331 ierr = KSPSetErrorIfNotConverged(check_ksp,pc->erroriffailure);CHKERRQ(ierr); 332 ierr = KSPSetOperators(check_ksp,local_mat,local_mat);CHKERRQ(ierr); 333 ierr = KSPSetTolerances(check_ksp,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 334 ierr = KSPSetComputeSingularValues(check_ksp,PETSC_TRUE);CHKERRQ(ierr); 335 ierr = KSPSetFromOptions(check_ksp);CHKERRQ(ierr); 336 ierr = KSPSetPC(check_ksp,check_pc);CHKERRQ(ierr); 337 ierr = KSPSetUp(check_ksp);CHKERRQ(ierr); 338 ierr = VecSetRandom(work1,NULL);CHKERRQ(ierr); 339 ierr = MatMult(local_mat,work1,work2);CHKERRQ(ierr); 340 ierr = KSPSolve(check_ksp,work2,work2);CHKERRQ(ierr); 341 ierr = VecAXPY(work2,-1.,work1);CHKERRQ(ierr); 342 ierr = VecNorm(work2,NORM_INFINITY,&test_err);CHKERRQ(ierr); 343 ierr = KSPComputeExtremeSingularValues(check_ksp,&lambda_max,&lambda_min);CHKERRQ(ierr); 344 ierr = KSPGetIterationNumber(check_ksp,&k);CHKERRQ(ierr); 345 if (isdir) { 346 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Dirichlet adapted KSP (scale %d) %1.14e (it %d, eigs %1.6e %1.6e)\n",PetscGlobalRank,shell_ctx->apply_scaling,test_err,k,lambda_min,lambda_max);CHKERRQ(ierr); 347 } else { 348 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Neumann adapted KSP (scale %d) %1.14e (it %d, eigs %1.6e %1.6e)\n",PetscGlobalRank,shell_ctx->apply_scaling,test_err,k,lambda_min,lambda_max);CHKERRQ(ierr); 349 } 350 ierr = KSPDestroy(&check_ksp);CHKERRQ(ierr); 351 ierr = VecDestroy(&work1);CHKERRQ(ierr); 352 ierr = VecDestroy(&work2);CHKERRQ(ierr); 353 ierr = VecDestroy(&work3);CHKERRQ(ierr); 354 PetscFunctionReturn(0); 355 } 356