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