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