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__ "PCBDDCNullSpaceAssembleCoarse" 6 PetscErrorCode PCBDDCNullSpaceAssembleCoarse(PC pc, Mat coarse_mat, MatNullSpace* CoarseNullSpace) 7 { 8 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 9 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 10 MatNullSpace tempCoarseNullSpace=NULL; 11 const Vec *nsp_vecs; 12 Vec *coarse_nsp_vecs,local_vec,local_primal_vec,wcoarse_vec,wcoarse_rhs; 13 PetscInt nsp_size,coarse_nsp_size,i; 14 PetscBool nsp_has_cnst; 15 PetscReal test_null; 16 PetscErrorCode ierr; 17 18 PetscFunctionBegin; 19 tempCoarseNullSpace = 0; 20 coarse_nsp_size = 0; 21 coarse_nsp_vecs = 0; 22 ierr = MatNullSpaceGetVecs(pcbddc->NullSpace,&nsp_has_cnst,&nsp_size,&nsp_vecs);CHKERRQ(ierr); 23 if (coarse_mat) { 24 ierr = PetscMalloc1(nsp_size+1,&coarse_nsp_vecs);CHKERRQ(ierr); 25 for (i=0;i<nsp_size+1;i++) { 26 ierr = MatCreateVecs(coarse_mat,&coarse_nsp_vecs[i],NULL);CHKERRQ(ierr); 27 } 28 if (pcbddc->dbg_flag) { 29 ierr = MatCreateVecs(coarse_mat,&wcoarse_vec,&wcoarse_rhs);CHKERRQ(ierr); 30 } 31 } 32 ierr = MatCreateVecs(pcbddc->ConstraintMatrix,&local_vec,&local_primal_vec);CHKERRQ(ierr); 33 if (nsp_has_cnst) { 34 ierr = VecSet(local_vec,1.0);CHKERRQ(ierr); 35 ierr = MatMult(pcbddc->ConstraintMatrix,local_vec,local_primal_vec);CHKERRQ(ierr); 36 ierr = VecScatterBegin(pcbddc->coarse_loc_to_glob,local_primal_vec,pcbddc->coarse_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 37 ierr = VecScatterEnd(pcbddc->coarse_loc_to_glob,local_primal_vec,pcbddc->coarse_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 38 if (coarse_mat) { 39 PetscScalar *array_out; 40 const PetscScalar *array_in; 41 PetscInt lsize; 42 if (pcbddc->dbg_flag) { 43 PetscViewer dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)coarse_mat)); 44 ierr = MatMult(coarse_mat,wcoarse_vec,wcoarse_rhs);CHKERRQ(ierr); 45 ierr = VecNorm(wcoarse_rhs,NORM_INFINITY,&test_null);CHKERRQ(ierr); 46 ierr = PetscViewerASCIIPrintf(dbg_viewer,"Constant coarse null space error % 1.14e\n",test_null);CHKERRQ(ierr); 47 ierr = PetscViewerFlush(dbg_viewer);CHKERRQ(ierr); 48 } 49 ierr = VecGetLocalSize(pcbddc->coarse_vec,&lsize);CHKERRQ(ierr); 50 ierr = VecGetArrayRead(pcbddc->coarse_vec,&array_in);CHKERRQ(ierr); 51 ierr = VecGetArray(coarse_nsp_vecs[coarse_nsp_size],&array_out);CHKERRQ(ierr); 52 ierr = PetscMemcpy(array_out,array_in,lsize*sizeof(PetscScalar));CHKERRQ(ierr); 53 ierr = VecRestoreArray(coarse_nsp_vecs[coarse_nsp_size],&array_out);CHKERRQ(ierr); 54 ierr = VecRestoreArrayRead(pcbddc->coarse_vec,&array_in);CHKERRQ(ierr); 55 coarse_nsp_size++; 56 } 57 } 58 for (i=0;i<nsp_size;i++) { 59 ierr = VecScatterBegin(matis->rctx,nsp_vecs[i],local_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 60 ierr = VecScatterEnd(matis->rctx,nsp_vecs[i],local_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 61 ierr = MatMult(pcbddc->ConstraintMatrix,local_vec,local_primal_vec);CHKERRQ(ierr); 62 ierr = VecScatterBegin(pcbddc->coarse_loc_to_glob,local_primal_vec,pcbddc->coarse_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 63 ierr = VecScatterEnd(pcbddc->coarse_loc_to_glob,local_primal_vec,pcbddc->coarse_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 64 if (coarse_mat) { 65 PetscScalar *array_out; 66 const PetscScalar *array_in; 67 PetscInt lsize; 68 if (pcbddc->dbg_flag) { 69 PetscViewer dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)coarse_mat)); 70 ierr = MatMult(coarse_mat,wcoarse_vec,wcoarse_rhs);CHKERRQ(ierr); 71 ierr = VecNorm(wcoarse_rhs,NORM_2,&test_null);CHKERRQ(ierr); 72 ierr = PetscViewerASCIIPrintf(dbg_viewer,"Vec %d coarse null space error % 1.14e\n",i,test_null);CHKERRQ(ierr); 73 ierr = PetscViewerFlush(dbg_viewer);CHKERRQ(ierr); 74 } 75 ierr = VecGetLocalSize(pcbddc->coarse_vec,&lsize);CHKERRQ(ierr); 76 ierr = VecGetArrayRead(pcbddc->coarse_vec,&array_in);CHKERRQ(ierr); 77 ierr = VecGetArray(coarse_nsp_vecs[coarse_nsp_size],&array_out);CHKERRQ(ierr); 78 ierr = PetscMemcpy(array_out,array_in,lsize*sizeof(PetscScalar));CHKERRQ(ierr); 79 ierr = VecRestoreArray(coarse_nsp_vecs[coarse_nsp_size],&array_out);CHKERRQ(ierr); 80 ierr = VecRestoreArrayRead(pcbddc->coarse_vec,&array_in);CHKERRQ(ierr); 81 coarse_nsp_size++; 82 } 83 } 84 if (coarse_nsp_size > 0) { 85 ierr = PCBDDCOrthonormalizeVecs(coarse_nsp_size,coarse_nsp_vecs);CHKERRQ(ierr); 86 ierr = MatNullSpaceCreate(PetscObjectComm((PetscObject)coarse_mat),PETSC_FALSE,coarse_nsp_size,coarse_nsp_vecs,&tempCoarseNullSpace);CHKERRQ(ierr); 87 for (i=0;i<nsp_size+1;i++) { 88 ierr = VecDestroy(&coarse_nsp_vecs[i]);CHKERRQ(ierr); 89 } 90 } 91 if (coarse_mat) { 92 ierr = PetscFree(coarse_nsp_vecs);CHKERRQ(ierr); 93 if (pcbddc->dbg_flag) { 94 ierr = VecDestroy(&wcoarse_vec);CHKERRQ(ierr); 95 ierr = VecDestroy(&wcoarse_rhs);CHKERRQ(ierr); 96 } 97 } 98 ierr = VecDestroy(&local_vec);CHKERRQ(ierr); 99 ierr = VecDestroy(&local_primal_vec);CHKERRQ(ierr); 100 *CoarseNullSpace = tempCoarseNullSpace; 101 PetscFunctionReturn(0); 102 } 103 104 105 #undef __FUNCT__ 106 #define __FUNCT__ "PCBDDCApplyNullSpaceCorrectionPC" 107 static PetscErrorCode PCBDDCApplyNullSpaceCorrectionPC(PC pc,Vec x,Vec y) 108 { 109 NullSpaceCorrection_ctx pc_ctx; 110 PetscErrorCode ierr; 111 112 PetscFunctionBegin; 113 ierr = PCShellGetContext(pc,(void**)&pc_ctx);CHKERRQ(ierr); 114 /* E */ 115 ierr = MatMultTranspose(pc_ctx->Lbasis_mat,x,pc_ctx->work_small_2);CHKERRQ(ierr); 116 ierr = MatMultAdd(pc_ctx->Kbasis_mat,pc_ctx->work_small_2,x,pc_ctx->work_full_1);CHKERRQ(ierr); 117 /* P^-1 */ 118 ierr = PCApply(pc_ctx->local_pc,pc_ctx->work_full_1,pc_ctx->work_full_2);CHKERRQ(ierr); 119 /* E^T */ 120 ierr = MatMultTranspose(pc_ctx->Kbasis_mat,pc_ctx->work_full_2,pc_ctx->work_small_1);CHKERRQ(ierr); 121 ierr = VecScale(pc_ctx->work_small_1,-1.0);CHKERRQ(ierr); 122 ierr = MatMultAdd(pc_ctx->Lbasis_mat,pc_ctx->work_small_1,pc_ctx->work_full_2,pc_ctx->work_full_1);CHKERRQ(ierr); 123 /* Sum contributions */ 124 ierr = MatMultAdd(pc_ctx->basis_mat,pc_ctx->work_small_2,pc_ctx->work_full_1,y);CHKERRQ(ierr); 125 PetscFunctionReturn(0); 126 } 127 128 #undef __FUNCT__ 129 #define __FUNCT__ "PCBDDCDestroyNullSpaceCorrectionPC" 130 static PetscErrorCode PCBDDCDestroyNullSpaceCorrectionPC(PC pc) 131 { 132 NullSpaceCorrection_ctx pc_ctx; 133 PetscErrorCode ierr; 134 135 PetscFunctionBegin; 136 ierr = PCShellGetContext(pc,(void**)&pc_ctx);CHKERRQ(ierr); 137 ierr = VecDestroy(&pc_ctx->work_small_1);CHKERRQ(ierr); 138 ierr = VecDestroy(&pc_ctx->work_small_2);CHKERRQ(ierr); 139 ierr = VecDestroy(&pc_ctx->work_full_1);CHKERRQ(ierr); 140 ierr = VecDestroy(&pc_ctx->work_full_2);CHKERRQ(ierr); 141 ierr = MatDestroy(&pc_ctx->basis_mat);CHKERRQ(ierr); 142 ierr = MatDestroy(&pc_ctx->Lbasis_mat);CHKERRQ(ierr); 143 ierr = MatDestroy(&pc_ctx->Kbasis_mat);CHKERRQ(ierr); 144 ierr = PCDestroy(&pc_ctx->local_pc);CHKERRQ(ierr); 145 ierr = PetscFree(pc_ctx);CHKERRQ(ierr); 146 PetscFunctionReturn(0); 147 } 148 149 /* 150 PETSC_EXTERN PetscErrorCode PCBDDCApplyNullSpaceCorrectionPC(PC,Vec,Vec); 151 PETSC_EXTERN PetscErrorCode PCBDDCDestroyNullSpaceCorrectionPC(PC); 152 */ 153 154 #undef __FUNCT__ 155 #define __FUNCT__ "PCBDDCNullSpaceAssembleCorrection" 156 PetscErrorCode PCBDDCNullSpaceAssembleCorrection(PC pc, PetscBool isdir, IS local_dofs) 157 { 158 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 159 PC_IS *pcis = (PC_IS*)pc->data; 160 Mat_IS* matis = (Mat_IS*)pc->pmat->data; 161 KSP local_ksp; 162 PC newpc; 163 NullSpaceCorrection_ctx shell_ctx; 164 Mat local_mat,local_pmat,small_mat,inv_small_mat; 165 Vec work1,work2; 166 const Vec *nullvecs; 167 VecScatter scatter_ctx; 168 IS is_aux; 169 MatFactorInfo matinfo; 170 PetscScalar *basis_mat,*Kbasis_mat,*array,*array_mat; 171 PetscScalar one = 1.0,zero = 0.0, m_one = -1.0; 172 PetscInt basis_dofs,basis_size,nnsp_size,i,k; 173 PetscBool nnsp_has_cnst; 174 PetscErrorCode ierr; 175 176 PetscFunctionBegin; 177 /* Infer the local solver */ 178 ierr = ISGetSize(local_dofs,&basis_dofs);CHKERRQ(ierr); 179 if (isdir) { 180 /* Dirichlet solver */ 181 local_ksp = pcbddc->ksp_D; 182 } else { 183 /* Neumann solver */ 184 local_ksp = pcbddc->ksp_R; 185 } 186 ierr = KSPGetOperators(local_ksp,&local_mat,&local_pmat);CHKERRQ(ierr); 187 188 /* Get null space vecs */ 189 ierr = MatNullSpaceGetVecs(pcbddc->NullSpace,&nnsp_has_cnst,&nnsp_size,&nullvecs);CHKERRQ(ierr); 190 basis_size = nnsp_size; 191 if (nnsp_has_cnst) { 192 basis_size++; 193 } 194 195 if (basis_dofs) { 196 /* Create shell ctx */ 197 ierr = PetscNew(&shell_ctx);CHKERRQ(ierr); 198 199 /* Create work vectors in shell context */ 200 ierr = VecCreate(PETSC_COMM_SELF,&shell_ctx->work_small_1);CHKERRQ(ierr); 201 ierr = VecSetSizes(shell_ctx->work_small_1,basis_size,basis_size);CHKERRQ(ierr); 202 ierr = VecSetType(shell_ctx->work_small_1,VECSEQ);CHKERRQ(ierr); 203 ierr = VecDuplicate(shell_ctx->work_small_1,&shell_ctx->work_small_2);CHKERRQ(ierr); 204 ierr = VecCreate(PETSC_COMM_SELF,&shell_ctx->work_full_1);CHKERRQ(ierr); 205 ierr = VecSetSizes(shell_ctx->work_full_1,basis_dofs,basis_dofs);CHKERRQ(ierr); 206 ierr = VecSetType(shell_ctx->work_full_1,VECSEQ);CHKERRQ(ierr); 207 ierr = VecDuplicate(shell_ctx->work_full_1,&shell_ctx->work_full_2);CHKERRQ(ierr); 208 209 /* Allocate workspace */ 210 ierr = MatCreateSeqDense(PETSC_COMM_SELF,basis_dofs,basis_size,NULL,&shell_ctx->basis_mat );CHKERRQ(ierr); 211 ierr = MatCreateSeqDense(PETSC_COMM_SELF,basis_dofs,basis_size,NULL,&shell_ctx->Kbasis_mat);CHKERRQ(ierr); 212 ierr = MatDenseGetArray(shell_ctx->basis_mat,&basis_mat);CHKERRQ(ierr); 213 ierr = MatDenseGetArray(shell_ctx->Kbasis_mat,&Kbasis_mat);CHKERRQ(ierr); 214 215 /* Restrict local null space on selected dofs (Dirichlet or Neumann) 216 and compute matrices N and K*N */ 217 ierr = VecDuplicate(shell_ctx->work_full_1,&work1);CHKERRQ(ierr); 218 ierr = VecDuplicate(shell_ctx->work_full_1,&work2);CHKERRQ(ierr); 219 ierr = VecScatterCreate(pcis->vec1_N,local_dofs,work1,(IS)0,&scatter_ctx);CHKERRQ(ierr); 220 } 221 222 for (k=0;k<nnsp_size;k++) { 223 ierr = VecScatterBegin(matis->rctx,nullvecs[k],pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 224 ierr = VecScatterEnd(matis->rctx,nullvecs[k],pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 225 if (basis_dofs) { 226 ierr = VecPlaceArray(work1,(const PetscScalar*)&basis_mat[k*basis_dofs]);CHKERRQ(ierr); 227 ierr = VecScatterBegin(scatter_ctx,pcis->vec1_N,work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 228 ierr = VecScatterEnd(scatter_ctx,pcis->vec1_N,work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 229 ierr = VecPlaceArray(work2,(const PetscScalar*)&Kbasis_mat[k*basis_dofs]);CHKERRQ(ierr); 230 ierr = MatMult(local_mat,work1,work2);CHKERRQ(ierr); 231 ierr = VecResetArray(work1);CHKERRQ(ierr); 232 ierr = VecResetArray(work2);CHKERRQ(ierr); 233 } 234 } 235 236 if (basis_dofs) { 237 if (nnsp_has_cnst) { 238 ierr = VecPlaceArray(work1,(const PetscScalar*)&basis_mat[k*basis_dofs]);CHKERRQ(ierr); 239 ierr = VecSet(work1,one);CHKERRQ(ierr); 240 ierr = VecPlaceArray(work2,(const PetscScalar*)&Kbasis_mat[k*basis_dofs]);CHKERRQ(ierr); 241 ierr = MatMult(local_mat,work1,work2);CHKERRQ(ierr); 242 ierr = VecResetArray(work1);CHKERRQ(ierr); 243 ierr = VecResetArray(work2);CHKERRQ(ierr); 244 } 245 ierr = VecDestroy(&work1);CHKERRQ(ierr); 246 ierr = VecDestroy(&work2);CHKERRQ(ierr); 247 ierr = VecScatterDestroy(&scatter_ctx);CHKERRQ(ierr); 248 ierr = MatDenseRestoreArray(shell_ctx->basis_mat,&basis_mat);CHKERRQ(ierr); 249 ierr = MatDenseRestoreArray(shell_ctx->Kbasis_mat,&Kbasis_mat);CHKERRQ(ierr); 250 251 /* Assemble another Mat object in shell context */ 252 ierr = MatTransposeMatMult(shell_ctx->basis_mat,shell_ctx->Kbasis_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&small_mat);CHKERRQ(ierr); 253 ierr = MatFactorInfoInitialize(&matinfo);CHKERRQ(ierr); 254 ierr = ISCreateStride(PETSC_COMM_SELF,basis_size,0,1,&is_aux);CHKERRQ(ierr); 255 ierr = MatLUFactor(small_mat,is_aux,is_aux,&matinfo);CHKERRQ(ierr); 256 ierr = ISDestroy(&is_aux);CHKERRQ(ierr); 257 ierr = PetscMalloc1(basis_size*basis_size,&array_mat);CHKERRQ(ierr); 258 for (k=0;k<basis_size;k++) { 259 ierr = VecSet(shell_ctx->work_small_1,zero);CHKERRQ(ierr); 260 ierr = VecSetValue(shell_ctx->work_small_1,k,one,INSERT_VALUES);CHKERRQ(ierr); 261 ierr = VecAssemblyBegin(shell_ctx->work_small_1);CHKERRQ(ierr); 262 ierr = VecAssemblyEnd(shell_ctx->work_small_1);CHKERRQ(ierr); 263 ierr = MatSolve(small_mat,shell_ctx->work_small_1,shell_ctx->work_small_2);CHKERRQ(ierr); 264 ierr = VecGetArrayRead(shell_ctx->work_small_2,(const PetscScalar**)&array);CHKERRQ(ierr); 265 for (i=0;i<basis_size;i++) { 266 array_mat[i*basis_size+k]=array[i]; 267 } 268 ierr = VecRestoreArrayRead(shell_ctx->work_small_2,(const PetscScalar**)&array);CHKERRQ(ierr); 269 } 270 ierr = MatCreateSeqDense(PETSC_COMM_SELF,basis_size,basis_size,array_mat,&inv_small_mat);CHKERRQ(ierr); 271 ierr = MatMatMult(shell_ctx->basis_mat,inv_small_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&shell_ctx->Lbasis_mat);CHKERRQ(ierr); 272 ierr = PetscFree(array_mat);CHKERRQ(ierr); 273 ierr = MatDestroy(&inv_small_mat);CHKERRQ(ierr); 274 ierr = MatDestroy(&small_mat);CHKERRQ(ierr); 275 ierr = MatScale(shell_ctx->Kbasis_mat,m_one);CHKERRQ(ierr); 276 277 /* Rebuild local PC */ 278 ierr = KSPGetPC(local_ksp,&shell_ctx->local_pc);CHKERRQ(ierr); 279 ierr = PetscObjectReference((PetscObject)shell_ctx->local_pc);CHKERRQ(ierr); 280 ierr = PCCreate(PETSC_COMM_SELF,&newpc);CHKERRQ(ierr); 281 ierr = PCSetOperators(newpc,local_mat,local_mat);CHKERRQ(ierr); 282 ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr); 283 ierr = PCShellSetContext(newpc,shell_ctx);CHKERRQ(ierr); 284 ierr = PCShellSetApply(newpc,PCBDDCApplyNullSpaceCorrectionPC);CHKERRQ(ierr); 285 ierr = PCShellSetDestroy(newpc,PCBDDCDestroyNullSpaceCorrectionPC);CHKERRQ(ierr); 286 ierr = PCSetUp(newpc);CHKERRQ(ierr); 287 ierr = KSPSetPC(local_ksp,newpc);CHKERRQ(ierr); 288 ierr = PCDestroy(&newpc);CHKERRQ(ierr); 289 ierr = KSPSetUp(local_ksp);CHKERRQ(ierr); 290 } 291 /* test */ 292 if (pcbddc->dbg_flag && basis_dofs) { 293 KSP check_ksp; 294 PC check_pc; 295 Mat test_mat; 296 Vec work3; 297 PetscReal test_err,lambda_min,lambda_max; 298 PetscBool setsym,issym=PETSC_FALSE; 299 PetscInt tabs; 300 301 ierr = PetscViewerASCIIGetTab(pcbddc->dbg_viewer,&tabs);CHKERRQ(ierr); 302 ierr = KSPGetPC(local_ksp,&check_pc);CHKERRQ(ierr); 303 ierr = VecDuplicate(shell_ctx->work_full_1,&work1);CHKERRQ(ierr); 304 ierr = VecDuplicate(shell_ctx->work_full_1,&work2);CHKERRQ(ierr); 305 ierr = VecDuplicate(shell_ctx->work_full_1,&work3);CHKERRQ(ierr); 306 ierr = VecSetRandom(shell_ctx->work_small_1,NULL);CHKERRQ(ierr); 307 ierr = MatMult(shell_ctx->basis_mat,shell_ctx->work_small_1,work1);CHKERRQ(ierr); 308 ierr = VecCopy(work1,work2);CHKERRQ(ierr); 309 ierr = MatMult(local_mat,work1,work3);CHKERRQ(ierr); 310 ierr = PCApply(check_pc,work3,work1);CHKERRQ(ierr); 311 ierr = VecAXPY(work1,m_one,work2);CHKERRQ(ierr); 312 ierr = VecNorm(work1,NORM_INFINITY,&test_err);CHKERRQ(ierr); 313 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d error for nullspace correction for ",PetscGlobalRank);CHKERRQ(ierr); 314 ierr = PetscViewerASCIIUseTabs(pcbddc->dbg_viewer,PETSC_FALSE);CHKERRQ(ierr); 315 if (isdir) { 316 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Dirichlet ");CHKERRQ(ierr); 317 } else { 318 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Neumann ");CHKERRQ(ierr); 319 } 320 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"solver is :%1.14e\n",test_err);CHKERRQ(ierr); 321 ierr = PetscViewerASCIISetTab(pcbddc->dbg_viewer,tabs);CHKERRQ(ierr); 322 ierr = PetscViewerASCIIUseTabs(pcbddc->dbg_viewer,PETSC_TRUE);CHKERRQ(ierr); 323 324 ierr = MatTransposeMatMult(shell_ctx->Lbasis_mat,shell_ctx->Kbasis_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&test_mat);CHKERRQ(ierr); 325 ierr = MatShift(test_mat,one);CHKERRQ(ierr); 326 ierr = MatNorm(test_mat,NORM_INFINITY,&test_err);CHKERRQ(ierr); 327 ierr = MatDestroy(&test_mat);CHKERRQ(ierr); 328 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d error for nullspace matrices is :%1.14e\n",PetscGlobalRank,test_err);CHKERRQ(ierr); 329 330 /* Create ksp object suitable for extreme eigenvalues' estimation */ 331 ierr = KSPCreate(PETSC_COMM_SELF,&check_ksp);CHKERRQ(ierr); 332 ierr = KSPSetErrorIfNotConverged(check_ksp,pc->erroriffailure);CHKERRQ(ierr); 333 ierr = KSPSetOperators(check_ksp,local_mat,local_mat);CHKERRQ(ierr); 334 ierr = KSPSetTolerances(check_ksp,1.e-8,1.e-8,PETSC_DEFAULT,basis_dofs);CHKERRQ(ierr); 335 ierr = KSPSetComputeSingularValues(check_ksp,PETSC_TRUE);CHKERRQ(ierr); 336 ierr = MatIsSymmetricKnown(pc->pmat,&setsym,&issym);CHKERRQ(ierr); 337 if (issym) { 338 ierr = KSPSetType(check_ksp,KSPCG);CHKERRQ(ierr); 339 } 340 ierr = KSPSetPC(check_ksp,check_pc);CHKERRQ(ierr); 341 ierr = KSPSetUp(check_ksp);CHKERRQ(ierr); 342 ierr = VecSetRandom(work1,NULL);CHKERRQ(ierr); 343 ierr = MatMult(local_mat,work1,work2);CHKERRQ(ierr); 344 ierr = KSPSolve(check_ksp,work2,work2);CHKERRQ(ierr); 345 ierr = VecAXPY(work2,m_one,work1);CHKERRQ(ierr); 346 ierr = VecNorm(work2,NORM_INFINITY,&test_err);CHKERRQ(ierr); 347 ierr = KSPComputeExtremeSingularValues(check_ksp,&lambda_max,&lambda_min);CHKERRQ(ierr); 348 ierr = KSPGetIterationNumber(check_ksp,&k);CHKERRQ(ierr); 349 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d error for adapted KSP %1.14e (it %d, eigs %1.6e %1.6e)\n",PetscGlobalRank,test_err,k,lambda_min,lambda_max);CHKERRQ(ierr); 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 } 355 /* all processes shoud call this, even the void ones */ 356 if (pcbddc->dbg_flag) { 357 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 358 } 359 PetscFunctionReturn(0); 360 } 361 362 #undef __FUNCT__ 363 #define __FUNCT__ "PCBDDCNullSpaceAdaptGlobal" 364 PetscErrorCode PCBDDCNullSpaceAdaptGlobal(PC pc) 365 { 366 PC_IS* pcis = (PC_IS*)(pc->data); 367 PC_BDDC* pcbddc = (PC_BDDC*)(pc->data); 368 KSP inv_change; 369 const Vec *nsp_vecs; 370 Vec *new_nsp_vecs; 371 PetscInt i,nsp_size,new_nsp_size,start_new; 372 PetscBool nsp_has_cnst; 373 MatNullSpace new_nsp; 374 PetscErrorCode ierr; 375 376 PetscFunctionBegin; 377 /* create KSP for change of basis */ 378 ierr = MatGetSize(pcbddc->ChangeOfBasisMatrix,&i,NULL);CHKERRQ(ierr); 379 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&inv_change);CHKERRQ(ierr); 380 ierr = KSPSetErrorIfNotConverged(inv_change,pc->erroriffailure);CHKERRQ(ierr); 381 ierr = KSPSetOperators(inv_change,pcbddc->ChangeOfBasisMatrix,pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr); 382 ierr = KSPSetTolerances(inv_change,1.e-8,1.e-8,PETSC_DEFAULT,2*i);CHKERRQ(ierr); 383 if (pcbddc->dbg_flag) { 384 ierr = KSPMonitorSet(inv_change,KSPMonitorDefault,pcbddc->dbg_viewer,NULL);CHKERRQ(ierr); 385 } 386 ierr = KSPSetUp(inv_change);CHKERRQ(ierr); 387 388 /* get nullspace and transform it */ 389 ierr = MatNullSpaceGetVecs(pcbddc->NullSpace,&nsp_has_cnst,&nsp_size,&nsp_vecs);CHKERRQ(ierr); 390 new_nsp_size = nsp_size; 391 if (nsp_has_cnst) { 392 new_nsp_size++; 393 } 394 ierr = VecDuplicateVecs(pcis->vec1_global,new_nsp_size,&new_nsp_vecs);CHKERRQ(ierr); 395 396 start_new = 0; 397 if (nsp_has_cnst) { 398 start_new = 1; 399 ierr = VecSet(new_nsp_vecs[0],1.0);CHKERRQ(ierr); 400 if (pcbddc->dbg_flag) { 401 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 402 ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"Mapping constant in nullspace\n");CHKERRQ(ierr); 403 } 404 ierr = KSPSolve(inv_change,new_nsp_vecs[0],new_nsp_vecs[0]);CHKERRQ(ierr); 405 } 406 for (i=0;i<nsp_size;i++) { 407 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 408 ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"Mapping %dth vector in nullspace\n",i);CHKERRQ(ierr); 409 ierr = KSPSolve(inv_change,nsp_vecs[i],new_nsp_vecs[i+start_new]);CHKERRQ(ierr); 410 } 411 ierr = PCBDDCOrthonormalizeVecs(new_nsp_size,new_nsp_vecs);CHKERRQ(ierr); 412 ierr = MatNullSpaceCreate(PetscObjectComm((PetscObject)pc),PETSC_FALSE,new_nsp_size,new_nsp_vecs,&new_nsp);CHKERRQ(ierr); 413 ierr = PCBDDCSetNullSpace(pc,new_nsp);CHKERRQ(ierr); 414 415 /* free */ 416 ierr = KSPDestroy(&inv_change);CHKERRQ(ierr); 417 ierr = MatNullSpaceDestroy(&new_nsp);CHKERRQ(ierr); 418 ierr = VecDestroyVecs(new_nsp_size,&new_nsp_vecs);CHKERRQ(ierr); 419 420 /* check */ 421 if (pcbddc->dbg_flag) { 422 PetscBool nsp_t=PETSC_FALSE; 423 Mat temp_mat; 424 Mat_IS* matis = (Mat_IS*)pc->pmat->data; 425 426 temp_mat = matis->A; 427 matis->A = pcbddc->local_mat; 428 pcbddc->local_mat = temp_mat; 429 ierr = MatNullSpaceTest(pcbddc->NullSpace,pc->pmat,&nsp_t);CHKERRQ(ierr); 430 ierr = PetscPrintf(PetscObjectComm((PetscObject)(pc->pmat)),"Check nullspace with change of basis: %d\n",nsp_t);CHKERRQ(ierr); 431 temp_mat = matis->A; 432 matis->A = pcbddc->local_mat; 433 pcbddc->local_mat = temp_mat; 434 } 435 PetscFunctionReturn(0); 436 } 437