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