1 #include "bddc.h" 2 #include "bddcprivate.h" 3 #include <petscblaslapack.h> 4 5 #undef __FUNCT__ 6 #define __FUNCT__ "PCBDDCResetCustomization" 7 PetscErrorCode PCBDDCResetCustomization(PC pc) 8 { 9 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 10 PetscInt i; 11 PetscErrorCode ierr; 12 13 PetscFunctionBegin; 14 ierr = PCBDDCGraphResetCSR(pcbddc->mat_graph);CHKERRQ(ierr); 15 ierr = ISDestroy(&pcbddc->user_primal_vertices);CHKERRQ(ierr); 16 ierr = MatNullSpaceDestroy(&pcbddc->NullSpace);CHKERRQ(ierr); 17 ierr = ISDestroy(&pcbddc->NeumannBoundaries);CHKERRQ(ierr); 18 ierr = ISDestroy(&pcbddc->DirichletBoundaries);CHKERRQ(ierr); 19 for (i=0;i<pcbddc->n_ISForDofs;i++) { 20 ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr); 21 } 22 ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr); 23 PetscFunctionReturn(0); 24 } 25 26 #undef __FUNCT__ 27 #define __FUNCT__ "PCBDDCResetTopography" 28 PetscErrorCode PCBDDCResetTopography(PC pc) 29 { 30 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 31 PetscErrorCode ierr; 32 33 PetscFunctionBegin; 34 ierr = MatDestroy(&pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr); 35 ierr = MatDestroy(&pcbddc->ConstraintMatrix);CHKERRQ(ierr); 36 ierr = PCBDDCGraphReset(pcbddc->mat_graph);CHKERRQ(ierr); 37 PetscFunctionReturn(0); 38 } 39 40 #undef __FUNCT__ 41 #define __FUNCT__ "PCBDDCResetSolvers" 42 PetscErrorCode PCBDDCResetSolvers(PC pc) 43 { 44 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 45 PetscErrorCode ierr; 46 47 PetscFunctionBegin; 48 ierr = VecDestroy(&pcbddc->temp_solution);CHKERRQ(ierr); 49 ierr = VecDestroy(&pcbddc->original_rhs);CHKERRQ(ierr); 50 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 51 ierr = VecDestroy(&pcbddc->coarse_vec);CHKERRQ(ierr); 52 ierr = VecDestroy(&pcbddc->coarse_rhs);CHKERRQ(ierr); 53 ierr = KSPDestroy(&pcbddc->coarse_ksp);CHKERRQ(ierr); 54 ierr = MatDestroy(&pcbddc->coarse_mat);CHKERRQ(ierr); 55 ierr = MatDestroy(&pcbddc->coarse_phi_B);CHKERRQ(ierr); 56 ierr = MatDestroy(&pcbddc->coarse_phi_D);CHKERRQ(ierr); 57 ierr = MatDestroy(&pcbddc->coarse_psi_B);CHKERRQ(ierr); 58 ierr = MatDestroy(&pcbddc->coarse_psi_D);CHKERRQ(ierr); 59 ierr = VecDestroy(&pcbddc->vec1_P);CHKERRQ(ierr); 60 ierr = VecDestroy(&pcbddc->vec1_C);CHKERRQ(ierr); 61 ierr = MatDestroy(&pcbddc->local_auxmat1);CHKERRQ(ierr); 62 ierr = MatDestroy(&pcbddc->local_auxmat2);CHKERRQ(ierr); 63 ierr = VecDestroy(&pcbddc->vec1_R);CHKERRQ(ierr); 64 ierr = VecDestroy(&pcbddc->vec2_R);CHKERRQ(ierr); 65 ierr = VecDestroy(&pcbddc->vec4_D);CHKERRQ(ierr); 66 ierr = VecScatterDestroy(&pcbddc->R_to_B);CHKERRQ(ierr); 67 ierr = VecScatterDestroy(&pcbddc->R_to_D);CHKERRQ(ierr); 68 ierr = VecScatterDestroy(&pcbddc->coarse_loc_to_glob);CHKERRQ(ierr); 69 ierr = PetscFree(pcbddc->local_primal_indices);CHKERRQ(ierr); 70 ierr = PetscFree(pcbddc->replicated_local_primal_indices);CHKERRQ(ierr); 71 ierr = PetscFree(pcbddc->replicated_local_primal_values);CHKERRQ(ierr); 72 ierr = PetscFree(pcbddc->local_primal_displacements);CHKERRQ(ierr); 73 ierr = PetscFree(pcbddc->local_primal_sizes);CHKERRQ(ierr); 74 PetscFunctionReturn(0); 75 } 76 77 #undef __FUNCT__ 78 #define __FUNCT__ "PCBDDCSetUseExactDirichlet" 79 PetscErrorCode PCBDDCSetUseExactDirichlet(PC pc,PetscBool use) 80 { 81 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 82 83 PetscFunctionBegin; 84 pcbddc->use_exact_dirichlet=use; 85 PetscFunctionReturn(0); 86 } 87 88 #undef __FUNCT__ 89 #define __FUNCT__ "PCBDDCSetUpLocalSolvers" 90 PetscErrorCode PCBDDCSetUpLocalSolvers(PC pc, IS is_I_local, IS is_R_local) 91 { 92 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 93 PC_IS *pcis = (PC_IS*)pc->data; 94 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 95 PC pc_temp; 96 Mat A_RR; 97 Vec vec1,vec2,vec3; 98 MatStructure matstruct; 99 PetscScalar m_one = -1.0; 100 PetscReal value; 101 PetscInt n_D,n_R,use_exact,use_exact_reduced; 102 PetscErrorCode ierr; 103 104 PetscFunctionBegin; 105 /* Creating PC contexts for local Dirichlet and Neumann problems */ 106 ierr = PCGetOperators(pc,NULL,NULL,&matstruct);CHKERRQ(ierr); 107 108 /* when matstruct is SAME_PRECONDITIONER, we shouldn't be here */ 109 110 /* DIRICHLET PROBLEM */ 111 /* Matrix for Dirichlet problem is A_II */ 112 /* HACK: A_II can be changed between nonlinear iterations */ 113 if (pc->setupcalled) { /* we dont need to rebuild dirichlet problem the first time we build BDDC */ 114 if (matstruct == SAME_PRECONDITIONER) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"This should not happen"); 115 if (matstruct == SAME_NONZERO_PATTERN) { 116 ierr = MatGetSubMatrix(matis->A,is_I_local,is_I_local,MAT_REUSE_MATRIX,&pcis->A_II);CHKERRQ(ierr); 117 } else { 118 ierr = MatDestroy(&pcis->A_II);CHKERRQ(ierr); 119 ierr = MatGetSubMatrix(matis->A,is_I_local,is_I_local,MAT_INITIAL_MATRIX,&pcis->A_II);CHKERRQ(ierr); 120 } 121 } 122 ierr = ISGetSize(is_I_local,&n_D);CHKERRQ(ierr); 123 if (!pcbddc->ksp_D) { /* create object if not yet build */ 124 ierr = KSPCreate(PETSC_COMM_SELF,&pcbddc->ksp_D);CHKERRQ(ierr); 125 ierr = PetscObjectIncrementTabLevel((PetscObject)pcbddc->ksp_D,(PetscObject)pc,1);CHKERRQ(ierr); 126 /* default */ 127 ierr = KSPSetType(pcbddc->ksp_D,KSPPREONLY);CHKERRQ(ierr); 128 ierr = KSPSetOptionsPrefix(pcbddc->ksp_D,"dirichlet_");CHKERRQ(ierr); 129 ierr = KSPGetPC(pcbddc->ksp_D,&pc_temp);CHKERRQ(ierr); 130 ierr = PCSetType(pc_temp,PCLU);CHKERRQ(ierr); 131 ierr = PCFactorSetReuseFill(pc_temp,PETSC_TRUE);CHKERRQ(ierr); 132 } 133 ierr = KSPSetOperators(pcbddc->ksp_D,pcis->A_II,pcis->A_II,matstruct);CHKERRQ(ierr); 134 /* Allow user's customization */ 135 ierr = KSPSetFromOptions(pcbddc->ksp_D);CHKERRQ(ierr); 136 /* umfpack interface has a bug when matrix dimension is zero. TODO solve from umfpack interface */ 137 if (!n_D) { 138 ierr = KSPGetPC(pcbddc->ksp_D,&pc_temp);CHKERRQ(ierr); 139 ierr = PCSetType(pc_temp,PCNONE);CHKERRQ(ierr); 140 } 141 /* Set Up KSP for Dirichlet problem of BDDC */ 142 ierr = KSPSetUp(pcbddc->ksp_D);CHKERRQ(ierr); 143 /* set ksp_D into pcis data */ 144 ierr = KSPDestroy(&pcis->ksp_D);CHKERRQ(ierr); 145 ierr = PetscObjectReference((PetscObject)pcbddc->ksp_D);CHKERRQ(ierr); 146 pcis->ksp_D = pcbddc->ksp_D; 147 148 /* NEUMANN PROBLEM */ 149 /* Matrix for Neumann problem is A_RR -> we need to create it */ 150 ierr = ISGetSize(is_R_local,&n_R);CHKERRQ(ierr); 151 ierr = MatGetSubMatrix(pcbddc->local_mat,is_R_local,is_R_local,MAT_INITIAL_MATRIX,&A_RR);CHKERRQ(ierr); 152 if (!pcbddc->ksp_R) { /* create object if not yet build */ 153 ierr = KSPCreate(PETSC_COMM_SELF,&pcbddc->ksp_R);CHKERRQ(ierr); 154 ierr = PetscObjectIncrementTabLevel((PetscObject)pcbddc->ksp_R,(PetscObject)pc,1);CHKERRQ(ierr); 155 /* default */ 156 ierr = KSPSetType(pcbddc->ksp_R,KSPPREONLY);CHKERRQ(ierr); 157 ierr = KSPSetOptionsPrefix(pcbddc->ksp_R,"neumann_");CHKERRQ(ierr); 158 ierr = KSPGetPC(pcbddc->ksp_R,&pc_temp);CHKERRQ(ierr); 159 ierr = PCSetType(pc_temp,PCLU);CHKERRQ(ierr); 160 ierr = PCFactorSetReuseFill(pc_temp,PETSC_TRUE);CHKERRQ(ierr); 161 } 162 ierr = KSPSetOperators(pcbddc->ksp_R,A_RR,A_RR,matstruct);CHKERRQ(ierr); 163 /* Allow user's customization */ 164 ierr = KSPSetFromOptions(pcbddc->ksp_R);CHKERRQ(ierr); 165 /* umfpack interface has a bug when matrix dimension is zero. TODO solve from umfpack interface */ 166 if (!n_R) { 167 ierr = KSPGetPC(pcbddc->ksp_R,&pc_temp);CHKERRQ(ierr); 168 ierr = PCSetType(pc_temp,PCNONE);CHKERRQ(ierr); 169 } 170 /* Set Up KSP for Neumann problem of BDDC */ 171 ierr = KSPSetUp(pcbddc->ksp_R);CHKERRQ(ierr); 172 173 /* check Dirichlet and Neumann solvers and adapt them if a nullspace correction is needed */ 174 175 /* Dirichlet */ 176 ierr = MatGetVecs(pcis->A_II,&vec1,&vec2);CHKERRQ(ierr); 177 ierr = VecDuplicate(vec1,&vec3);CHKERRQ(ierr); 178 ierr = VecSetRandom(vec1,NULL);CHKERRQ(ierr); 179 ierr = MatMult(pcis->A_II,vec1,vec2);CHKERRQ(ierr); 180 ierr = KSPSolve(pcbddc->ksp_D,vec2,vec3);CHKERRQ(ierr); 181 ierr = VecAXPY(vec3,m_one,vec1);CHKERRQ(ierr); 182 ierr = VecNorm(vec3,NORM_INFINITY,&value);CHKERRQ(ierr); 183 ierr = VecDestroy(&vec1);CHKERRQ(ierr); 184 ierr = VecDestroy(&vec2);CHKERRQ(ierr); 185 ierr = VecDestroy(&vec3);CHKERRQ(ierr); 186 /* need to be adapted? */ 187 use_exact = (PetscAbsReal(value) > 1.e-4 ? 0 : 1); 188 ierr = MPI_Allreduce(&use_exact,&use_exact_reduced,1,MPIU_INT,MPI_LAND,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr); 189 ierr = PCBDDCSetUseExactDirichlet(pc,(PetscBool)use_exact_reduced);CHKERRQ(ierr); 190 /* print info */ 191 if (pcbddc->dbg_flag) { 192 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 193 ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"--------------------------------------------------\n");CHKERRQ(ierr); 194 ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"Checking solution of Dirichlet and Neumann problems\n");CHKERRQ(ierr); 195 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Dirichlet solve = % 1.14e \n",PetscGlobalRank,value);CHKERRQ(ierr); 196 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 197 } 198 if (n_D && pcbddc->NullSpace && !use_exact_reduced && !pcbddc->inexact_prec_type) { 199 ierr = PCBDDCNullSpaceAssembleCorrection(pc,is_I_local);CHKERRQ(ierr); 200 } 201 202 /* Neumann */ 203 ierr = MatGetVecs(A_RR,&vec1,&vec2);CHKERRQ(ierr); 204 ierr = VecDuplicate(vec1,&vec3);CHKERRQ(ierr); 205 ierr = VecSetRandom(vec1,NULL);CHKERRQ(ierr); 206 ierr = MatMult(A_RR,vec1,vec2);CHKERRQ(ierr); 207 ierr = KSPSolve(pcbddc->ksp_R,vec2,vec3);CHKERRQ(ierr); 208 ierr = VecAXPY(vec3,m_one,vec1);CHKERRQ(ierr); 209 ierr = VecNorm(vec3,NORM_INFINITY,&value);CHKERRQ(ierr); 210 ierr = VecDestroy(&vec1);CHKERRQ(ierr); 211 ierr = VecDestroy(&vec2);CHKERRQ(ierr); 212 ierr = VecDestroy(&vec3);CHKERRQ(ierr); 213 /* need to be adapted? */ 214 use_exact = (PetscAbsReal(value) > 1.e-4 ? 0 : 1); 215 if (PetscAbsReal(value) > 1.e-4) use_exact = 0; 216 ierr = MPI_Allreduce(&use_exact,&use_exact_reduced,1,MPIU_INT,MPI_LAND,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr); 217 /* print info */ 218 if (pcbddc->dbg_flag) { 219 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Neumann solve = % 1.14e \n",PetscGlobalRank,value);CHKERRQ(ierr); 220 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 221 } 222 if (n_R && pcbddc->NullSpace && !use_exact_reduced) { /* is it the right logic? */ 223 ierr = PCBDDCNullSpaceAssembleCorrection(pc,is_R_local);CHKERRQ(ierr); 224 } 225 226 /* free Neumann problem's matrix */ 227 ierr = MatDestroy(&A_RR);CHKERRQ(ierr); 228 PetscFunctionReturn(0); 229 } 230 231 #undef __FUNCT__ 232 #define __FUNCT__ "PCBDDCSolveSaddlePoint" 233 static PetscErrorCode PCBDDCSolveSaddlePoint(PC pc) 234 { 235 PetscErrorCode ierr; 236 PC_BDDC* pcbddc = (PC_BDDC*)(pc->data); 237 238 PetscFunctionBegin; 239 ierr = KSPSolve(pcbddc->ksp_R,pcbddc->vec1_R,pcbddc->vec2_R);CHKERRQ(ierr); 240 if (pcbddc->local_auxmat1) { 241 ierr = MatMult(pcbddc->local_auxmat1,pcbddc->vec2_R,pcbddc->vec1_C);CHKERRQ(ierr); 242 ierr = MatMultAdd(pcbddc->local_auxmat2,pcbddc->vec1_C,pcbddc->vec2_R,pcbddc->vec2_R);CHKERRQ(ierr); 243 } 244 PetscFunctionReturn(0); 245 } 246 247 #undef __FUNCT__ 248 #define __FUNCT__ "PCBDDCApplyInterfacePreconditioner" 249 PetscErrorCode PCBDDCApplyInterfacePreconditioner(PC pc) 250 { 251 PetscErrorCode ierr; 252 PC_BDDC* pcbddc = (PC_BDDC*)(pc->data); 253 PC_IS* pcis = (PC_IS*) (pc->data); 254 const PetscScalar zero = 0.0; 255 256 PetscFunctionBegin; 257 /* Application of PHI^T (or PSI^T) */ 258 if (pcbddc->coarse_psi_B) { 259 ierr = MatMultTranspose(pcbddc->coarse_psi_B,pcis->vec1_B,pcbddc->vec1_P);CHKERRQ(ierr); 260 if (pcbddc->inexact_prec_type) { ierr = MatMultTransposeAdd(pcbddc->coarse_psi_D,pcis->vec1_D,pcbddc->vec1_P,pcbddc->vec1_P);CHKERRQ(ierr); } 261 } else { 262 ierr = MatMultTranspose(pcbddc->coarse_phi_B,pcis->vec1_B,pcbddc->vec1_P);CHKERRQ(ierr); 263 if (pcbddc->inexact_prec_type) { ierr = MatMultTransposeAdd(pcbddc->coarse_phi_D,pcis->vec1_D,pcbddc->vec1_P,pcbddc->vec1_P);CHKERRQ(ierr); } 264 } 265 /* Scatter data of coarse_rhs */ 266 if (pcbddc->coarse_rhs) { ierr = VecSet(pcbddc->coarse_rhs,zero);CHKERRQ(ierr); } 267 ierr = PCBDDCScatterCoarseDataBegin(pc,pcbddc->vec1_P,pcbddc->coarse_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 268 269 /* Local solution on R nodes */ 270 ierr = VecSet(pcbddc->vec1_R,zero);CHKERRQ(ierr); 271 ierr = VecScatterBegin(pcbddc->R_to_B,pcis->vec1_B,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 272 ierr = VecScatterEnd (pcbddc->R_to_B,pcis->vec1_B,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 273 if (pcbddc->inexact_prec_type) { 274 ierr = VecScatterBegin(pcbddc->R_to_D,pcis->vec1_D,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 275 ierr = VecScatterEnd (pcbddc->R_to_D,pcis->vec1_D,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 276 } 277 ierr = PCBDDCSolveSaddlePoint(pc);CHKERRQ(ierr); 278 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 279 ierr = VecScatterBegin(pcbddc->R_to_B,pcbddc->vec2_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 280 ierr = VecScatterEnd (pcbddc->R_to_B,pcbddc->vec2_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 281 if (pcbddc->inexact_prec_type) { 282 ierr = VecScatterBegin(pcbddc->R_to_D,pcbddc->vec2_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 283 ierr = VecScatterEnd (pcbddc->R_to_D,pcbddc->vec2_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 284 } 285 286 /* Coarse solution */ 287 ierr = PCBDDCScatterCoarseDataEnd(pc,pcbddc->vec1_P,pcbddc->coarse_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 288 if (pcbddc->coarse_rhs) { /* TODO remove null space when doing multilevel */ 289 ierr = KSPSolve(pcbddc->coarse_ksp,pcbddc->coarse_rhs,pcbddc->coarse_vec);CHKERRQ(ierr); 290 } 291 ierr = PCBDDCScatterCoarseDataBegin(pc,pcbddc->coarse_vec,pcbddc->vec1_P,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 292 ierr = PCBDDCScatterCoarseDataEnd (pc,pcbddc->coarse_vec,pcbddc->vec1_P,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 293 294 /* Sum contributions from two levels */ 295 ierr = MatMultAdd(pcbddc->coarse_phi_B,pcbddc->vec1_P,pcis->vec1_B,pcis->vec1_B);CHKERRQ(ierr); 296 if (pcbddc->inexact_prec_type) { ierr = MatMultAdd(pcbddc->coarse_phi_D,pcbddc->vec1_P,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); } 297 PetscFunctionReturn(0); 298 } 299 300 #undef __FUNCT__ 301 #define __FUNCT__ "PCBDDCScatterCoarseDataBegin" 302 PetscErrorCode PCBDDCScatterCoarseDataBegin(PC pc,Vec vec_from, Vec vec_to, InsertMode imode, ScatterMode smode) 303 { 304 PetscErrorCode ierr; 305 PC_BDDC* pcbddc = (PC_BDDC*)(pc->data); 306 307 PetscFunctionBegin; 308 switch (pcbddc->coarse_communications_type) { 309 case SCATTERS_BDDC: 310 ierr = VecScatterBegin(pcbddc->coarse_loc_to_glob,vec_from,vec_to,imode,smode);CHKERRQ(ierr); 311 break; 312 case GATHERS_BDDC: 313 break; 314 } 315 PetscFunctionReturn(0); 316 } 317 318 #undef __FUNCT__ 319 #define __FUNCT__ "PCBDDCScatterCoarseDataEnd" 320 PetscErrorCode PCBDDCScatterCoarseDataEnd(PC pc,Vec vec_from, Vec vec_to, InsertMode imode, ScatterMode smode) 321 { 322 PetscErrorCode ierr; 323 PC_BDDC* pcbddc = (PC_BDDC*)(pc->data); 324 PetscScalar* array_to; 325 PetscScalar* array_from; 326 MPI_Comm comm; 327 PetscInt i; 328 329 PetscFunctionBegin; 330 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 331 switch (pcbddc->coarse_communications_type) { 332 case SCATTERS_BDDC: 333 ierr = VecScatterEnd(pcbddc->coarse_loc_to_glob,vec_from,vec_to,imode,smode);CHKERRQ(ierr); 334 break; 335 case GATHERS_BDDC: 336 if (vec_from) { 337 ierr = VecGetArray(vec_from,&array_from);CHKERRQ(ierr); 338 } 339 if (vec_to) { 340 ierr = VecGetArray(vec_to,&array_to);CHKERRQ(ierr); 341 } 342 switch(pcbddc->coarse_problem_type){ 343 case SEQUENTIAL_BDDC: 344 if (smode == SCATTER_FORWARD) { 345 ierr = MPI_Gatherv(&array_from[0],pcbddc->local_primal_size,MPIU_SCALAR,&pcbddc->replicated_local_primal_values[0],pcbddc->local_primal_sizes,pcbddc->local_primal_displacements,MPIU_SCALAR,0,comm);CHKERRQ(ierr); 346 if (vec_to) { 347 if (imode == ADD_VALUES) { 348 for (i=0;i<pcbddc->replicated_primal_size;i++) { 349 array_to[pcbddc->replicated_local_primal_indices[i]]+=pcbddc->replicated_local_primal_values[i]; 350 } 351 } else { 352 for (i=0;i<pcbddc->replicated_primal_size;i++) { 353 array_to[pcbddc->replicated_local_primal_indices[i]]=pcbddc->replicated_local_primal_values[i]; 354 } 355 } 356 } 357 } else { 358 if (vec_from) { 359 if (imode == ADD_VALUES) { 360 MPI_Comm vec_from_comm; 361 ierr = PetscObjectGetComm((PetscObject)(vec_from),&vec_from_comm);CHKERRQ(ierr); 362 SETERRQ2(vec_from_comm,PETSC_ERR_SUP,"Unsupported insert mode ADD_VALUES for SCATTER_REVERSE in %s for case %d\n",__FUNCT__,pcbddc->coarse_problem_type); 363 } 364 for (i=0;i<pcbddc->replicated_primal_size;i++) { 365 pcbddc->replicated_local_primal_values[i]=array_from[pcbddc->replicated_local_primal_indices[i]]; 366 } 367 } 368 ierr = MPI_Scatterv(&pcbddc->replicated_local_primal_values[0],pcbddc->local_primal_sizes,pcbddc->local_primal_displacements,MPIU_SCALAR,&array_to[0],pcbddc->local_primal_size,MPIU_SCALAR,0,comm);CHKERRQ(ierr); 369 } 370 break; 371 case REPLICATED_BDDC: 372 if (smode == SCATTER_FORWARD) { 373 ierr = MPI_Allgatherv(&array_from[0],pcbddc->local_primal_size,MPIU_SCALAR,&pcbddc->replicated_local_primal_values[0],pcbddc->local_primal_sizes,pcbddc->local_primal_displacements,MPIU_SCALAR,comm);CHKERRQ(ierr); 374 if (imode == ADD_VALUES) { 375 for (i=0;i<pcbddc->replicated_primal_size;i++) { 376 array_to[pcbddc->replicated_local_primal_indices[i]]+=pcbddc->replicated_local_primal_values[i]; 377 } 378 } else { 379 for (i=0;i<pcbddc->replicated_primal_size;i++) { 380 array_to[pcbddc->replicated_local_primal_indices[i]]=pcbddc->replicated_local_primal_values[i]; 381 } 382 } 383 } else { /* no communications needed for SCATTER_REVERSE since needed data is already present */ 384 if (imode == ADD_VALUES) { 385 for (i=0;i<pcbddc->local_primal_size;i++) { 386 array_to[i]+=array_from[pcbddc->local_primal_indices[i]]; 387 } 388 } else { 389 for (i=0;i<pcbddc->local_primal_size;i++) { 390 array_to[i]=array_from[pcbddc->local_primal_indices[i]]; 391 } 392 } 393 } 394 break; 395 case MULTILEVEL_BDDC: 396 break; 397 case PARALLEL_BDDC: 398 break; 399 } 400 if (vec_from) { 401 ierr = VecRestoreArray(vec_from,&array_from);CHKERRQ(ierr); 402 } 403 if (vec_to) { 404 ierr = VecRestoreArray(vec_to,&array_to);CHKERRQ(ierr); 405 } 406 break; 407 } 408 PetscFunctionReturn(0); 409 } 410 411 #undef __FUNCT__ 412 #define __FUNCT__ "PCBDDCConstraintsSetUp" 413 PetscErrorCode PCBDDCConstraintsSetUp(PC pc) 414 { 415 PetscErrorCode ierr; 416 PC_IS* pcis = (PC_IS*)(pc->data); 417 PC_BDDC* pcbddc = (PC_BDDC*)pc->data; 418 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 419 PetscInt *nnz,*is_indices; 420 PetscScalar *temp_quadrature_constraint; 421 PetscInt *temp_indices,*temp_indices_to_constraint,*temp_indices_to_constraint_B,*local_to_B; 422 PetscInt local_primal_size,i,j,k,total_counts,max_size_of_constraint; 423 PetscInt n_vertices,size_of_constraint; 424 PetscReal real_value; 425 PetscBool nnsp_has_cnst=PETSC_FALSE,use_nnsp_true=pcbddc->use_nnsp_true; 426 PetscInt nnsp_size=0,nnsp_addone=0,temp_constraints,temp_start_ptr,n_ISForFaces,n_ISForEdges; 427 IS *used_IS,ISForVertices,*ISForFaces,*ISForEdges; 428 MatType impMatType=MATSEQAIJ; 429 PetscBLASInt Bs,Bt,lwork,lierr; 430 PetscReal tol=1.0e-8; 431 MatNullSpace nearnullsp; 432 const Vec *nearnullvecs; 433 Vec *localnearnullsp; 434 PetscScalar *work,*temp_basis,*array_vector,*correlation_mat; 435 PetscReal *rwork,*singular_vals; 436 PetscBLASInt Bone=1,*ipiv; 437 Vec temp_vec; 438 Mat temp_mat; 439 KSP temp_ksp; 440 PC temp_pc; 441 PetscInt s,start_constraint,dual_dofs; 442 PetscBool compute_submatrix,useksp=PETSC_FALSE; 443 PetscInt *aux_primal_permutation,*aux_primal_numbering; 444 PetscBool boolforchange,*change_basis; 445 /* some ugly conditional declarations */ 446 #if defined(PETSC_MISSING_LAPACK_GESVD) 447 PetscScalar one=1.0,zero=0.0; 448 PetscInt ii; 449 PetscScalar *singular_vectors; 450 PetscBLASInt *iwork,*ifail; 451 PetscReal dummy_real,abs_tol; 452 PetscBLASInt eigs_found; 453 #endif 454 PetscBLASInt dummy_int; 455 PetscScalar dummy_scalar; 456 PetscBool used_vertex,get_faces,get_edges,get_vertices; 457 458 PetscFunctionBegin; 459 /* Get index sets for faces, edges and vertices from graph */ 460 get_faces = PETSC_TRUE; 461 get_edges = PETSC_TRUE; 462 get_vertices = PETSC_TRUE; 463 if (pcbddc->vertices_flag) { 464 get_faces = PETSC_FALSE; 465 get_edges = PETSC_FALSE; 466 } 467 if (pcbddc->constraints_flag) { 468 get_vertices = PETSC_FALSE; 469 } 470 if (pcbddc->faces_flag) { 471 get_edges = PETSC_FALSE; 472 } 473 if (pcbddc->edges_flag) { 474 get_faces = PETSC_FALSE; 475 } 476 /* default */ 477 if (!get_faces && !get_edges && !get_vertices) { 478 get_vertices = PETSC_TRUE; 479 } 480 ierr = PCBDDCGraphGetCandidatesIS(pcbddc->mat_graph,get_faces,get_edges,get_vertices,&n_ISForFaces,&ISForFaces,&n_ISForEdges,&ISForEdges,&ISForVertices); 481 if (pcbddc->dbg_flag) { 482 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"--------------------------------------------------------------\n");CHKERRQ(ierr); 483 i = 0; 484 if (ISForVertices) { 485 ierr = ISGetSize(ISForVertices,&i);CHKERRQ(ierr); 486 } 487 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d got %02d local candidate vertices\n",PetscGlobalRank,i);CHKERRQ(ierr); 488 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d got %02d local candidate edges\n",PetscGlobalRank,n_ISForEdges);CHKERRQ(ierr); 489 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d got %02d local candidate faces\n",PetscGlobalRank,n_ISForFaces);CHKERRQ(ierr); 490 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 491 } 492 /* check if near null space is attached to global mat */ 493 ierr = MatGetNearNullSpace(pc->pmat,&nearnullsp);CHKERRQ(ierr); 494 if (nearnullsp) { 495 ierr = MatNullSpaceGetVecs(nearnullsp,&nnsp_has_cnst,&nnsp_size,&nearnullvecs);CHKERRQ(ierr); 496 } else { /* if near null space is not provided it uses constants */ 497 nnsp_has_cnst = PETSC_TRUE; 498 use_nnsp_true = PETSC_TRUE; 499 } 500 if (nnsp_has_cnst) { 501 nnsp_addone = 1; 502 } 503 /* 504 Evaluate maximum storage size needed by the procedure 505 - temp_indices will contain start index of each constraint stored as follows 506 - temp_indices_to_constraint [temp_indices[i],...,temp[indices[i+1]-1] will contain the indices (in local numbering) on which the constraint acts 507 - temp_indices_to_constraint_B[temp_indices[i],...,temp[indices[i+1]-1] will contain the indices (in boundary numbering) on which the constraint acts 508 - temp_quadrature_constraint [temp_indices[i],...,temp[indices[i+1]-1] will contain the scalars representing the constraint itself 509 */ 510 total_counts = n_ISForFaces+n_ISForEdges; 511 total_counts *= (nnsp_addone+nnsp_size); 512 n_vertices = 0; 513 if (ISForVertices) { 514 ierr = ISGetSize(ISForVertices,&n_vertices);CHKERRQ(ierr); 515 } 516 total_counts += n_vertices; 517 ierr = PetscMalloc((total_counts+1)*sizeof(PetscInt),&temp_indices);CHKERRQ(ierr); 518 ierr = PetscMalloc((total_counts+1)*sizeof(PetscBool),&change_basis);CHKERRQ(ierr); 519 total_counts = 0; 520 max_size_of_constraint = 0; 521 for (i=0;i<n_ISForEdges+n_ISForFaces;i++) { 522 if (i<n_ISForEdges) { 523 used_IS = &ISForEdges[i]; 524 } else { 525 used_IS = &ISForFaces[i-n_ISForEdges]; 526 } 527 ierr = ISGetSize(*used_IS,&j);CHKERRQ(ierr); 528 total_counts += j; 529 max_size_of_constraint = PetscMax(j,max_size_of_constraint); 530 } 531 total_counts *= (nnsp_addone+nnsp_size); 532 total_counts += n_vertices; 533 ierr = PetscMalloc(total_counts*sizeof(PetscScalar),&temp_quadrature_constraint);CHKERRQ(ierr); 534 ierr = PetscMalloc(total_counts*sizeof(PetscInt),&temp_indices_to_constraint);CHKERRQ(ierr); 535 ierr = PetscMalloc(total_counts*sizeof(PetscInt),&temp_indices_to_constraint_B);CHKERRQ(ierr); 536 ierr = PetscMalloc(pcis->n*sizeof(PetscInt),&local_to_B);CHKERRQ(ierr); 537 ierr = ISGetIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 538 for (i=0;i<pcis->n;i++) { 539 local_to_B[i]=-1; 540 } 541 for (i=0;i<pcis->n_B;i++) { 542 local_to_B[is_indices[i]]=i; 543 } 544 ierr = ISRestoreIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 545 546 /* First we issue queries to allocate optimal workspace for LAPACKgesvd or LAPACKsyev/LAPACKheev */ 547 rwork = 0; 548 work = 0; 549 singular_vals = 0; 550 temp_basis = 0; 551 correlation_mat = 0; 552 if (!pcbddc->use_nnsp_true) { 553 PetscScalar temp_work; 554 #if defined(PETSC_MISSING_LAPACK_GESVD) 555 /* POD */ 556 PetscInt max_n; 557 max_n = nnsp_addone+nnsp_size; 558 /* using some techniques borrowed from Proper Orthogonal Decomposition */ 559 ierr = PetscMalloc(max_n*max_n*sizeof(PetscScalar),&correlation_mat);CHKERRQ(ierr); 560 ierr = PetscMalloc(max_n*max_n*sizeof(PetscScalar),&singular_vectors);CHKERRQ(ierr); 561 ierr = PetscMalloc(max_n*sizeof(PetscReal),&singular_vals);CHKERRQ(ierr); 562 ierr = PetscMalloc(max_size_of_constraint*(nnsp_addone+nnsp_size)*sizeof(PetscScalar),&temp_basis);CHKERRQ(ierr); 563 #if defined(PETSC_USE_COMPLEX) 564 ierr = PetscMalloc(3*max_n*sizeof(PetscReal),&rwork);CHKERRQ(ierr); 565 #endif 566 ierr = PetscMalloc(5*max_n*sizeof(PetscBLASInt),&iwork);CHKERRQ(ierr); 567 ierr = PetscMalloc(max_n*sizeof(PetscBLASInt),&ifail);CHKERRQ(ierr); 568 /* now we evaluate the optimal workspace using query with lwork=-1 */ 569 ierr = PetscBLASIntCast(max_n,&Bt);CHKERRQ(ierr); 570 lwork=-1; 571 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 572 #if !defined(PETSC_USE_COMPLEX) 573 abs_tol=1.e-8; 574 /* LAPACKsyev_("V","U",&Bt,correlation_mat,&Bt,singular_vals,&temp_work,&lwork,&lierr); */ 575 PetscStackCallBLAS("LAPACKsyevx",LAPACKsyevx_("V","A","U",&Bt,correlation_mat,&Bt,&dummy_real,&dummy_real,&dummy_int,&dummy_int,&abs_tol,&eigs_found,singular_vals,singular_vectors,&Bt,&temp_work,&lwork,iwork,ifail,&lierr)); 576 #else 577 /* LAPACKsyev_("V","U",&Bt,correlation_mat,&Bt,singular_vals,&temp_work,&lwork,rwork,&lierr); */ 578 /* LAPACK call is missing here! TODO */ 579 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not yet implemented for complexes when PETSC_MISSING_GESVD = 1"); 580 #endif 581 if ( lierr ) { 582 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in query to SYEVX Lapack routine %d",(int)lierr); 583 } 584 ierr = PetscFPTrapPop();CHKERRQ(ierr); 585 #else /* on missing GESVD */ 586 /* SVD */ 587 PetscInt max_n,min_n; 588 max_n = max_size_of_constraint; 589 min_n = nnsp_addone+nnsp_size; 590 if (max_size_of_constraint < ( nnsp_addone+nnsp_size ) ) { 591 min_n = max_size_of_constraint; 592 max_n = nnsp_addone+nnsp_size; 593 } 594 ierr = PetscMalloc(min_n*sizeof(PetscReal),&singular_vals);CHKERRQ(ierr); 595 #if defined(PETSC_USE_COMPLEX) 596 ierr = PetscMalloc(5*min_n*sizeof(PetscReal),&rwork);CHKERRQ(ierr); 597 #endif 598 /* now we evaluate the optimal workspace using query with lwork=-1 */ 599 lwork=-1; 600 ierr = PetscBLASIntCast(max_n,&Bs);CHKERRQ(ierr); 601 ierr = PetscBLASIntCast(min_n,&Bt);CHKERRQ(ierr); 602 dummy_int = Bs; 603 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 604 #if !defined(PETSC_USE_COMPLEX) 605 PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&Bs,&Bt,&temp_quadrature_constraint[0],&Bs,singular_vals,&dummy_scalar,&dummy_int,&dummy_scalar,&dummy_int,&temp_work,&lwork,&lierr)); 606 #else 607 PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&Bs,&Bt,&temp_quadrature_constraint[0],&Bs,singular_vals,&dummy_scalar,&dummy_int,&dummy_scalar,&dummy_int,&temp_work,&lwork,rwork,&lierr)); 608 #endif 609 if ( lierr ) { 610 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in query to SVD Lapack routine %d",(int)lierr); 611 } 612 ierr = PetscFPTrapPop();CHKERRQ(ierr); 613 #endif 614 /* Allocate optimal workspace */ 615 ierr = PetscBLASIntCast((PetscInt)PetscRealPart(temp_work),&lwork);CHKERRQ(ierr); 616 total_counts = (PetscInt)lwork; 617 ierr = PetscMalloc(total_counts*sizeof(PetscScalar),&work);CHKERRQ(ierr); 618 } 619 /* get local part of global near null space vectors */ 620 ierr = PetscMalloc(nnsp_size*sizeof(Vec),&localnearnullsp);CHKERRQ(ierr); 621 for (k=0;k<nnsp_size;k++) { 622 ierr = VecDuplicate(pcis->vec1_N,&localnearnullsp[k]);CHKERRQ(ierr); 623 ierr = VecScatterBegin(matis->ctx,nearnullvecs[k],localnearnullsp[k],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 624 ierr = VecScatterEnd(matis->ctx,nearnullvecs[k],localnearnullsp[k],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 625 } 626 /* Now we can loop on constraining sets */ 627 total_counts = 0; 628 temp_indices[0] = 0; 629 /* vertices */ 630 if (ISForVertices) { 631 ierr = ISGetIndices(ISForVertices,(const PetscInt**)&is_indices);CHKERRQ(ierr); 632 if (nnsp_has_cnst) { /* consider all vertices */ 633 for (i=0;i<n_vertices;i++) { 634 temp_indices_to_constraint[temp_indices[total_counts]]=is_indices[i]; 635 temp_indices_to_constraint_B[temp_indices[total_counts]]=local_to_B[is_indices[i]]; 636 temp_quadrature_constraint[temp_indices[total_counts]]=1.0; 637 temp_indices[total_counts+1]=temp_indices[total_counts]+1; 638 change_basis[total_counts]=PETSC_FALSE; 639 total_counts++; 640 } 641 } else { /* consider vertices for which exist at least a localnearnullsp which is not null there */ 642 for (i=0;i<n_vertices;i++) { 643 used_vertex=PETSC_FALSE; 644 k=0; 645 while (!used_vertex && k<nnsp_size) { 646 ierr = VecGetArrayRead(localnearnullsp[k],(const PetscScalar**)&array_vector);CHKERRQ(ierr); 647 if (PetscAbsScalar(array_vector[is_indices[i]])>0.0) { 648 temp_indices_to_constraint[temp_indices[total_counts]]=is_indices[i]; 649 temp_indices_to_constraint_B[temp_indices[total_counts]]=local_to_B[is_indices[i]]; 650 temp_quadrature_constraint[temp_indices[total_counts]]=1.0; 651 temp_indices[total_counts+1]=temp_indices[total_counts]+1; 652 change_basis[total_counts]=PETSC_FALSE; 653 total_counts++; 654 used_vertex=PETSC_TRUE; 655 } 656 ierr = VecRestoreArrayRead(localnearnullsp[k],(const PetscScalar**)&array_vector);CHKERRQ(ierr); 657 k++; 658 } 659 } 660 } 661 ierr = ISRestoreIndices(ISForVertices,(const PetscInt**)&is_indices);CHKERRQ(ierr); 662 n_vertices = total_counts; 663 } 664 /* edges and faces */ 665 for (i=0;i<n_ISForEdges+n_ISForFaces;i++) { 666 if (i<n_ISForEdges) { 667 used_IS = &ISForEdges[i]; 668 boolforchange = pcbddc->use_change_of_basis; 669 } else { 670 used_IS = &ISForFaces[i-n_ISForEdges]; 671 boolforchange = pcbddc->use_change_on_faces; 672 } 673 temp_constraints = 0; /* zero the number of constraints I have on this conn comp */ 674 temp_start_ptr = total_counts; /* need to know the starting index of constraints stored */ 675 ierr = ISGetSize(*used_IS,&size_of_constraint);CHKERRQ(ierr); 676 ierr = ISGetIndices(*used_IS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 677 /* HACK: change of basis should not performed on local periodic nodes */ 678 if (pcbddc->mat_graph->mirrors && pcbddc->mat_graph->mirrors[is_indices[0]]) { 679 boolforchange = PETSC_FALSE; 680 } 681 if (nnsp_has_cnst) { 682 PetscScalar quad_value; 683 temp_constraints++; 684 quad_value = (PetscScalar)(1.0/PetscSqrtReal((PetscReal)size_of_constraint)); 685 for (j=0;j<size_of_constraint;j++) { 686 temp_indices_to_constraint[temp_indices[total_counts]+j]=is_indices[j]; 687 temp_indices_to_constraint_B[temp_indices[total_counts]+j]=local_to_B[is_indices[j]]; 688 temp_quadrature_constraint[temp_indices[total_counts]+j]=quad_value; 689 } 690 temp_indices[total_counts+1]=temp_indices[total_counts]+size_of_constraint; /* store new starting point */ 691 change_basis[total_counts]=boolforchange; 692 total_counts++; 693 } 694 for (k=0;k<nnsp_size;k++) { 695 ierr = VecGetArrayRead(localnearnullsp[k],(const PetscScalar**)&array_vector);CHKERRQ(ierr); 696 for (j=0;j<size_of_constraint;j++) { 697 temp_indices_to_constraint[temp_indices[total_counts]+j]=is_indices[j]; 698 temp_indices_to_constraint_B[temp_indices[total_counts]+j]=local_to_B[is_indices[j]]; 699 temp_quadrature_constraint[temp_indices[total_counts]+j]=array_vector[is_indices[j]]; 700 } 701 ierr = VecRestoreArrayRead(localnearnullsp[k],(const PetscScalar**)&array_vector);CHKERRQ(ierr); 702 real_value = 1.0; 703 if (use_nnsp_true) { /* check if array is null on the connected component in case use_nnsp_true has been requested */ 704 ierr = PetscBLASIntCast(size_of_constraint,&Bs);CHKERRQ(ierr); 705 PetscStackCallBLAS("BLASasum",real_value = BLASasum_(&Bs,&temp_quadrature_constraint[temp_indices[total_counts]],&Bone)); 706 } 707 if (real_value > 0.0) { /* keep indices and values */ 708 temp_constraints++; 709 temp_indices[total_counts+1]=temp_indices[total_counts]+size_of_constraint; /* store new starting point */ 710 change_basis[total_counts]=boolforchange; 711 total_counts++; 712 } 713 } 714 ierr = ISRestoreIndices(*used_IS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 715 /* perform SVD on the constraint if use_nnsp_true has not be requested by the user */ 716 if (!use_nnsp_true) { 717 ierr = PetscBLASIntCast(size_of_constraint,&Bs);CHKERRQ(ierr); 718 ierr = PetscBLASIntCast(temp_constraints,&Bt);CHKERRQ(ierr); 719 720 #if defined(PETSC_MISSING_LAPACK_GESVD) 721 ierr = PetscMemzero(correlation_mat,Bt*Bt*sizeof(PetscScalar));CHKERRQ(ierr); 722 /* Store upper triangular part of correlation matrix */ 723 for (j=0;j<temp_constraints;j++) { 724 for (k=0;k<j+1;k++) { 725 PetscStackCallBLAS("BLASdot",correlation_mat[j*temp_constraints+k]=BLASdot_(&Bs,&temp_quadrature_constraint[temp_indices[temp_start_ptr+j]],&Bone,&temp_quadrature_constraint[temp_indices[temp_start_ptr+k]],&Bone)); 726 727 } 728 } 729 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 730 #if !defined(PETSC_USE_COMPLEX) 731 /* LAPACKsyev_("V","U",&Bt,correlation_mat,&Bt,singular_vals,work,&lwork,&lierr); */ 732 PetscStackCallBLAS("LAPACKsyevx",LAPACKsyevx_("V","A","U",&Bt,correlation_mat,&Bt,&dummy_real,&dummy_real,&dummy_int,&dummy_int,&abs_tol,&eigs_found,singular_vals,singular_vectors,&Bt,work,&lwork,iwork,ifail,&lierr)); 733 #else 734 /* LAPACK call is missing here! TODO */ 735 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not yet implemented for complexes when PETSC_MISSING_GESVD = 1"); 736 #endif 737 if (lierr) { 738 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SYEVX Lapack routine %d",(int)lierr); 739 } 740 ierr = PetscFPTrapPop();CHKERRQ(ierr); 741 /* retain eigenvalues greater than tol: note that lapack SYEV gives eigs in ascending order */ 742 j=0; 743 while (j < Bt && singular_vals[j] < tol) j++; 744 total_counts=total_counts-j; 745 if (j<temp_constraints) { 746 for (k=j;k<Bt;k++) { 747 singular_vals[k]=1.0/PetscSqrtReal(singular_vals[k]); 748 } 749 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 750 PetscStackCallBLAS("BLASgemm_",BLASgemm_("N","N",&Bs,&Bt,&Bt,&one,&temp_quadrature_constraint[temp_indices[temp_start_ptr]],&Bs,correlation_mat,&Bt,&zero,temp_basis,&Bs)); 751 ierr = PetscFPTrapPop();CHKERRQ(ierr); 752 /* copy POD basis into used quadrature memory */ 753 for (k=0;k<Bt-j;k++) { 754 for (ii=0;ii<size_of_constraint;ii++) { 755 temp_quadrature_constraint[temp_indices[temp_start_ptr+k]+ii]=singular_vals[Bt-1-k]*temp_basis[(Bt-1-k)*size_of_constraint+ii]; 756 } 757 } 758 } 759 760 #else /* on missing GESVD */ 761 PetscInt min_n = temp_constraints; 762 if (min_n > size_of_constraint) min_n = size_of_constraint; 763 dummy_int = Bs; 764 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 765 #if !defined(PETSC_USE_COMPLEX) 766 PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&Bs,&Bt,&temp_quadrature_constraint[temp_indices[temp_start_ptr]],&Bs,singular_vals,&dummy_scalar,&dummy_int,&dummy_scalar,&dummy_int,work,&lwork,&lierr)); 767 #else 768 PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&Bs,&Bt,&temp_quadrature_constraint[temp_indices[temp_start_ptr]],&Bs,singular_vals,&dummy_scalar,&dummy_int,&dummy_scalar,&dummy_int,work,&lwork,rwork,&lierr)); 769 #endif 770 if (lierr) { 771 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SVD Lapack routine %d",(int)lierr); 772 } 773 ierr = PetscFPTrapPop();CHKERRQ(ierr); 774 /* retain eigenvalues greater than tol: note that lapack SVD gives eigs in descending order */ 775 j = 0; 776 while (j < min_n && singular_vals[min_n-j-1] < tol) j++; 777 total_counts = total_counts-(PetscInt)Bt+(min_n-j); 778 #endif 779 } 780 } 781 /* free index sets of faces, edges and vertices */ 782 for (i=0;i<n_ISForFaces;i++) { 783 ierr = ISDestroy(&ISForFaces[i]);CHKERRQ(ierr); 784 } 785 ierr = PetscFree(ISForFaces);CHKERRQ(ierr); 786 for (i=0;i<n_ISForEdges;i++) { 787 ierr = ISDestroy(&ISForEdges[i]);CHKERRQ(ierr); 788 } 789 ierr = PetscFree(ISForEdges);CHKERRQ(ierr); 790 ierr = ISDestroy(&ISForVertices);CHKERRQ(ierr); 791 792 /* set quantities in pcbddc data structure */ 793 /* n_vertices defines the number of point primal dofs */ 794 /* n_constraints defines the number of averages (they can be point primal dofs if change of basis is requested) */ 795 local_primal_size = total_counts; 796 pcbddc->n_vertices = n_vertices; 797 pcbddc->n_constraints = total_counts-n_vertices; 798 pcbddc->local_primal_size = local_primal_size; 799 800 /* Create constraint matrix */ 801 /* The constraint matrix is used to compute the l2g map of primal dofs */ 802 /* so we need to set it up properly either with or without change of basis */ 803 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->ConstraintMatrix);CHKERRQ(ierr); 804 ierr = MatSetType(pcbddc->ConstraintMatrix,impMatType);CHKERRQ(ierr); 805 ierr = MatSetSizes(pcbddc->ConstraintMatrix,local_primal_size,pcis->n,local_primal_size,pcis->n);CHKERRQ(ierr); 806 /* compute a local numbering of constraints : vertices first then constraints */ 807 ierr = VecSet(pcis->vec1_N,0.0);CHKERRQ(ierr); 808 ierr = VecGetArray(pcis->vec1_N,&array_vector);CHKERRQ(ierr); 809 ierr = PetscMalloc(local_primal_size*sizeof(PetscInt),&aux_primal_numbering);CHKERRQ(ierr); 810 ierr = PetscMalloc(local_primal_size*sizeof(PetscInt),&aux_primal_permutation);CHKERRQ(ierr); 811 total_counts=0; 812 /* find vertices: subdomain corners plus dofs with basis changed */ 813 for (i=0;i<local_primal_size;i++) { 814 size_of_constraint=temp_indices[i+1]-temp_indices[i]; 815 if (change_basis[i] || size_of_constraint == 1) { 816 k=0; 817 while(k < size_of_constraint && array_vector[temp_indices_to_constraint[temp_indices[i]+size_of_constraint-k-1]] != 0.0) { 818 k=k+1; 819 } 820 j=temp_indices_to_constraint[temp_indices[i]+size_of_constraint-k-1]; 821 array_vector[j] = 1.0; 822 aux_primal_numbering[total_counts]=j; 823 aux_primal_permutation[total_counts]=total_counts; 824 total_counts++; 825 } 826 } 827 ierr = VecRestoreArray(pcis->vec1_N,&array_vector);CHKERRQ(ierr); 828 /* permute indices in order to have a sorted set of vertices */ 829 ierr = PetscSortIntWithPermutation(total_counts,aux_primal_numbering,aux_primal_permutation); 830 /* nonzero structure */ 831 ierr = PetscMalloc(local_primal_size*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 832 for (i=0;i<total_counts;i++) { 833 nnz[i]=1; 834 } 835 j=total_counts; 836 for (i=n_vertices;i<local_primal_size;i++) { 837 if (!change_basis[i]) { 838 nnz[j]=temp_indices[i+1]-temp_indices[i]; 839 j++; 840 } 841 } 842 ierr = MatSeqAIJSetPreallocation(pcbddc->ConstraintMatrix,0,nnz);CHKERRQ(ierr); 843 ierr = PetscFree(nnz);CHKERRQ(ierr); 844 /* set values in constraint matrix */ 845 for (i=0;i<total_counts;i++) { 846 j = aux_primal_permutation[i]; 847 k = aux_primal_numbering[j]; 848 ierr = MatSetValue(pcbddc->ConstraintMatrix,i,k,1.0,INSERT_VALUES);CHKERRQ(ierr); 849 } 850 for (i=n_vertices;i<local_primal_size;i++) { 851 if (!change_basis[i]) { 852 size_of_constraint=temp_indices[i+1]-temp_indices[i]; 853 ierr = MatSetValues(pcbddc->ConstraintMatrix,1,&total_counts,size_of_constraint,&temp_indices_to_constraint[temp_indices[i]],&temp_quadrature_constraint[temp_indices[i]],INSERT_VALUES);CHKERRQ(ierr); 854 total_counts++; 855 } 856 } 857 ierr = PetscFree(aux_primal_numbering);CHKERRQ(ierr); 858 ierr = PetscFree(aux_primal_permutation);CHKERRQ(ierr); 859 /* assembling */ 860 ierr = MatAssemblyBegin(pcbddc->ConstraintMatrix,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 861 ierr = MatAssemblyEnd(pcbddc->ConstraintMatrix,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 862 863 /* Create matrix for change of basis. We don't need it in case pcbddc->use_change_of_basis is FALSE */ 864 if (pcbddc->use_change_of_basis) { 865 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr); 866 ierr = MatSetType(pcbddc->ChangeOfBasisMatrix,impMatType);CHKERRQ(ierr); 867 ierr = MatSetSizes(pcbddc->ChangeOfBasisMatrix,pcis->n_B,pcis->n_B,pcis->n_B,pcis->n_B);CHKERRQ(ierr); 868 /* work arrays */ 869 /* we need to reuse these arrays, so we free them */ 870 ierr = PetscFree(temp_basis);CHKERRQ(ierr); 871 ierr = PetscFree(work);CHKERRQ(ierr); 872 ierr = PetscMalloc(pcis->n_B*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 873 ierr = PetscMalloc((nnsp_addone+nnsp_size)*(nnsp_addone+nnsp_size)*sizeof(PetscScalar),&temp_basis);CHKERRQ(ierr); 874 ierr = PetscMalloc((nnsp_addone+nnsp_size)*sizeof(PetscScalar),&work);CHKERRQ(ierr); 875 ierr = PetscMalloc((nnsp_addone+nnsp_size)*sizeof(PetscBLASInt),&ipiv);CHKERRQ(ierr); 876 for (i=0;i<pcis->n_B;i++) { 877 nnz[i]=1; 878 } 879 /* Overestimated nonzeros per row */ 880 k=1; 881 for (i=pcbddc->n_vertices;i<local_primal_size;i++) { 882 if (change_basis[i]) { 883 size_of_constraint = temp_indices[i+1]-temp_indices[i]; 884 if (k < size_of_constraint) { 885 k = size_of_constraint; 886 } 887 for (j=0;j<size_of_constraint;j++) { 888 nnz[temp_indices_to_constraint_B[temp_indices[i]+j]] = size_of_constraint; 889 } 890 } 891 } 892 ierr = MatSeqAIJSetPreallocation(pcbddc->ChangeOfBasisMatrix,0,nnz);CHKERRQ(ierr); 893 ierr = PetscFree(nnz);CHKERRQ(ierr); 894 /* Temporary array to store indices */ 895 ierr = PetscMalloc(k*sizeof(PetscInt),&is_indices);CHKERRQ(ierr); 896 /* Set initial identity in the matrix */ 897 for (i=0;i<pcis->n_B;i++) { 898 ierr = MatSetValue(pcbddc->ChangeOfBasisMatrix,i,i,1.0,INSERT_VALUES);CHKERRQ(ierr); 899 } 900 /* Now we loop on the constraints which need a change of basis */ 901 /* Change of basis matrix is evaluated as the FIRST APPROACH in */ 902 /* Klawonn and Widlund, Dual-primal FETI-DP methods for linear elasticity, (6.2.1) */ 903 temp_constraints = 0; 904 if (pcbddc->n_vertices < local_primal_size) { 905 temp_start_ptr = temp_indices_to_constraint_B[temp_indices[pcbddc->n_vertices]]; 906 } 907 for (i=pcbddc->n_vertices;i<local_primal_size;i++) { 908 if (change_basis[i]) { 909 compute_submatrix = PETSC_FALSE; 910 useksp = PETSC_FALSE; 911 if (temp_start_ptr == temp_indices_to_constraint_B[temp_indices[i]]) { 912 temp_constraints++; 913 if (i == local_primal_size -1 || temp_start_ptr != temp_indices_to_constraint_B[temp_indices[i+1]]) { 914 compute_submatrix = PETSC_TRUE; 915 } 916 } 917 if (compute_submatrix) { 918 if (temp_constraints > 1 || pcbddc->use_nnsp_true) { 919 useksp = PETSC_TRUE; 920 } 921 size_of_constraint = temp_indices[i+1]-temp_indices[i]; 922 if (useksp) { /* experimental TODO: reuse KSP and MAT instead of creating them each time */ 923 ierr = MatCreate(PETSC_COMM_SELF,&temp_mat);CHKERRQ(ierr); 924 ierr = MatSetType(temp_mat,impMatType);CHKERRQ(ierr); 925 ierr = MatSetSizes(temp_mat,size_of_constraint,size_of_constraint,size_of_constraint,size_of_constraint);CHKERRQ(ierr); 926 ierr = MatSeqAIJSetPreallocation(temp_mat,size_of_constraint,NULL);CHKERRQ(ierr); 927 } 928 /* First _size_of_constraint-temp_constraints_ columns */ 929 dual_dofs = size_of_constraint-temp_constraints; 930 start_constraint = i+1-temp_constraints; 931 for (s=0;s<dual_dofs;s++) { 932 is_indices[0] = s; 933 for (j=0;j<temp_constraints;j++) { 934 for (k=0;k<temp_constraints;k++) { 935 temp_basis[j*temp_constraints+k]=temp_quadrature_constraint[temp_indices[start_constraint+k]+s+j+1]; 936 } 937 work[j]=-temp_quadrature_constraint[temp_indices[start_constraint+j]+s]; 938 is_indices[j+1]=s+j+1; 939 } 940 Bt = temp_constraints; 941 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 942 PetscStackCallBLAS("LAPACKgesv",LAPACKgesv_(&Bt,&Bone,temp_basis,&Bt,ipiv,work,&Bt,&lierr)); 943 if ( lierr ) { 944 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GESV Lapack routine %d",(int)lierr); 945 } 946 ierr = PetscFPTrapPop();CHKERRQ(ierr); 947 j = temp_indices_to_constraint_B[temp_indices[start_constraint]+s]; 948 ierr = MatSetValues(pcbddc->ChangeOfBasisMatrix,temp_constraints,&temp_indices_to_constraint_B[temp_indices[start_constraint]+s+1],1,&j,work,INSERT_VALUES);CHKERRQ(ierr); 949 if (useksp) { 950 /* temp mat with transposed rows and columns */ 951 ierr = MatSetValues(temp_mat,1,&s,temp_constraints,&is_indices[1],work,INSERT_VALUES);CHKERRQ(ierr); 952 ierr = MatSetValue(temp_mat,is_indices[0],is_indices[0],1.0,INSERT_VALUES);CHKERRQ(ierr); 953 } 954 } 955 if (useksp) { 956 /* last rows of temp_mat */ 957 for (j=0;j<size_of_constraint;j++) { 958 is_indices[j] = j; 959 } 960 for (s=0;s<temp_constraints;s++) { 961 k = s + dual_dofs; 962 ierr = MatSetValues(temp_mat,1,&k,size_of_constraint,is_indices,&temp_quadrature_constraint[temp_indices[start_constraint+s]],INSERT_VALUES);CHKERRQ(ierr); 963 } 964 ierr = MatAssemblyBegin(temp_mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 965 ierr = MatAssemblyEnd(temp_mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 966 ierr = MatGetVecs(temp_mat,&temp_vec,NULL);CHKERRQ(ierr); 967 ierr = KSPCreate(PETSC_COMM_SELF,&temp_ksp);CHKERRQ(ierr); 968 ierr = KSPSetOperators(temp_ksp,temp_mat,temp_mat,SAME_PRECONDITIONER);CHKERRQ(ierr); 969 ierr = KSPSetType(temp_ksp,KSPPREONLY);CHKERRQ(ierr); 970 ierr = KSPGetPC(temp_ksp,&temp_pc);CHKERRQ(ierr); 971 ierr = PCSetType(temp_pc,PCLU);CHKERRQ(ierr); 972 ierr = KSPSetUp(temp_ksp);CHKERRQ(ierr); 973 for (s=0;s<temp_constraints;s++) { 974 ierr = VecSet(temp_vec,0.0);CHKERRQ(ierr); 975 ierr = VecSetValue(temp_vec,s+dual_dofs,1.0,INSERT_VALUES);CHKERRQ(ierr); 976 ierr = VecAssemblyBegin(temp_vec);CHKERRQ(ierr); 977 ierr = VecAssemblyEnd(temp_vec);CHKERRQ(ierr); 978 ierr = KSPSolve(temp_ksp,temp_vec,temp_vec);CHKERRQ(ierr); 979 ierr = VecGetArray(temp_vec,&array_vector);CHKERRQ(ierr); 980 j = temp_indices_to_constraint_B[temp_indices[start_constraint+s]+size_of_constraint-s-1]; 981 /* last columns of change of basis matrix associated to new primal dofs */ 982 ierr = MatSetValues(pcbddc->ChangeOfBasisMatrix,size_of_constraint,&temp_indices_to_constraint_B[temp_indices[start_constraint+s]],1,&j,array_vector,INSERT_VALUES);CHKERRQ(ierr); 983 ierr = VecRestoreArray(temp_vec,&array_vector);CHKERRQ(ierr); 984 } 985 ierr = MatDestroy(&temp_mat);CHKERRQ(ierr); 986 ierr = KSPDestroy(&temp_ksp);CHKERRQ(ierr); 987 ierr = VecDestroy(&temp_vec);CHKERRQ(ierr); 988 } else { 989 /* last columns of change of basis matrix associated to new primal dofs */ 990 for (s=0;s<temp_constraints;s++) { 991 j = temp_indices_to_constraint_B[temp_indices[start_constraint+s]+size_of_constraint-s-1]; 992 ierr = MatSetValues(pcbddc->ChangeOfBasisMatrix,size_of_constraint,&temp_indices_to_constraint_B[temp_indices[start_constraint+s]],1,&j,&temp_quadrature_constraint[temp_indices[start_constraint+s]],INSERT_VALUES);CHKERRQ(ierr); 993 } 994 } 995 /* prepare for the next cycle */ 996 temp_constraints = 0; 997 if (i != local_primal_size -1 ) { 998 temp_start_ptr = temp_indices_to_constraint_B[temp_indices[i+1]]; 999 } 1000 } 1001 } 1002 } 1003 /* assembling */ 1004 ierr = MatAssemblyBegin(pcbddc->ChangeOfBasisMatrix,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1005 ierr = MatAssemblyEnd(pcbddc->ChangeOfBasisMatrix,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1006 ierr = PetscFree(ipiv);CHKERRQ(ierr); 1007 ierr = PetscFree(is_indices);CHKERRQ(ierr); 1008 } 1009 /* free workspace no longer needed */ 1010 ierr = PetscFree(rwork);CHKERRQ(ierr); 1011 ierr = PetscFree(work);CHKERRQ(ierr); 1012 ierr = PetscFree(temp_basis);CHKERRQ(ierr); 1013 ierr = PetscFree(singular_vals);CHKERRQ(ierr); 1014 ierr = PetscFree(correlation_mat);CHKERRQ(ierr); 1015 ierr = PetscFree(temp_indices);CHKERRQ(ierr); 1016 ierr = PetscFree(change_basis);CHKERRQ(ierr); 1017 ierr = PetscFree(temp_indices_to_constraint);CHKERRQ(ierr); 1018 ierr = PetscFree(temp_indices_to_constraint_B);CHKERRQ(ierr); 1019 ierr = PetscFree(local_to_B);CHKERRQ(ierr); 1020 ierr = PetscFree(temp_quadrature_constraint);CHKERRQ(ierr); 1021 #if defined(PETSC_MISSING_LAPACK_GESVD) 1022 ierr = PetscFree(iwork);CHKERRQ(ierr); 1023 ierr = PetscFree(ifail);CHKERRQ(ierr); 1024 ierr = PetscFree(singular_vectors);CHKERRQ(ierr); 1025 #endif 1026 for (k=0;k<nnsp_size;k++) { 1027 ierr = VecDestroy(&localnearnullsp[k]);CHKERRQ(ierr); 1028 } 1029 ierr = PetscFree(localnearnullsp);CHKERRQ(ierr); 1030 PetscFunctionReturn(0); 1031 } 1032 1033 #undef __FUNCT__ 1034 #define __FUNCT__ "PCBDDCAnalyzeInterface" 1035 PetscErrorCode PCBDDCAnalyzeInterface(PC pc) 1036 { 1037 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1038 PC_IS *pcis = (PC_IS*)pc->data; 1039 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 1040 PetscInt bs,ierr,i,vertex_size; 1041 PetscViewer viewer=pcbddc->dbg_viewer; 1042 1043 PetscFunctionBegin; 1044 /* Init local Graph struct */ 1045 ierr = PCBDDCGraphInit(pcbddc->mat_graph,matis->mapping);CHKERRQ(ierr); 1046 1047 /* Check validity of the csr graph passed in by the user */ 1048 if (pcbddc->mat_graph->nvtxs_csr != pcbddc->mat_graph->nvtxs) { 1049 ierr = PCBDDCGraphResetCSR(pcbddc->mat_graph);CHKERRQ(ierr); 1050 } 1051 /* Set default CSR adjacency of local dofs if not provided by the user with PCBDDCSetLocalAdjacencyGraph */ 1052 if (!pcbddc->mat_graph->xadj || !pcbddc->mat_graph->adjncy) { 1053 Mat mat_adj; 1054 const PetscInt *xadj,*adjncy; 1055 PetscBool flg_row=PETSC_TRUE; 1056 1057 ierr = MatConvert(matis->A,MATMPIADJ,MAT_INITIAL_MATRIX,&mat_adj);CHKERRQ(ierr); 1058 ierr = MatGetRowIJ(mat_adj,0,PETSC_TRUE,PETSC_FALSE,&i,&xadj,&adjncy,&flg_row);CHKERRQ(ierr); 1059 if (!flg_row) { 1060 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in MatGetRowIJ called in %s\n",__FUNCT__); 1061 } 1062 ierr = PCBDDCSetLocalAdjacencyGraph(pc,i,xadj,adjncy,PETSC_COPY_VALUES);CHKERRQ(ierr); 1063 ierr = MatRestoreRowIJ(mat_adj,0,PETSC_TRUE,PETSC_FALSE,&i,&xadj,&adjncy,&flg_row);CHKERRQ(ierr); 1064 if (!flg_row) { 1065 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in MatRestoreRowIJ called in %s\n",__FUNCT__); 1066 } 1067 ierr = MatDestroy(&mat_adj);CHKERRQ(ierr); 1068 } 1069 1070 /* Set default dofs' splitting if no information has been provided by the user with PCBDDCSetDofsSplitting */ 1071 vertex_size = 1; 1072 if (!pcbddc->n_ISForDofs) { 1073 IS *custom_ISForDofs; 1074 1075 ierr = MatGetBlockSize(matis->A,&bs);CHKERRQ(ierr); 1076 ierr = PetscMalloc(bs*sizeof(IS),&custom_ISForDofs);CHKERRQ(ierr); 1077 for (i=0;i<bs;i++) { 1078 ierr = ISCreateStride(PETSC_COMM_SELF,pcis->n/bs,i,bs,&custom_ISForDofs[i]);CHKERRQ(ierr); 1079 } 1080 ierr = PCBDDCSetDofsSplitting(pc,bs,custom_ISForDofs);CHKERRQ(ierr); 1081 /* remove my references to IS objects */ 1082 for (i=0;i<bs;i++) { 1083 ierr = ISDestroy(&custom_ISForDofs[i]);CHKERRQ(ierr); 1084 } 1085 ierr = PetscFree(custom_ISForDofs);CHKERRQ(ierr); 1086 } else { /* mat block size as vertex size (used for elasticity) */ 1087 ierr = MatGetBlockSize(matis->A,&vertex_size);CHKERRQ(ierr); 1088 } 1089 1090 /* Setup of Graph */ 1091 ierr = PCBDDCGraphSetUp(pcbddc->mat_graph,vertex_size,pcbddc->NeumannBoundaries,pcbddc->DirichletBoundaries,pcbddc->n_ISForDofs,pcbddc->ISForDofs,pcbddc->user_primal_vertices); 1092 1093 /* Graph's connected components analysis */ 1094 ierr = PCBDDCGraphComputeConnectedComponents(pcbddc->mat_graph);CHKERRQ(ierr); 1095 1096 /* print some info to stdout */ 1097 if (pcbddc->dbg_flag) { 1098 ierr = PCBDDCGraphASCIIView(pcbddc->mat_graph,pcbddc->dbg_flag,viewer); 1099 } 1100 PetscFunctionReturn(0); 1101 } 1102 1103 #undef __FUNCT__ 1104 #define __FUNCT__ "PCBDDCGetPrimalVerticesLocalIdx" 1105 PetscErrorCode PCBDDCGetPrimalVerticesLocalIdx(PC pc, PetscInt *n_vertices, PetscInt *vertices_idx[]) 1106 { 1107 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1108 PetscInt *vertices,*row_cmat_indices,n,i,size_of_constraint,local_primal_size; 1109 PetscErrorCode ierr; 1110 1111 PetscFunctionBegin; 1112 n = 0; 1113 vertices = 0; 1114 if (pcbddc->ConstraintMatrix) { 1115 ierr = MatGetSize(pcbddc->ConstraintMatrix,&local_primal_size,&i);CHKERRQ(ierr); 1116 for (i=0;i<local_primal_size;i++) { 1117 ierr = MatGetRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,NULL,NULL);CHKERRQ(ierr); 1118 if (size_of_constraint == 1) n++; 1119 ierr = MatRestoreRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,NULL,NULL);CHKERRQ(ierr); 1120 } 1121 ierr = PetscMalloc(n*sizeof(PetscInt),&vertices);CHKERRQ(ierr); 1122 n = 0; 1123 for (i=0;i<local_primal_size;i++) { 1124 ierr = MatGetRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,(const PetscInt**)&row_cmat_indices,NULL);CHKERRQ(ierr); 1125 if (size_of_constraint == 1) { 1126 vertices[n++]=row_cmat_indices[0]; 1127 } 1128 ierr = MatRestoreRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,(const PetscInt**)&row_cmat_indices,NULL);CHKERRQ(ierr); 1129 } 1130 } 1131 *n_vertices = n; 1132 *vertices_idx = vertices; 1133 PetscFunctionReturn(0); 1134 } 1135 1136 #undef __FUNCT__ 1137 #define __FUNCT__ "PCBDDCGetPrimalConstraintsLocalIdx" 1138 PetscErrorCode PCBDDCGetPrimalConstraintsLocalIdx(PC pc, PetscInt *n_constraints, PetscInt *constraints_idx[]) 1139 { 1140 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1141 PetscInt *constraints_index,*row_cmat_indices,*row_cmat_global_indices; 1142 PetscInt n,i,j,size_of_constraint,local_primal_size,local_size,max_size_of_constraint,min_index,min_loc; 1143 PetscBool *touched; 1144 PetscErrorCode ierr; 1145 1146 PetscFunctionBegin; 1147 n = 0; 1148 constraints_index = 0; 1149 if (pcbddc->ConstraintMatrix) { 1150 ierr = MatGetSize(pcbddc->ConstraintMatrix,&local_primal_size,&local_size);CHKERRQ(ierr); 1151 max_size_of_constraint = 0; 1152 for (i=0;i<local_primal_size;i++) { 1153 ierr = MatGetRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,NULL,NULL);CHKERRQ(ierr); 1154 if (size_of_constraint > 1) { 1155 n++; 1156 } 1157 max_size_of_constraint = PetscMax(size_of_constraint,max_size_of_constraint); 1158 ierr = MatRestoreRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,NULL,NULL);CHKERRQ(ierr); 1159 } 1160 ierr = PetscMalloc(n*sizeof(PetscInt),&constraints_index);CHKERRQ(ierr); 1161 ierr = PetscMalloc(max_size_of_constraint*sizeof(PetscInt),&row_cmat_global_indices);CHKERRQ(ierr); 1162 ierr = PetscMalloc(local_size*sizeof(PetscBool),&touched);CHKERRQ(ierr); 1163 ierr = PetscMemzero(touched,local_size*sizeof(PetscBool));CHKERRQ(ierr); 1164 n = 0; 1165 for (i=0;i<local_primal_size;i++) { 1166 ierr = MatGetRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,(const PetscInt**)&row_cmat_indices,NULL);CHKERRQ(ierr); 1167 if (size_of_constraint > 1) { 1168 ierr = ISLocalToGlobalMappingApply(pcbddc->mat_graph->l2gmap,size_of_constraint,row_cmat_indices,row_cmat_global_indices);CHKERRQ(ierr); 1169 min_index = row_cmat_global_indices[0]; 1170 min_loc = 0; 1171 for (j=1;j<size_of_constraint;j++) { 1172 /* there can be more than one constraint on a single connected component */ 1173 if (min_index > row_cmat_global_indices[j] && !touched[row_cmat_indices[j]]) { 1174 min_index = row_cmat_global_indices[j]; 1175 min_loc = j; 1176 } 1177 } 1178 touched[row_cmat_indices[min_loc]] = PETSC_TRUE; 1179 constraints_index[n++] = row_cmat_indices[min_loc]; 1180 } 1181 ierr = MatRestoreRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,(const PetscInt**)&row_cmat_indices,NULL);CHKERRQ(ierr); 1182 } 1183 } 1184 ierr = PetscFree(touched);CHKERRQ(ierr); 1185 ierr = PetscFree(row_cmat_global_indices);CHKERRQ(ierr); 1186 *n_constraints = n; 1187 *constraints_idx = constraints_index; 1188 PetscFunctionReturn(0); 1189 } 1190 1191 /* the next two functions has been adapted from pcis.c */ 1192 #undef __FUNCT__ 1193 #define __FUNCT__ "PCBDDCApplySchur" 1194 PetscErrorCode PCBDDCApplySchur(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D) 1195 { 1196 PetscErrorCode ierr; 1197 PC_IS *pcis = (PC_IS*)(pc->data); 1198 1199 PetscFunctionBegin; 1200 if (!vec2_B) { vec2_B = v; } 1201 ierr = MatMult(pcis->A_BB,v,vec1_B);CHKERRQ(ierr); 1202 ierr = MatMult(pcis->A_IB,v,vec1_D);CHKERRQ(ierr); 1203 ierr = KSPSolve(pcis->ksp_D,vec1_D,vec2_D);CHKERRQ(ierr); 1204 ierr = MatMult(pcis->A_BI,vec2_D,vec2_B);CHKERRQ(ierr); 1205 ierr = VecAXPY(vec1_B,-1.0,vec2_B);CHKERRQ(ierr); 1206 PetscFunctionReturn(0); 1207 } 1208 1209 #undef __FUNCT__ 1210 #define __FUNCT__ "PCBDDCApplySchurTranspose" 1211 PetscErrorCode PCBDDCApplySchurTranspose(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D) 1212 { 1213 PetscErrorCode ierr; 1214 PC_IS *pcis = (PC_IS*)(pc->data); 1215 1216 PetscFunctionBegin; 1217 if (!vec2_B) { vec2_B = v; } 1218 ierr = MatMultTranspose(pcis->A_BB,v,vec1_B);CHKERRQ(ierr); 1219 ierr = MatMultTranspose(pcis->A_BI,v,vec1_D);CHKERRQ(ierr); 1220 ierr = KSPSolveTranspose(pcis->ksp_D,vec1_D,vec2_D);CHKERRQ(ierr); 1221 ierr = MatMultTranspose(pcis->A_IB,vec2_D,vec2_B);CHKERRQ(ierr); 1222 ierr = VecAXPY(vec1_B,-1.0,vec2_B);CHKERRQ(ierr); 1223 PetscFunctionReturn(0); 1224 } 1225 1226 #undef __FUNCT__ 1227 #define __FUNCT__ "PCBDDCSubsetNumbering" 1228 PetscErrorCode PCBDDCSubsetNumbering(MPI_Comm comm,ISLocalToGlobalMapping l2gmap, PetscInt n_local_dofs, PetscInt local_dofs[], PetscInt local_dofs_mult[], PetscInt* n_global_subset, PetscInt* global_numbering_subset[]) 1229 { 1230 Vec local_vec,global_vec; 1231 IS seqis,paris; 1232 VecScatter scatter_ctx; 1233 PetscScalar *array; 1234 PetscInt *temp_global_dofs; 1235 PetscScalar globalsum; 1236 PetscInt i,j,s; 1237 PetscInt nlocals,first_index,old_index,max_local; 1238 PetscMPIInt rank_prec_comm,size_prec_comm,max_global; 1239 PetscMPIInt *dof_sizes,*dof_displs; 1240 PetscBool first_found; 1241 PetscErrorCode ierr; 1242 1243 PetscFunctionBegin; 1244 /* mpi buffers */ 1245 MPI_Comm_size(comm,&size_prec_comm); 1246 MPI_Comm_rank(comm,&rank_prec_comm); 1247 j = ( !rank_prec_comm ? size_prec_comm : 0); 1248 ierr = PetscMalloc(j*sizeof(*dof_sizes),&dof_sizes);CHKERRQ(ierr); 1249 ierr = PetscMalloc(j*sizeof(*dof_displs),&dof_displs);CHKERRQ(ierr); 1250 /* get maximum size of subset */ 1251 ierr = PetscMalloc(n_local_dofs*sizeof(PetscInt),&temp_global_dofs);CHKERRQ(ierr); 1252 ierr = ISLocalToGlobalMappingApply(l2gmap,n_local_dofs,local_dofs,temp_global_dofs);CHKERRQ(ierr); 1253 max_local = 0; 1254 if (n_local_dofs) { 1255 max_local = temp_global_dofs[0]; 1256 for (i=1;i<n_local_dofs;i++) { 1257 if (max_local < temp_global_dofs[i] ) { 1258 max_local = temp_global_dofs[i]; 1259 } 1260 } 1261 } 1262 ierr = MPI_Allreduce(&max_local,&max_global,1,MPIU_INT,MPI_MAX,comm); 1263 max_global++; 1264 max_local = 0; 1265 if (n_local_dofs) { 1266 max_local = local_dofs[0]; 1267 for (i=1;i<n_local_dofs;i++) { 1268 if (max_local < local_dofs[i] ) { 1269 max_local = local_dofs[i]; 1270 } 1271 } 1272 } 1273 max_local++; 1274 /* allocate workspace */ 1275 ierr = VecCreate(PETSC_COMM_SELF,&local_vec);CHKERRQ(ierr); 1276 ierr = VecSetSizes(local_vec,PETSC_DECIDE,max_local);CHKERRQ(ierr); 1277 ierr = VecSetType(local_vec,VECSEQ);CHKERRQ(ierr); 1278 ierr = VecCreate(comm,&global_vec);CHKERRQ(ierr); 1279 ierr = VecSetSizes(global_vec,PETSC_DECIDE,max_global);CHKERRQ(ierr); 1280 ierr = VecSetType(global_vec,VECMPI);CHKERRQ(ierr); 1281 /* create scatter */ 1282 ierr = ISCreateGeneral(PETSC_COMM_SELF,n_local_dofs,local_dofs,PETSC_COPY_VALUES,&seqis);CHKERRQ(ierr); 1283 ierr = ISCreateGeneral(comm,n_local_dofs,temp_global_dofs,PETSC_COPY_VALUES,&paris);CHKERRQ(ierr); 1284 ierr = VecScatterCreate(local_vec,seqis,global_vec,paris,&scatter_ctx);CHKERRQ(ierr); 1285 ierr = ISDestroy(&seqis);CHKERRQ(ierr); 1286 ierr = ISDestroy(&paris);CHKERRQ(ierr); 1287 /* init array */ 1288 ierr = VecSet(global_vec,0.0);CHKERRQ(ierr); 1289 ierr = VecSet(local_vec,0.0);CHKERRQ(ierr); 1290 ierr = VecGetArray(local_vec,&array);CHKERRQ(ierr); 1291 if (local_dofs_mult) { 1292 for (i=0;i<n_local_dofs;i++) { 1293 array[local_dofs[i]]=(PetscScalar)local_dofs_mult[i]; 1294 } 1295 } else { 1296 for (i=0;i<n_local_dofs;i++) { 1297 array[local_dofs[i]]=1.0; 1298 } 1299 } 1300 ierr = VecRestoreArray(local_vec,&array);CHKERRQ(ierr); 1301 /* scatter into global vec and get total number of global dofs */ 1302 ierr = VecScatterBegin(scatter_ctx,local_vec,global_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1303 ierr = VecScatterEnd(scatter_ctx,local_vec,global_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1304 ierr = VecSum(global_vec,&globalsum);CHKERRQ(ierr); 1305 *n_global_subset = (PetscInt)PetscRealPart(globalsum); 1306 /* Fill global_vec with cumulative function for global numbering */ 1307 ierr = VecGetArray(global_vec,&array);CHKERRQ(ierr); 1308 ierr = VecGetLocalSize(global_vec,&s);CHKERRQ(ierr); 1309 nlocals = 0; 1310 first_index = -1; 1311 first_found = PETSC_FALSE; 1312 for (i=0;i<s;i++) { 1313 if (!first_found && PetscRealPart(array[i]) > 0.0) { 1314 first_found = PETSC_TRUE; 1315 first_index = i; 1316 } 1317 nlocals += (PetscInt)PetscRealPart(array[i]); 1318 } 1319 ierr = MPI_Gather(&nlocals,1,MPIU_INT,dof_sizes,1,MPIU_INT,0,comm);CHKERRQ(ierr); 1320 if (!rank_prec_comm) { 1321 dof_displs[0]=0; 1322 for (i=1;i<size_prec_comm;i++) { 1323 dof_displs[i] = dof_displs[i-1]+dof_sizes[i-1]; 1324 } 1325 } 1326 ierr = MPI_Scatter(dof_displs,1,MPIU_INT,&nlocals,1,MPIU_INT,0,comm);CHKERRQ(ierr); 1327 if (first_found) { 1328 array[first_index] += (PetscScalar)nlocals; 1329 old_index = first_index; 1330 for (i=first_index+1;i<s;i++) { 1331 if (PetscRealPart(array[i]) > 0.0) { 1332 array[i] += array[old_index]; 1333 old_index = i; 1334 } 1335 } 1336 } 1337 ierr = VecRestoreArray(global_vec,&array);CHKERRQ(ierr); 1338 ierr = VecSet(local_vec,0.0);CHKERRQ(ierr); 1339 ierr = VecScatterBegin(scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1340 ierr = VecScatterEnd (scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1341 /* get global ordering of local dofs */ 1342 ierr = VecGetArray(local_vec,&array);CHKERRQ(ierr); 1343 if (local_dofs_mult) { 1344 for (i=0;i<n_local_dofs;i++) { 1345 temp_global_dofs[i] = (PetscInt)PetscRealPart(array[local_dofs[i]])-local_dofs_mult[i]; 1346 } 1347 } else { 1348 for (i=0;i<n_local_dofs;i++) { 1349 temp_global_dofs[i] = (PetscInt)PetscRealPart(array[local_dofs[i]])-1; 1350 } 1351 } 1352 ierr = VecRestoreArray(local_vec,&array);CHKERRQ(ierr); 1353 /* free workspace */ 1354 ierr = VecScatterDestroy(&scatter_ctx);CHKERRQ(ierr); 1355 ierr = VecDestroy(&local_vec);CHKERRQ(ierr); 1356 ierr = VecDestroy(&global_vec);CHKERRQ(ierr); 1357 ierr = PetscFree(dof_sizes);CHKERRQ(ierr); 1358 ierr = PetscFree(dof_displs);CHKERRQ(ierr); 1359 /* return pointer to global ordering of local dofs */ 1360 *global_numbering_subset = temp_global_dofs; 1361 PetscFunctionReturn(0); 1362 } 1363