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