1 #include <petsc/private/pcbddcimpl.h> 2 #include <petsc/private/pcbddcprivateimpl.h> 3 #include <petscblaslapack.h> 4 5 static PetscErrorCode MatMult_BDdelta_deluxe_nonred(Mat A, Vec x, Vec y) 6 { 7 BDdelta_DN ctx; 8 9 PetscFunctionBegin; 10 PetscCall(MatShellGetContext(A,&ctx)); 11 PetscCall(MatMultTranspose(ctx->BD,x,ctx->work)); 12 PetscCall(KSPSolveTranspose(ctx->kBD,ctx->work,y)); 13 /* No PC so cannot propagate up failure in KSPSolveTranspose() */ 14 PetscFunctionReturn(0); 15 } 16 17 static PetscErrorCode MatMultTranspose_BDdelta_deluxe_nonred(Mat A, Vec x, Vec y) 18 { 19 BDdelta_DN ctx; 20 21 PetscFunctionBegin; 22 PetscCall(MatShellGetContext(A,&ctx)); 23 PetscCall(KSPSolve(ctx->kBD,x,ctx->work)); 24 /* No PC so cannot propagate up failure in KSPSolve() */ 25 PetscCall(MatMult(ctx->BD,ctx->work,y)); 26 PetscFunctionReturn(0); 27 } 28 29 static PetscErrorCode MatDestroy_BDdelta_deluxe_nonred(Mat A) 30 { 31 BDdelta_DN ctx; 32 33 PetscFunctionBegin; 34 PetscCall(MatShellGetContext(A,&ctx)); 35 PetscCall(MatDestroy(&ctx->BD)); 36 PetscCall(KSPDestroy(&ctx->kBD)); 37 PetscCall(VecDestroy(&ctx->work)); 38 PetscCall(PetscFree(ctx)); 39 PetscFunctionReturn(0); 40 } 41 42 PetscErrorCode PCBDDCCreateFETIDPMatContext(PC pc, FETIDPMat_ctx *fetidpmat_ctx) 43 { 44 FETIDPMat_ctx newctx; 45 46 PetscFunctionBegin; 47 PetscCall(PetscNew(&newctx)); 48 /* increase the reference count for BDDC preconditioner */ 49 PetscCall(PetscObjectReference((PetscObject)pc)); 50 newctx->pc = pc; 51 *fetidpmat_ctx = newctx; 52 PetscFunctionReturn(0); 53 } 54 55 PetscErrorCode PCBDDCCreateFETIDPPCContext(PC pc, FETIDPPC_ctx *fetidppc_ctx) 56 { 57 FETIDPPC_ctx newctx; 58 59 PetscFunctionBegin; 60 PetscCall(PetscNew(&newctx)); 61 /* increase the reference count for BDDC preconditioner */ 62 PetscCall(PetscObjectReference((PetscObject)pc)); 63 newctx->pc = pc; 64 *fetidppc_ctx = newctx; 65 PetscFunctionReturn(0); 66 } 67 68 PetscErrorCode PCBDDCDestroyFETIDPMat(Mat A) 69 { 70 FETIDPMat_ctx mat_ctx; 71 72 PetscFunctionBegin; 73 PetscCall(MatShellGetContext(A,&mat_ctx)); 74 PetscCall(VecDestroy(&mat_ctx->lambda_local)); 75 PetscCall(VecDestroy(&mat_ctx->temp_solution_D)); 76 PetscCall(VecDestroy(&mat_ctx->temp_solution_B)); 77 PetscCall(MatDestroy(&mat_ctx->B_delta)); 78 PetscCall(MatDestroy(&mat_ctx->B_Ddelta)); 79 PetscCall(MatDestroy(&mat_ctx->B_BB)); 80 PetscCall(MatDestroy(&mat_ctx->B_BI)); 81 PetscCall(MatDestroy(&mat_ctx->Bt_BB)); 82 PetscCall(MatDestroy(&mat_ctx->Bt_BI)); 83 PetscCall(MatDestroy(&mat_ctx->C)); 84 PetscCall(VecDestroy(&mat_ctx->rhs_flip)); 85 PetscCall(VecDestroy(&mat_ctx->vP)); 86 PetscCall(VecDestroy(&mat_ctx->xPg)); 87 PetscCall(VecDestroy(&mat_ctx->yPg)); 88 PetscCall(VecScatterDestroy(&mat_ctx->l2g_lambda)); 89 PetscCall(VecScatterDestroy(&mat_ctx->l2g_lambda_only)); 90 PetscCall(VecScatterDestroy(&mat_ctx->l2g_p)); 91 PetscCall(VecScatterDestroy(&mat_ctx->g2g_p)); 92 PetscCall(PCDestroy(&mat_ctx->pc)); /* decrease PCBDDC reference count */ 93 PetscCall(ISDestroy(&mat_ctx->pressure)); 94 PetscCall(ISDestroy(&mat_ctx->lagrange)); 95 PetscCall(PetscFree(mat_ctx)); 96 PetscFunctionReturn(0); 97 } 98 99 PetscErrorCode PCBDDCDestroyFETIDPPC(PC pc) 100 { 101 FETIDPPC_ctx pc_ctx; 102 103 PetscFunctionBegin; 104 PetscCall(PCShellGetContext(pc,&pc_ctx)); 105 PetscCall(VecDestroy(&pc_ctx->lambda_local)); 106 PetscCall(MatDestroy(&pc_ctx->B_Ddelta)); 107 PetscCall(VecScatterDestroy(&pc_ctx->l2g_lambda)); 108 PetscCall(MatDestroy(&pc_ctx->S_j)); 109 PetscCall(PCDestroy(&pc_ctx->pc)); /* decrease PCBDDC reference count */ 110 PetscCall(VecDestroy(&pc_ctx->xPg)); 111 PetscCall(VecDestroy(&pc_ctx->yPg)); 112 PetscCall(PetscFree(pc_ctx)); 113 PetscFunctionReturn(0); 114 } 115 116 PetscErrorCode PCBDDCSetupFETIDPMatContext(FETIDPMat_ctx fetidpmat_ctx) 117 { 118 PC_IS *pcis=(PC_IS*)fetidpmat_ctx->pc->data; 119 PC_BDDC *pcbddc=(PC_BDDC*)fetidpmat_ctx->pc->data; 120 PCBDDCGraph mat_graph=pcbddc->mat_graph; 121 Mat_IS *matis = (Mat_IS*)fetidpmat_ctx->pc->pmat->data; 122 MPI_Comm comm; 123 Mat ScalingMat,BD1,BD2; 124 Vec fetidp_global; 125 IS IS_l2g_lambda; 126 IS subset,subset_mult,subset_n,isvert; 127 PetscBool skip_node,fully_redundant; 128 PetscInt i,j,k,s,n_boundary_dofs,n_global_lambda,n_vertices,partial_sum; 129 PetscInt cum,n_local_lambda,n_lambda_for_dof,dual_size,n_neg_values,n_pos_values; 130 PetscMPIInt rank,size,buf_size,neigh; 131 PetscScalar scalar_value; 132 const PetscInt *vertex_indices; 133 PetscInt *dual_dofs_boundary_indices,*aux_local_numbering_1; 134 const PetscInt *aux_global_numbering; 135 PetscInt *aux_sums,*cols_B_delta,*l2g_indices; 136 PetscScalar *array,*scaling_factors,*vals_B_delta; 137 PetscScalar **all_factors; 138 PetscInt *aux_local_numbering_2; 139 PetscLayout llay; 140 141 /* saddlepoint */ 142 ISLocalToGlobalMapping l2gmap_p; 143 PetscLayout play; 144 IS gP,pP; 145 PetscInt nPl,nPg,nPgl; 146 147 PetscFunctionBegin; 148 PetscCall(PetscObjectGetComm((PetscObject)(fetidpmat_ctx->pc),&comm)); 149 PetscCallMPI(MPI_Comm_rank(comm,&rank)); 150 PetscCallMPI(MPI_Comm_size(comm,&size)); 151 152 /* saddlepoint */ 153 nPl = 0; 154 nPg = 0; 155 nPgl = 0; 156 gP = NULL; 157 pP = NULL; 158 l2gmap_p = NULL; 159 play = NULL; 160 PetscCall(PetscObjectQuery((PetscObject)fetidpmat_ctx->pc,"__KSPFETIDP_pP",(PetscObject*)&pP)); 161 if (pP) { /* saddle point */ 162 /* subdomain pressures in global numbering */ 163 PetscCall(PetscObjectQuery((PetscObject)fetidpmat_ctx->pc,"__KSPFETIDP_gP",(PetscObject*)&gP)); 164 PetscCheck(gP,PETSC_COMM_SELF,PETSC_ERR_PLIB,"gP not present"); 165 PetscCall(ISGetLocalSize(gP,&nPl)); 166 PetscCall(VecCreate(PETSC_COMM_SELF,&fetidpmat_ctx->vP)); 167 PetscCall(VecSetSizes(fetidpmat_ctx->vP,nPl,nPl)); 168 PetscCall(VecSetType(fetidpmat_ctx->vP,VECSTANDARD)); 169 PetscCall(VecSetUp(fetidpmat_ctx->vP)); 170 171 /* pressure matrix */ 172 PetscCall(PetscObjectQuery((PetscObject)fetidpmat_ctx->pc,"__KSPFETIDP_C",(PetscObject*)&fetidpmat_ctx->C)); 173 if (!fetidpmat_ctx->C) { /* null pressure block, compute layout and global numbering for pressures */ 174 IS Pg; 175 176 PetscCall(ISRenumber(gP,NULL,&nPg,&Pg)); 177 PetscCall(ISLocalToGlobalMappingCreateIS(Pg,&l2gmap_p)); 178 PetscCall(ISDestroy(&Pg)); 179 PetscCall(PetscLayoutCreate(comm,&play)); 180 PetscCall(PetscLayoutSetBlockSize(play,1)); 181 PetscCall(PetscLayoutSetSize(play,nPg)); 182 PetscCall(ISGetLocalSize(pP,&nPgl)); 183 PetscCall(PetscLayoutSetLocalSize(play,nPgl)); 184 PetscCall(PetscLayoutSetUp(play)); 185 } else { 186 PetscCall(PetscObjectReference((PetscObject)fetidpmat_ctx->C)); 187 PetscCall(MatISGetLocalToGlobalMapping(fetidpmat_ctx->C,&l2gmap_p,NULL)); 188 PetscCall(PetscObjectReference((PetscObject)l2gmap_p)); 189 PetscCall(MatGetSize(fetidpmat_ctx->C,&nPg,NULL)); 190 PetscCall(MatGetLocalSize(fetidpmat_ctx->C,NULL,&nPgl)); 191 PetscCall(MatGetLayouts(fetidpmat_ctx->C,NULL,&llay)); 192 PetscCall(PetscLayoutReference(llay,&play)); 193 } 194 PetscCall(VecCreateMPIWithArray(comm,1,nPgl,nPg,NULL,&fetidpmat_ctx->xPg)); 195 PetscCall(VecCreateMPIWithArray(comm,1,nPgl,nPg,NULL,&fetidpmat_ctx->yPg)); 196 197 /* import matrices for pressures coupling */ 198 PetscCall(PetscObjectQuery((PetscObject)fetidpmat_ctx->pc,"__KSPFETIDP_B_BI",(PetscObject*)&fetidpmat_ctx->B_BI)); 199 PetscCheck(fetidpmat_ctx->B_BI,PETSC_COMM_SELF,PETSC_ERR_PLIB,"B_BI not present"); 200 PetscCall(PetscObjectReference((PetscObject)fetidpmat_ctx->B_BI)); 201 202 PetscCall(PetscObjectQuery((PetscObject)fetidpmat_ctx->pc,"__KSPFETIDP_B_BB",(PetscObject*)&fetidpmat_ctx->B_BB)); 203 PetscCheck(fetidpmat_ctx->B_BB,PETSC_COMM_SELF,PETSC_ERR_PLIB,"B_BB not present"); 204 PetscCall(PetscObjectReference((PetscObject)fetidpmat_ctx->B_BB)); 205 206 PetscCall(PetscObjectQuery((PetscObject)fetidpmat_ctx->pc,"__KSPFETIDP_Bt_BI",(PetscObject*)&fetidpmat_ctx->Bt_BI)); 207 PetscCheck(fetidpmat_ctx->Bt_BI,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Bt_BI not present"); 208 PetscCall(PetscObjectReference((PetscObject)fetidpmat_ctx->Bt_BI)); 209 210 PetscCall(PetscObjectQuery((PetscObject)fetidpmat_ctx->pc,"__KSPFETIDP_Bt_BB",(PetscObject*)&fetidpmat_ctx->Bt_BB)); 211 PetscCheck(fetidpmat_ctx->Bt_BB,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Bt_BB not present"); 212 PetscCall(PetscObjectReference((PetscObject)fetidpmat_ctx->Bt_BB)); 213 214 PetscCall(PetscObjectQuery((PetscObject)fetidpmat_ctx->pc,"__KSPFETIDP_flip" ,(PetscObject*)&fetidpmat_ctx->rhs_flip)); 215 if (fetidpmat_ctx->rhs_flip) PetscCall(PetscObjectReference((PetscObject)fetidpmat_ctx->rhs_flip)); 216 } 217 218 /* Default type of lagrange multipliers is non-redundant */ 219 fully_redundant = fetidpmat_ctx->fully_redundant; 220 221 /* Evaluate local and global number of lagrange multipliers */ 222 PetscCall(VecSet(pcis->vec1_N,0.0)); 223 n_local_lambda = 0; 224 partial_sum = 0; 225 n_boundary_dofs = 0; 226 s = 0; 227 228 /* Get Vertices used to define the BDDC */ 229 PetscCall(PCBDDCGraphGetCandidatesIS(pcbddc->mat_graph,NULL,NULL,NULL,NULL,&isvert)); 230 PetscCall(ISGetLocalSize(isvert,&n_vertices)); 231 PetscCall(ISGetIndices(isvert,&vertex_indices)); 232 233 dual_size = pcis->n_B-n_vertices; 234 PetscCall(PetscMalloc1(dual_size,&dual_dofs_boundary_indices)); 235 PetscCall(PetscMalloc1(dual_size,&aux_local_numbering_1)); 236 PetscCall(PetscMalloc1(dual_size,&aux_local_numbering_2)); 237 238 PetscCall(VecGetArray(pcis->vec1_N,&array)); 239 for (i=0;i<pcis->n;i++) { 240 j = mat_graph->count[i]; /* RECALL: mat_graph->count[i] does not count myself */ 241 if (j > 0) n_boundary_dofs++; 242 skip_node = PETSC_FALSE; 243 if (s < n_vertices && vertex_indices[s] == i) { /* it works for a sorted set of vertices */ 244 skip_node = PETSC_TRUE; 245 s++; 246 } 247 if (j < 1) skip_node = PETSC_TRUE; 248 if (mat_graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK) skip_node = PETSC_TRUE; 249 if (!skip_node) { 250 if (fully_redundant) { 251 /* fully redundant set of lagrange multipliers */ 252 n_lambda_for_dof = (j*(j+1))/2; 253 } else { 254 n_lambda_for_dof = j; 255 } 256 n_local_lambda += j; 257 /* needed to evaluate global number of lagrange multipliers */ 258 array[i]=(1.0*n_lambda_for_dof)/(j+1.0); /* already scaled for the next global sum */ 259 /* store some data needed */ 260 dual_dofs_boundary_indices[partial_sum] = n_boundary_dofs-1; 261 aux_local_numbering_1[partial_sum] = i; 262 aux_local_numbering_2[partial_sum] = n_lambda_for_dof; 263 partial_sum++; 264 } 265 } 266 PetscCall(VecRestoreArray(pcis->vec1_N,&array)); 267 PetscCall(ISRestoreIndices(isvert,&vertex_indices)); 268 PetscCall(PCBDDCGraphRestoreCandidatesIS(pcbddc->mat_graph,NULL,NULL,NULL,NULL,&isvert)); 269 dual_size = partial_sum; 270 271 /* compute global ordering of lagrange multipliers and associate l2g map */ 272 PetscCall(ISCreateGeneral(comm,partial_sum,aux_local_numbering_1,PETSC_COPY_VALUES,&subset_n)); 273 PetscCall(ISLocalToGlobalMappingApplyIS(pcis->mapping,subset_n,&subset)); 274 PetscCall(ISDestroy(&subset_n)); 275 PetscCall(ISCreateGeneral(comm,partial_sum,aux_local_numbering_2,PETSC_OWN_POINTER,&subset_mult)); 276 PetscCall(ISRenumber(subset,subset_mult,&fetidpmat_ctx->n_lambda,&subset_n)); 277 PetscCall(ISDestroy(&subset)); 278 279 if (PetscDefined(USE_DEBUG)) { 280 PetscCall(VecSet(pcis->vec1_global,0.0)); 281 PetscCall(VecScatterBegin(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE)); 282 PetscCall(VecScatterEnd(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE)); 283 PetscCall(VecSum(pcis->vec1_global,&scalar_value)); 284 i = (PetscInt)PetscRealPart(scalar_value); 285 PetscCheck(i == fetidpmat_ctx->n_lambda,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Global number of multipliers mismatch! (%" PetscInt_FMT " != %" PetscInt_FMT ")",fetidpmat_ctx->n_lambda,i); 286 } 287 288 /* init data for scaling factors exchange */ 289 if (!pcbddc->use_deluxe_scaling) { 290 PetscInt *ptrs_buffer,neigh_position; 291 PetscScalar *send_buffer,*recv_buffer; 292 MPI_Request *send_reqs,*recv_reqs; 293 294 partial_sum = 0; 295 PetscCall(PetscMalloc1(pcis->n_neigh,&ptrs_buffer)); 296 PetscCall(PetscMalloc1(PetscMax(pcis->n_neigh-1,0),&send_reqs)); 297 PetscCall(PetscMalloc1(PetscMax(pcis->n_neigh-1,0),&recv_reqs)); 298 PetscCall(PetscMalloc1(pcis->n+1,&all_factors)); 299 if (pcis->n_neigh > 0) ptrs_buffer[0]=0; 300 for (i=1;i<pcis->n_neigh;i++) { 301 partial_sum += pcis->n_shared[i]; 302 ptrs_buffer[i] = ptrs_buffer[i-1]+pcis->n_shared[i]; 303 } 304 PetscCall(PetscMalloc1(partial_sum,&send_buffer)); 305 PetscCall(PetscMalloc1(partial_sum,&recv_buffer)); 306 PetscCall(PetscMalloc1(partial_sum,&all_factors[0])); 307 for (i=0;i<pcis->n-1;i++) { 308 j = mat_graph->count[i]; 309 all_factors[i+1]=all_factors[i]+j; 310 } 311 312 /* scatter B scaling to N vec */ 313 PetscCall(VecScatterBegin(pcis->N_to_B,pcis->D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE)); 314 PetscCall(VecScatterEnd(pcis->N_to_B,pcis->D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE)); 315 /* communications */ 316 PetscCall(VecGetArrayRead(pcis->vec1_N,(const PetscScalar**)&array)); 317 for (i=1;i<pcis->n_neigh;i++) { 318 for (j=0;j<pcis->n_shared[i];j++) { 319 send_buffer[ptrs_buffer[i-1]+j]=array[pcis->shared[i][j]]; 320 } 321 PetscCall(PetscMPIIntCast(ptrs_buffer[i]-ptrs_buffer[i-1],&buf_size)); 322 PetscCall(PetscMPIIntCast(pcis->neigh[i],&neigh)); 323 PetscCallMPI(MPI_Isend(&send_buffer[ptrs_buffer[i-1]],buf_size,MPIU_SCALAR,neigh,0,comm,&send_reqs[i-1])); 324 PetscCallMPI(MPI_Irecv(&recv_buffer[ptrs_buffer[i-1]],buf_size,MPIU_SCALAR,neigh,0,comm,&recv_reqs[i-1])); 325 } 326 PetscCall(VecRestoreArrayRead(pcis->vec1_N,(const PetscScalar**)&array)); 327 if (pcis->n_neigh > 0) { 328 PetscCallMPI(MPI_Waitall(pcis->n_neigh-1,recv_reqs,MPI_STATUSES_IGNORE)); 329 } 330 /* put values in correct places */ 331 for (i=1;i<pcis->n_neigh;i++) { 332 for (j=0;j<pcis->n_shared[i];j++) { 333 k = pcis->shared[i][j]; 334 neigh_position = 0; 335 while (mat_graph->neighbours_set[k][neigh_position] != pcis->neigh[i]) {neigh_position++;} 336 all_factors[k][neigh_position]=recv_buffer[ptrs_buffer[i-1]+j]; 337 } 338 } 339 if (pcis->n_neigh > 0) { 340 PetscCallMPI(MPI_Waitall(pcis->n_neigh-1,send_reqs,MPI_STATUSES_IGNORE)); 341 } 342 PetscCall(PetscFree(send_reqs)); 343 PetscCall(PetscFree(recv_reqs)); 344 PetscCall(PetscFree(send_buffer)); 345 PetscCall(PetscFree(recv_buffer)); 346 PetscCall(PetscFree(ptrs_buffer)); 347 } 348 349 /* Compute B and B_delta (local actions) */ 350 PetscCall(PetscMalloc1(pcis->n_neigh,&aux_sums)); 351 PetscCall(PetscMalloc1(n_local_lambda,&l2g_indices)); 352 PetscCall(PetscMalloc1(n_local_lambda,&vals_B_delta)); 353 PetscCall(PetscMalloc1(n_local_lambda,&cols_B_delta)); 354 if (!pcbddc->use_deluxe_scaling) { 355 PetscCall(PetscMalloc1(n_local_lambda,&scaling_factors)); 356 } else { 357 scaling_factors = NULL; 358 all_factors = NULL; 359 } 360 PetscCall(ISGetIndices(subset_n,&aux_global_numbering)); 361 partial_sum=0; 362 cum = 0; 363 for (i=0;i<dual_size;i++) { 364 n_global_lambda = aux_global_numbering[cum]; 365 j = mat_graph->count[aux_local_numbering_1[i]]; 366 aux_sums[0]=0; 367 for (s=1;s<j;s++) { 368 aux_sums[s]=aux_sums[s-1]+j-s+1; 369 } 370 if (all_factors) array = all_factors[aux_local_numbering_1[i]]; 371 n_neg_values = 0; 372 while (n_neg_values < j && mat_graph->neighbours_set[aux_local_numbering_1[i]][n_neg_values] < rank) {n_neg_values++;} 373 n_pos_values = j - n_neg_values; 374 if (fully_redundant) { 375 for (s=0;s<n_neg_values;s++) { 376 l2g_indices [partial_sum+s]=aux_sums[s]+n_neg_values-s-1+n_global_lambda; 377 cols_B_delta [partial_sum+s]=dual_dofs_boundary_indices[i]; 378 vals_B_delta [partial_sum+s]=-1.0; 379 if (!pcbddc->use_deluxe_scaling) scaling_factors[partial_sum+s]=array[s]; 380 } 381 for (s=0;s<n_pos_values;s++) { 382 l2g_indices [partial_sum+s+n_neg_values]=aux_sums[n_neg_values]+s+n_global_lambda; 383 cols_B_delta [partial_sum+s+n_neg_values]=dual_dofs_boundary_indices[i]; 384 vals_B_delta [partial_sum+s+n_neg_values]=1.0; 385 if (!pcbddc->use_deluxe_scaling) scaling_factors[partial_sum+s+n_neg_values]=array[s+n_neg_values]; 386 } 387 partial_sum += j; 388 } else { 389 /* l2g_indices and default cols and vals of B_delta */ 390 for (s=0;s<j;s++) { 391 l2g_indices [partial_sum+s]=n_global_lambda+s; 392 cols_B_delta [partial_sum+s]=dual_dofs_boundary_indices[i]; 393 vals_B_delta [partial_sum+s]=0.0; 394 } 395 /* B_delta */ 396 if (n_neg_values > 0) { /* there's a rank next to me to the left */ 397 vals_B_delta [partial_sum+n_neg_values-1]=-1.0; 398 } 399 if (n_neg_values < j) { /* there's a rank next to me to the right */ 400 vals_B_delta [partial_sum+n_neg_values]=1.0; 401 } 402 /* scaling as in Klawonn-Widlund 1999 */ 403 if (!pcbddc->use_deluxe_scaling) { 404 for (s=0;s<n_neg_values;s++) { 405 scalar_value = 0.0; 406 for (k=0;k<s+1;k++) scalar_value += array[k]; 407 scaling_factors[partial_sum+s] = -scalar_value; 408 } 409 for (s=0;s<n_pos_values;s++) { 410 scalar_value = 0.0; 411 for (k=s+n_neg_values;k<j;k++) scalar_value += array[k]; 412 scaling_factors[partial_sum+s+n_neg_values] = scalar_value; 413 } 414 } 415 partial_sum += j; 416 } 417 cum += aux_local_numbering_2[i]; 418 } 419 PetscCall(ISRestoreIndices(subset_n,&aux_global_numbering)); 420 PetscCall(ISDestroy(&subset_mult)); 421 PetscCall(ISDestroy(&subset_n)); 422 PetscCall(PetscFree(aux_sums)); 423 PetscCall(PetscFree(aux_local_numbering_1)); 424 PetscCall(PetscFree(dual_dofs_boundary_indices)); 425 if (all_factors) { 426 PetscCall(PetscFree(all_factors[0])); 427 PetscCall(PetscFree(all_factors)); 428 } 429 430 /* Create local part of B_delta */ 431 PetscCall(MatCreate(PETSC_COMM_SELF,&fetidpmat_ctx->B_delta)); 432 PetscCall(MatSetSizes(fetidpmat_ctx->B_delta,n_local_lambda,pcis->n_B,n_local_lambda,pcis->n_B)); 433 PetscCall(MatSetType(fetidpmat_ctx->B_delta,MATSEQAIJ)); 434 PetscCall(MatSeqAIJSetPreallocation(fetidpmat_ctx->B_delta,1,NULL)); 435 PetscCall(MatSetOption(fetidpmat_ctx->B_delta,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE)); 436 for (i=0;i<n_local_lambda;i++) { 437 PetscCall(MatSetValue(fetidpmat_ctx->B_delta,i,cols_B_delta[i],vals_B_delta[i],INSERT_VALUES)); 438 } 439 PetscCall(PetscFree(vals_B_delta)); 440 PetscCall(MatAssemblyBegin(fetidpmat_ctx->B_delta,MAT_FINAL_ASSEMBLY)); 441 PetscCall(MatAssemblyEnd(fetidpmat_ctx->B_delta,MAT_FINAL_ASSEMBLY)); 442 443 BD1 = NULL; 444 BD2 = NULL; 445 if (fully_redundant) { 446 PetscCheck(!pcbddc->use_deluxe_scaling,comm,PETSC_ERR_SUP,"Deluxe FETIDP with fully-redundant multipliers to be implemented"); 447 PetscCall(MatCreate(PETSC_COMM_SELF,&ScalingMat)); 448 PetscCall(MatSetSizes(ScalingMat,n_local_lambda,n_local_lambda,n_local_lambda,n_local_lambda)); 449 PetscCall(MatSetType(ScalingMat,MATSEQAIJ)); 450 PetscCall(MatSeqAIJSetPreallocation(ScalingMat,1,NULL)); 451 for (i=0;i<n_local_lambda;i++) { 452 PetscCall(MatSetValue(ScalingMat,i,i,scaling_factors[i],INSERT_VALUES)); 453 } 454 PetscCall(MatAssemblyBegin(ScalingMat,MAT_FINAL_ASSEMBLY)); 455 PetscCall(MatAssemblyEnd(ScalingMat,MAT_FINAL_ASSEMBLY)); 456 PetscCall(MatMatMult(ScalingMat,fetidpmat_ctx->B_delta,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&fetidpmat_ctx->B_Ddelta)); 457 PetscCall(MatDestroy(&ScalingMat)); 458 } else { 459 PetscCall(MatCreate(PETSC_COMM_SELF,&fetidpmat_ctx->B_Ddelta)); 460 PetscCall(MatSetSizes(fetidpmat_ctx->B_Ddelta,n_local_lambda,pcis->n_B,n_local_lambda,pcis->n_B)); 461 if (!pcbddc->use_deluxe_scaling || !pcbddc->sub_schurs) { 462 PetscCall(MatSetType(fetidpmat_ctx->B_Ddelta,MATSEQAIJ)); 463 PetscCall(MatSeqAIJSetPreallocation(fetidpmat_ctx->B_Ddelta,1,NULL)); 464 for (i=0;i<n_local_lambda;i++) { 465 PetscCall(MatSetValue(fetidpmat_ctx->B_Ddelta,i,cols_B_delta[i],scaling_factors[i],INSERT_VALUES)); 466 } 467 PetscCall(MatAssemblyBegin(fetidpmat_ctx->B_Ddelta,MAT_FINAL_ASSEMBLY)); 468 PetscCall(MatAssemblyEnd(fetidpmat_ctx->B_Ddelta,MAT_FINAL_ASSEMBLY)); 469 } else { 470 /* scaling as in Klawonn-Widlund 1999 */ 471 PCBDDCDeluxeScaling deluxe_ctx = pcbddc->deluxe_ctx; 472 PCBDDCSubSchurs sub_schurs = pcbddc->sub_schurs; 473 Mat T; 474 PetscScalar *W,lwork,*Bwork; 475 const PetscInt *idxs = NULL; 476 PetscInt cum,mss,*nnz; 477 PetscBLASInt *pivots,B_lwork,B_N,B_ierr; 478 479 PetscCheck(pcbddc->deluxe_singlemat,comm,PETSC_ERR_USER,"Cannot compute B_Ddelta! rerun with -pc_bddc_deluxe_singlemat"); 480 mss = 0; 481 PetscCall(PetscCalloc1(pcis->n_B,&nnz)); 482 if (sub_schurs->is_Ej_all) { 483 PetscCall(ISGetIndices(sub_schurs->is_Ej_all,&idxs)); 484 for (i=0,cum=0;i<sub_schurs->n_subs;i++) { 485 PetscInt subset_size; 486 487 PetscCall(ISGetLocalSize(sub_schurs->is_subs[i],&subset_size)); 488 for (j=0;j<subset_size;j++) nnz[idxs[j+cum]] = subset_size; 489 mss = PetscMax(mss,subset_size); 490 cum += subset_size; 491 } 492 } 493 PetscCall(MatCreate(PETSC_COMM_SELF,&T)); 494 PetscCall(MatSetSizes(T,pcis->n_B,pcis->n_B,pcis->n_B,pcis->n_B)); 495 PetscCall(MatSetType(T,MATSEQAIJ)); 496 PetscCall(MatSeqAIJSetPreallocation(T,0,nnz)); 497 PetscCall(PetscFree(nnz)); 498 499 /* workspace allocation */ 500 B_lwork = 0; 501 if (mss) { 502 PetscScalar dummy = 1; 503 504 B_lwork = -1; 505 PetscCall(PetscBLASIntCast(mss,&B_N)); 506 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 507 PetscCallBLAS("LAPACKgetri",LAPACKgetri_(&B_N,&dummy,&B_N,&B_N,&lwork,&B_lwork,&B_ierr)); 508 PetscCall(PetscFPTrapPop()); 509 PetscCheck(!B_ierr,PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in query to GETRI Lapack routine %d",(int)B_ierr); 510 PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(lwork),&B_lwork)); 511 } 512 PetscCall(PetscMalloc3(mss*mss,&W,mss,&pivots,B_lwork,&Bwork)); 513 514 for (i=0,cum=0;i<sub_schurs->n_subs;i++) { 515 const PetscScalar *M; 516 PetscInt subset_size; 517 518 PetscCall(ISGetLocalSize(sub_schurs->is_subs[i],&subset_size)); 519 PetscCall(PetscBLASIntCast(subset_size,&B_N)); 520 PetscCall(MatDenseGetArrayRead(deluxe_ctx->seq_mat[i],&M)); 521 PetscCall(PetscArraycpy(W,M,subset_size*subset_size)); 522 PetscCall(MatDenseRestoreArrayRead(deluxe_ctx->seq_mat[i],&M)); 523 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 524 PetscCallBLAS("LAPACKgetrf",LAPACKgetrf_(&B_N,&B_N,W,&B_N,pivots,&B_ierr)); 525 PetscCheck(!B_ierr,PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRF Lapack routine %d",(int)B_ierr); 526 PetscCallBLAS("LAPACKgetri",LAPACKgetri_(&B_N,W,&B_N,pivots,Bwork,&B_lwork,&B_ierr)); 527 PetscCheck(!B_ierr,PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRI Lapack routine %d",(int)B_ierr); 528 PetscCall(PetscFPTrapPop()); 529 /* silent static analyzer */ 530 PetscCheck(idxs,PETSC_COMM_SELF,PETSC_ERR_PLIB,"IDXS not present"); 531 PetscCall(MatSetValues(T,subset_size,idxs+cum,subset_size,idxs+cum,W,INSERT_VALUES)); 532 cum += subset_size; 533 } 534 PetscCall(MatAssemblyBegin(T,MAT_FINAL_ASSEMBLY)); 535 PetscCall(MatAssemblyEnd(T,MAT_FINAL_ASSEMBLY)); 536 PetscCall(MatMatTransposeMult(T,fetidpmat_ctx->B_delta,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&BD1)); 537 PetscCall(MatMatMult(fetidpmat_ctx->B_delta,BD1,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&BD2)); 538 PetscCall(MatDestroy(&T)); 539 PetscCall(PetscFree3(W,pivots,Bwork)); 540 if (sub_schurs->is_Ej_all) { 541 PetscCall(ISRestoreIndices(sub_schurs->is_Ej_all,&idxs)); 542 } 543 } 544 } 545 PetscCall(PetscFree(scaling_factors)); 546 PetscCall(PetscFree(cols_B_delta)); 547 548 /* Layout of multipliers */ 549 PetscCall(PetscLayoutCreate(comm,&llay)); 550 PetscCall(PetscLayoutSetBlockSize(llay,1)); 551 PetscCall(PetscLayoutSetSize(llay,fetidpmat_ctx->n_lambda)); 552 PetscCall(PetscLayoutSetUp(llay)); 553 PetscCall(PetscLayoutGetLocalSize(llay,&fetidpmat_ctx->n)); 554 555 /* Local work vector of multipliers */ 556 PetscCall(VecCreate(PETSC_COMM_SELF,&fetidpmat_ctx->lambda_local)); 557 PetscCall(VecSetSizes(fetidpmat_ctx->lambda_local,n_local_lambda,n_local_lambda)); 558 PetscCall(VecSetType(fetidpmat_ctx->lambda_local,VECSEQ)); 559 560 if (BD2) { 561 ISLocalToGlobalMapping l2g; 562 Mat T,TA,*pT; 563 IS is; 564 PetscInt nl,N; 565 BDdelta_DN ctx; 566 567 PetscCall(PetscLayoutGetLocalSize(llay,&nl)); 568 PetscCall(PetscLayoutGetSize(llay,&N)); 569 PetscCall(MatCreate(comm,&T)); 570 PetscCall(MatSetSizes(T,nl,nl,N,N)); 571 PetscCall(MatSetType(T,MATIS)); 572 PetscCall(ISLocalToGlobalMappingCreate(comm,1,n_local_lambda,l2g_indices,PETSC_COPY_VALUES,&l2g)); 573 PetscCall(MatSetLocalToGlobalMapping(T,l2g,l2g)); 574 PetscCall(ISLocalToGlobalMappingDestroy(&l2g)); 575 PetscCall(MatISSetLocalMat(T,BD2)); 576 PetscCall(MatAssemblyBegin(T,MAT_FINAL_ASSEMBLY)); 577 PetscCall(MatAssemblyEnd(T,MAT_FINAL_ASSEMBLY)); 578 PetscCall(MatDestroy(&BD2)); 579 PetscCall(MatConvert(T,MATAIJ,MAT_INITIAL_MATRIX,&TA)); 580 PetscCall(MatDestroy(&T)); 581 PetscCall(ISCreateGeneral(comm,n_local_lambda,l2g_indices,PETSC_USE_POINTER,&is)); 582 PetscCall(MatCreateSubMatrices(TA,1,&is,&is,MAT_INITIAL_MATRIX,&pT)); 583 PetscCall(MatDestroy(&TA)); 584 PetscCall(ISDestroy(&is)); 585 BD2 = pT[0]; 586 PetscCall(PetscFree(pT)); 587 588 /* B_Ddelta for non-redundant multipliers with deluxe scaling */ 589 PetscCall(PetscNew(&ctx)); 590 PetscCall(MatSetType(fetidpmat_ctx->B_Ddelta,MATSHELL)); 591 PetscCall(MatShellSetContext(fetidpmat_ctx->B_Ddelta,ctx)); 592 PetscCall(MatShellSetOperation(fetidpmat_ctx->B_Ddelta,MATOP_MULT,(void (*)(void))MatMult_BDdelta_deluxe_nonred)); 593 PetscCall(MatShellSetOperation(fetidpmat_ctx->B_Ddelta,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_BDdelta_deluxe_nonred)); 594 PetscCall(MatShellSetOperation(fetidpmat_ctx->B_Ddelta,MATOP_DESTROY,(void (*)(void))MatDestroy_BDdelta_deluxe_nonred)); 595 PetscCall(MatSetUp(fetidpmat_ctx->B_Ddelta)); 596 597 PetscCall(PetscObjectReference((PetscObject)BD1)); 598 ctx->BD = BD1; 599 PetscCall(KSPCreate(PETSC_COMM_SELF,&ctx->kBD)); 600 PetscCall(KSPSetOperators(ctx->kBD,BD2,BD2)); 601 PetscCall(VecDuplicate(fetidpmat_ctx->lambda_local,&ctx->work)); 602 fetidpmat_ctx->deluxe_nonred = PETSC_TRUE; 603 } 604 PetscCall(MatDestroy(&BD1)); 605 PetscCall(MatDestroy(&BD2)); 606 607 /* fetidpmat sizes */ 608 fetidpmat_ctx->n += nPgl; 609 fetidpmat_ctx->N = fetidpmat_ctx->n_lambda+nPg; 610 611 /* Global vector for FETI-DP linear system */ 612 PetscCall(VecCreate(comm,&fetidp_global)); 613 PetscCall(VecSetSizes(fetidp_global,fetidpmat_ctx->n,fetidpmat_ctx->N)); 614 PetscCall(VecSetType(fetidp_global,VECMPI)); 615 PetscCall(VecSetUp(fetidp_global)); 616 617 /* Decide layout for fetidp dofs: if it is a saddle point problem 618 pressure is ordered first in the local part of the global vector 619 of the FETI-DP linear system */ 620 if (nPg) { 621 Vec v; 622 IS IS_l2g_p,ais; 623 PetscLayout alay; 624 const PetscInt *idxs,*pranges,*aranges,*lranges; 625 PetscInt *l2g_indices_p,rst; 626 PetscMPIInt rank; 627 628 PetscCall(PetscMalloc1(nPl,&l2g_indices_p)); 629 PetscCall(VecGetLayout(fetidp_global,&alay)); 630 PetscCall(PetscLayoutGetRanges(alay,&aranges)); 631 PetscCall(PetscLayoutGetRanges(play,&pranges)); 632 PetscCall(PetscLayoutGetRanges(llay,&lranges)); 633 634 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)fetidp_global),&rank)); 635 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)fetidp_global),pranges[rank+1]-pranges[rank],aranges[rank],1,&fetidpmat_ctx->pressure)); 636 PetscCall(PetscObjectSetName((PetscObject)fetidpmat_ctx->pressure,"F_P")); 637 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)fetidp_global),lranges[rank+1]-lranges[rank],aranges[rank]+pranges[rank+1]-pranges[rank],1,&fetidpmat_ctx->lagrange)); 638 PetscCall(PetscObjectSetName((PetscObject)fetidpmat_ctx->lagrange,"F_L")); 639 PetscCall(ISLocalToGlobalMappingGetIndices(l2gmap_p,&idxs)); 640 /* shift local to global indices for pressure */ 641 for (i=0;i<nPl;i++) { 642 PetscMPIInt owner; 643 644 PetscCall(PetscLayoutFindOwner(play,idxs[i],&owner)); 645 l2g_indices_p[i] = idxs[i]-pranges[owner]+aranges[owner]; 646 } 647 PetscCall(ISLocalToGlobalMappingRestoreIndices(l2gmap_p,&idxs)); 648 PetscCall(ISCreateGeneral(comm,nPl,l2g_indices_p,PETSC_OWN_POINTER,&IS_l2g_p)); 649 650 /* local to global scatter for pressure */ 651 PetscCall(VecScatterCreate(fetidpmat_ctx->vP,NULL,fetidp_global,IS_l2g_p,&fetidpmat_ctx->l2g_p)); 652 PetscCall(ISDestroy(&IS_l2g_p)); 653 654 /* scatter for lagrange multipliers only */ 655 PetscCall(VecCreate(comm,&v)); 656 PetscCall(VecSetType(v,VECSTANDARD)); 657 PetscCall(VecSetLayout(v,llay)); 658 PetscCall(VecSetUp(v)); 659 PetscCall(ISCreateGeneral(comm,n_local_lambda,l2g_indices,PETSC_COPY_VALUES,&ais)); 660 PetscCall(VecScatterCreate(fetidpmat_ctx->lambda_local,NULL,v,ais,&fetidpmat_ctx->l2g_lambda_only)); 661 PetscCall(ISDestroy(&ais)); 662 PetscCall(VecDestroy(&v)); 663 664 /* shift local to global indices for multipliers */ 665 for (i=0;i<n_local_lambda;i++) { 666 PetscInt ps; 667 PetscMPIInt owner; 668 669 PetscCall(PetscLayoutFindOwner(llay,l2g_indices[i],&owner)); 670 ps = pranges[owner+1]-pranges[owner]; 671 l2g_indices[i] = l2g_indices[i]-lranges[owner]+aranges[owner]+ps; 672 } 673 674 /* scatter from alldofs to pressures global fetidp vector */ 675 PetscCall(PetscLayoutGetRange(alay,&rst,NULL)); 676 PetscCall(ISCreateStride(comm,nPgl,rst,1,&ais)); 677 PetscCall(VecScatterCreate(pcis->vec1_global,pP,fetidp_global,ais,&fetidpmat_ctx->g2g_p)); 678 PetscCall(ISDestroy(&ais)); 679 } 680 PetscCall(PetscLayoutDestroy(&llay)); 681 PetscCall(PetscLayoutDestroy(&play)); 682 PetscCall(ISCreateGeneral(comm,n_local_lambda,l2g_indices,PETSC_OWN_POINTER,&IS_l2g_lambda)); 683 684 /* scatter from local to global multipliers */ 685 PetscCall(VecScatterCreate(fetidpmat_ctx->lambda_local,NULL,fetidp_global,IS_l2g_lambda,&fetidpmat_ctx->l2g_lambda)); 686 PetscCall(ISDestroy(&IS_l2g_lambda)); 687 PetscCall(ISLocalToGlobalMappingDestroy(&l2gmap_p)); 688 PetscCall(VecDestroy(&fetidp_global)); 689 690 /* Create some work vectors needed by fetidp */ 691 PetscCall(VecDuplicate(pcis->vec1_B,&fetidpmat_ctx->temp_solution_B)); 692 PetscCall(VecDuplicate(pcis->vec1_D,&fetidpmat_ctx->temp_solution_D)); 693 PetscFunctionReturn(0); 694 } 695 696 PetscErrorCode PCBDDCSetupFETIDPPCContext(Mat fetimat, FETIDPPC_ctx fetidppc_ctx) 697 { 698 FETIDPMat_ctx mat_ctx; 699 PC_BDDC *pcbddc = (PC_BDDC*)fetidppc_ctx->pc->data; 700 PC_IS *pcis = (PC_IS*)fetidppc_ctx->pc->data; 701 PetscBool lumped = PETSC_FALSE; 702 703 PetscFunctionBegin; 704 PetscCall(MatShellGetContext(fetimat,&mat_ctx)); 705 /* get references from objects created when setting up feti mat context */ 706 PetscCall(PetscObjectReference((PetscObject)mat_ctx->lambda_local)); 707 fetidppc_ctx->lambda_local = mat_ctx->lambda_local; 708 PetscCall(PetscObjectReference((PetscObject)mat_ctx->B_Ddelta)); 709 fetidppc_ctx->B_Ddelta = mat_ctx->B_Ddelta; 710 if (mat_ctx->deluxe_nonred) { 711 PC pc,mpc; 712 BDdelta_DN ctx; 713 MatSolverType solver; 714 const char *prefix; 715 716 PetscCall(MatShellGetContext(mat_ctx->B_Ddelta,&ctx)); 717 PetscCall(KSPSetType(ctx->kBD,KSPPREONLY)); 718 PetscCall(KSPGetPC(ctx->kBD,&mpc)); 719 PetscCall(KSPGetPC(pcbddc->ksp_D,&pc)); 720 PetscCall(PCSetType(mpc,PCLU)); 721 PetscCall(PCFactorGetMatSolverType(pc,(MatSolverType*)&solver)); 722 if (solver) PetscCall(PCFactorSetMatSolverType(mpc,solver)); 723 PetscCall(MatGetOptionsPrefix(fetimat,&prefix)); 724 PetscCall(KSPSetOptionsPrefix(ctx->kBD,prefix)); 725 PetscCall(KSPAppendOptionsPrefix(ctx->kBD,"bddelta_")); 726 PetscCall(KSPSetFromOptions(ctx->kBD)); 727 } 728 729 if (mat_ctx->l2g_lambda_only) { 730 PetscCall(PetscObjectReference((PetscObject)mat_ctx->l2g_lambda_only)); 731 fetidppc_ctx->l2g_lambda = mat_ctx->l2g_lambda_only; 732 } else { 733 PetscCall(PetscObjectReference((PetscObject)mat_ctx->l2g_lambda)); 734 fetidppc_ctx->l2g_lambda = mat_ctx->l2g_lambda; 735 } 736 /* Dirichlet preconditioner */ 737 PetscCall(PetscOptionsGetBool(NULL,((PetscObject)fetimat)->prefix,"-pc_lumped",&lumped,NULL)); 738 if (!lumped) { 739 IS iV; 740 PetscBool discrete_harmonic = PETSC_FALSE; 741 742 PetscCall(PetscObjectQuery((PetscObject)fetidppc_ctx->pc,"__KSPFETIDP_iV",(PetscObject*)&iV)); 743 if (iV) { 744 PetscCall(PetscOptionsGetBool(NULL,((PetscObject)fetimat)->prefix,"-pc_discrete_harmonic",&discrete_harmonic,NULL)); 745 } 746 if (discrete_harmonic) { 747 KSP sksp; 748 PC pc; 749 PCBDDCSubSchurs sub_schurs = pcbddc->sub_schurs; 750 Mat A_II,A_IB,A_BI; 751 IS iP = NULL; 752 PetscBool isshell,reuse = PETSC_FALSE; 753 KSPType ksptype; 754 const char *prefix; 755 756 /* 757 We constructs a Schur complement for 758 759 | A_II A_ID | 760 | A_DI A_DD | 761 762 instead of 763 764 | A_II B^t_II A_ID | 765 | B_II -C_II B_ID | 766 | A_DI B^t_ID A_DD | 767 768 */ 769 if (sub_schurs && sub_schurs->reuse_solver) { 770 PetscCall(PetscObjectQuery((PetscObject)sub_schurs->A,"__KSPFETIDP_iP",(PetscObject*)&iP)); 771 if (iP) reuse = PETSC_TRUE; 772 } 773 if (!reuse) { 774 IS aB; 775 PetscInt nb; 776 PetscCall(ISGetLocalSize(pcis->is_B_local,&nb)); 777 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pcis->A_II),nb,0,1,&aB)); 778 PetscCall(MatCreateSubMatrix(pcis->A_II,iV,iV,MAT_INITIAL_MATRIX,&A_II)); 779 PetscCall(MatCreateSubMatrix(pcis->A_IB,iV,aB,MAT_INITIAL_MATRIX,&A_IB)); 780 PetscCall(MatCreateSubMatrix(pcis->A_BI,aB,iV,MAT_INITIAL_MATRIX,&A_BI)); 781 PetscCall(ISDestroy(&aB)); 782 } else { 783 PetscCall(MatCreateSubMatrix(sub_schurs->A,pcis->is_I_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&A_IB)); 784 PetscCall(MatCreateSubMatrix(sub_schurs->A,pcis->is_B_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&A_BI)); 785 PetscCall(PetscObjectReference((PetscObject)pcis->A_II)); 786 A_II = pcis->A_II; 787 } 788 PetscCall(MatCreateSchurComplement(A_II,A_II,A_IB,A_BI,pcis->A_BB,&fetidppc_ctx->S_j)); 789 790 /* propagate settings of solver */ 791 PetscCall(MatSchurComplementGetKSP(fetidppc_ctx->S_j,&sksp)); 792 PetscCall(KSPGetType(pcis->ksp_D,&ksptype)); 793 PetscCall(KSPSetType(sksp,ksptype)); 794 PetscCall(KSPGetPC(pcis->ksp_D,&pc)); 795 PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCSHELL,&isshell)); 796 if (!isshell) { 797 MatSolverType solver; 798 PCType pctype; 799 800 PetscCall(PCGetType(pc,&pctype)); 801 PetscCall(PCFactorGetMatSolverType(pc,(MatSolverType*)&solver)); 802 PetscCall(KSPGetPC(sksp,&pc)); 803 PetscCall(PCSetType(pc,pctype)); 804 if (solver) PetscCall(PCFactorSetMatSolverType(pc,solver)); 805 } else { 806 PetscCall(KSPGetPC(sksp,&pc)); 807 PetscCall(PCSetType(pc,PCLU)); 808 } 809 PetscCall(MatDestroy(&A_II)); 810 PetscCall(MatDestroy(&A_IB)); 811 PetscCall(MatDestroy(&A_BI)); 812 PetscCall(MatGetOptionsPrefix(fetimat,&prefix)); 813 PetscCall(KSPSetOptionsPrefix(sksp,prefix)); 814 PetscCall(KSPAppendOptionsPrefix(sksp,"harmonic_")); 815 PetscCall(KSPSetFromOptions(sksp)); 816 if (reuse) { 817 PetscCall(KSPSetPC(sksp,sub_schurs->reuse_solver->interior_solver)); 818 PetscCall(PetscObjectIncrementTabLevel((PetscObject)sub_schurs->reuse_solver->interior_solver,(PetscObject)sksp,0)); 819 } 820 } else { /* default Dirichlet preconditioner is pde-harmonic */ 821 PetscCall(MatCreateSchurComplement(pcis->A_II,pcis->A_II,pcis->A_IB,pcis->A_BI,pcis->A_BB,&fetidppc_ctx->S_j)); 822 PetscCall(MatSchurComplementSetKSP(fetidppc_ctx->S_j,pcis->ksp_D)); 823 } 824 } else { 825 PetscCall(PetscObjectReference((PetscObject)pcis->A_BB)); 826 fetidppc_ctx->S_j = pcis->A_BB; 827 } 828 /* saddle-point */ 829 if (mat_ctx->xPg) { 830 PetscCall(PetscObjectReference((PetscObject)mat_ctx->xPg)); 831 fetidppc_ctx->xPg = mat_ctx->xPg; 832 PetscCall(PetscObjectReference((PetscObject)mat_ctx->yPg)); 833 fetidppc_ctx->yPg = mat_ctx->yPg; 834 } 835 PetscFunctionReturn(0); 836 } 837 838 PetscErrorCode FETIDPMatMult_Kernel(Mat fetimat, Vec x, Vec y, PetscBool trans) 839 { 840 FETIDPMat_ctx mat_ctx; 841 PC_BDDC *pcbddc; 842 PC_IS *pcis; 843 844 PetscFunctionBegin; 845 PetscCall(MatShellGetContext(fetimat,&mat_ctx)); 846 pcis = (PC_IS*)mat_ctx->pc->data; 847 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 848 /* Application of B_delta^T */ 849 PetscCall(VecSet(pcis->vec1_B,0.)); 850 PetscCall(VecScatterBegin(mat_ctx->l2g_lambda,x,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE)); 851 PetscCall(VecScatterEnd(mat_ctx->l2g_lambda,x,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE)); 852 PetscCall(MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B)); 853 854 /* Add contribution from saddle point */ 855 if (mat_ctx->l2g_p) { 856 PetscCall(VecScatterBegin(mat_ctx->l2g_p,x,mat_ctx->vP,INSERT_VALUES,SCATTER_REVERSE)); 857 PetscCall(VecScatterEnd(mat_ctx->l2g_p,x,mat_ctx->vP,INSERT_VALUES,SCATTER_REVERSE)); 858 if (pcbddc->switch_static) { 859 if (trans) { 860 PetscCall(MatMultTranspose(mat_ctx->B_BI,mat_ctx->vP,pcis->vec1_D)); 861 } else { 862 PetscCall(MatMult(mat_ctx->Bt_BI,mat_ctx->vP,pcis->vec1_D)); 863 } 864 } 865 if (trans) { 866 PetscCall(MatMultTransposeAdd(mat_ctx->B_BB,mat_ctx->vP,pcis->vec1_B,pcis->vec1_B)); 867 } else { 868 PetscCall(MatMultAdd(mat_ctx->Bt_BB,mat_ctx->vP,pcis->vec1_B,pcis->vec1_B)); 869 } 870 } else { 871 if (pcbddc->switch_static) PetscCall(VecSet(pcis->vec1_D,0.0)); 872 } 873 /* Application of \widetilde{S}^-1 */ 874 PetscCall(PetscArrayzero(pcbddc->benign_p0,pcbddc->benign_n)); 875 PetscCall(PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,trans)); 876 PetscCall(PetscArrayzero(pcbddc->benign_p0,pcbddc->benign_n)); 877 PetscCall(VecSet(y,0.0)); 878 /* Application of B_delta */ 879 PetscCall(MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local)); 880 /* Contribution from boundary pressures */ 881 if (mat_ctx->C) { 882 const PetscScalar *lx; 883 PetscScalar *ly; 884 885 /* pressure ordered first in the local part of x and y */ 886 PetscCall(VecGetArrayRead(x,&lx)); 887 PetscCall(VecGetArray(y,&ly)); 888 PetscCall(VecPlaceArray(mat_ctx->xPg,lx)); 889 PetscCall(VecPlaceArray(mat_ctx->yPg,ly)); 890 if (trans) { 891 PetscCall(MatMultTranspose(mat_ctx->C,mat_ctx->xPg,mat_ctx->yPg)); 892 } else { 893 PetscCall(MatMult(mat_ctx->C,mat_ctx->xPg,mat_ctx->yPg)); 894 } 895 PetscCall(VecResetArray(mat_ctx->xPg)); 896 PetscCall(VecResetArray(mat_ctx->yPg)); 897 PetscCall(VecRestoreArrayRead(x,&lx)); 898 PetscCall(VecRestoreArray(y,&ly)); 899 } 900 /* Add contribution from saddle point */ 901 if (mat_ctx->l2g_p) { 902 if (trans) { 903 PetscCall(MatMultTranspose(mat_ctx->Bt_BB,pcis->vec1_B,mat_ctx->vP)); 904 } else { 905 PetscCall(MatMult(mat_ctx->B_BB,pcis->vec1_B,mat_ctx->vP)); 906 } 907 if (pcbddc->switch_static) { 908 if (trans) { 909 PetscCall(MatMultTransposeAdd(mat_ctx->Bt_BI,pcis->vec1_D,mat_ctx->vP,mat_ctx->vP)); 910 } else { 911 PetscCall(MatMultAdd(mat_ctx->B_BI,pcis->vec1_D,mat_ctx->vP,mat_ctx->vP)); 912 } 913 } 914 PetscCall(VecScatterBegin(mat_ctx->l2g_p,mat_ctx->vP,y,ADD_VALUES,SCATTER_FORWARD)); 915 PetscCall(VecScatterEnd(mat_ctx->l2g_p,mat_ctx->vP,y,ADD_VALUES,SCATTER_FORWARD)); 916 } 917 PetscCall(VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,y,ADD_VALUES,SCATTER_FORWARD)); 918 PetscCall(VecScatterEnd(mat_ctx->l2g_lambda,mat_ctx->lambda_local,y,ADD_VALUES,SCATTER_FORWARD)); 919 PetscFunctionReturn(0); 920 } 921 922 PetscErrorCode FETIDPMatMult(Mat fetimat, Vec x, Vec y) 923 { 924 PetscFunctionBegin; 925 PetscCall(FETIDPMatMult_Kernel(fetimat,x,y,PETSC_FALSE)); 926 PetscFunctionReturn(0); 927 } 928 929 PetscErrorCode FETIDPMatMultTranspose(Mat fetimat, Vec x, Vec y) 930 { 931 PetscFunctionBegin; 932 PetscCall(FETIDPMatMult_Kernel(fetimat,x,y,PETSC_TRUE)); 933 PetscFunctionReturn(0); 934 } 935 936 PetscErrorCode FETIDPPCApply_Kernel(PC fetipc, Vec x, Vec y, PetscBool trans) 937 { 938 FETIDPPC_ctx pc_ctx; 939 PC_IS *pcis; 940 941 PetscFunctionBegin; 942 PetscCall(PCShellGetContext(fetipc,&pc_ctx)); 943 pcis = (PC_IS*)pc_ctx->pc->data; 944 /* Application of B_Ddelta^T */ 945 PetscCall(VecScatterBegin(pc_ctx->l2g_lambda,x,pc_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE)); 946 PetscCall(VecScatterEnd(pc_ctx->l2g_lambda,x,pc_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE)); 947 PetscCall(VecSet(pcis->vec2_B,0.0)); 948 PetscCall(MatMultTranspose(pc_ctx->B_Ddelta,pc_ctx->lambda_local,pcis->vec2_B)); 949 /* Application of local Schur complement */ 950 if (trans) { 951 PetscCall(MatMultTranspose(pc_ctx->S_j,pcis->vec2_B,pcis->vec1_B)); 952 } else { 953 PetscCall(MatMult(pc_ctx->S_j,pcis->vec2_B,pcis->vec1_B)); 954 } 955 /* Application of B_Ddelta */ 956 PetscCall(MatMult(pc_ctx->B_Ddelta,pcis->vec1_B,pc_ctx->lambda_local)); 957 PetscCall(VecSet(y,0.0)); 958 PetscCall(VecScatterBegin(pc_ctx->l2g_lambda,pc_ctx->lambda_local,y,ADD_VALUES,SCATTER_FORWARD)); 959 PetscCall(VecScatterEnd(pc_ctx->l2g_lambda,pc_ctx->lambda_local,y,ADD_VALUES,SCATTER_FORWARD)); 960 PetscFunctionReturn(0); 961 } 962 963 PetscErrorCode FETIDPPCApply(PC pc, Vec x, Vec y) 964 { 965 PetscFunctionBegin; 966 PetscCall(FETIDPPCApply_Kernel(pc,x,y,PETSC_FALSE)); 967 PetscFunctionReturn(0); 968 } 969 970 PetscErrorCode FETIDPPCApplyTranspose(PC pc, Vec x, Vec y) 971 { 972 PetscFunctionBegin; 973 PetscCall(FETIDPPCApply_Kernel(pc,x,y,PETSC_TRUE)); 974 PetscFunctionReturn(0); 975 } 976 977 PetscErrorCode FETIDPPCView(PC pc, PetscViewer viewer) 978 { 979 FETIDPPC_ctx pc_ctx; 980 PetscBool iascii; 981 PetscViewer sviewer; 982 983 PetscFunctionBegin; 984 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 985 if (iascii) { 986 PetscMPIInt rank; 987 PetscBool isschur,isshell; 988 989 PetscCall(PCShellGetContext(pc,&pc_ctx)); 990 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank)); 991 PetscCall(PetscObjectTypeCompare((PetscObject)pc_ctx->S_j,MATSCHURCOMPLEMENT,&isschur)); 992 if (isschur) { 993 PetscCall(PetscViewerASCIIPrintf(viewer," Dirichlet preconditioner (just from rank 0)\n")); 994 } else { 995 PetscCall(PetscViewerASCIIPrintf(viewer," Lumped preconditioner (just from rank 0)\n")); 996 } 997 PetscCall(PetscViewerGetSubViewer(viewer,PetscObjectComm((PetscObject)pc_ctx->S_j),&sviewer)); 998 if (rank == 0) { 999 PetscCall(PetscViewerPushFormat(sviewer,PETSC_VIEWER_ASCII_INFO)); 1000 PetscCall(PetscViewerASCIIPushTab(sviewer)); 1001 PetscCall(MatView(pc_ctx->S_j,sviewer)); 1002 PetscCall(PetscViewerASCIIPopTab(sviewer)); 1003 PetscCall(PetscViewerPopFormat(sviewer)); 1004 } 1005 PetscCall(PetscViewerRestoreSubViewer(viewer,PetscObjectComm((PetscObject)pc_ctx->S_j),&sviewer)); 1006 PetscCall(PetscObjectTypeCompare((PetscObject)pc_ctx->B_Ddelta,MATSHELL,&isshell)); 1007 if (isshell) { 1008 BDdelta_DN ctx; 1009 PetscCall(PetscViewerASCIIPrintf(viewer," FETI-DP BDdelta: DB^t * (B D^-1 B^t)^-1 for deluxe scaling (just from rank 0)\n")); 1010 PetscCall(MatShellGetContext(pc_ctx->B_Ddelta,&ctx)); 1011 PetscCall(PetscViewerGetSubViewer(viewer,PetscObjectComm((PetscObject)pc_ctx->S_j),&sviewer)); 1012 if (rank == 0) { 1013 PetscInt tl; 1014 1015 PetscCall(PetscViewerASCIIGetTab(sviewer,&tl)); 1016 PetscCall(PetscObjectSetTabLevel((PetscObject)ctx->kBD,tl)); 1017 PetscCall(KSPView(ctx->kBD,sviewer)); 1018 PetscCall(PetscViewerPushFormat(sviewer,PETSC_VIEWER_ASCII_INFO)); 1019 PetscCall(MatView(ctx->BD,sviewer)); 1020 PetscCall(PetscViewerPopFormat(sviewer)); 1021 } 1022 PetscCall(PetscViewerRestoreSubViewer(viewer,PetscObjectComm((PetscObject)pc_ctx->S_j),&sviewer)); 1023 } 1024 PetscCall(PetscViewerFlush(viewer)); 1025 } 1026 PetscFunctionReturn(0); 1027 } 1028