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