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