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