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