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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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++) send_buffer[ptrs_buffer[i - 1] + j] = array[pcis->shared[i][j]]; 319 PetscCall(PetscMPIIntCast(ptrs_buffer[i] - ptrs_buffer[i - 1], &buf_size)); 320 PetscCall(PetscMPIIntCast(pcis->neigh[i], &neigh)); 321 PetscCallMPI(MPI_Isend(&send_buffer[ptrs_buffer[i - 1]], buf_size, MPIU_SCALAR, neigh, 0, comm, &send_reqs[i - 1])); 322 PetscCallMPI(MPI_Irecv(&recv_buffer[ptrs_buffer[i - 1]], buf_size, MPIU_SCALAR, neigh, 0, comm, &recv_reqs[i - 1])); 323 } 324 PetscCall(VecRestoreArrayRead(pcis->vec1_N, (const PetscScalar **)&array)); 325 if (pcis->n_neigh > 0) PetscCallMPI(MPI_Waitall(pcis->n_neigh - 1, recv_reqs, MPI_STATUSES_IGNORE)); 326 /* put values in correct places */ 327 for (i = 1; i < pcis->n_neigh; i++) { 328 for (j = 0; j < pcis->n_shared[i]; j++) { 329 k = pcis->shared[i][j]; 330 neigh_position = 0; 331 while (mat_graph->neighbours_set[k][neigh_position] != pcis->neigh[i]) neigh_position++; 332 all_factors[k][neigh_position] = recv_buffer[ptrs_buffer[i - 1] + j]; 333 } 334 } 335 if (pcis->n_neigh > 0) PetscCallMPI(MPI_Waitall(pcis->n_neigh - 1, send_reqs, MPI_STATUSES_IGNORE)); 336 PetscCall(PetscFree(send_reqs)); 337 PetscCall(PetscFree(recv_reqs)); 338 PetscCall(PetscFree(send_buffer)); 339 PetscCall(PetscFree(recv_buffer)); 340 PetscCall(PetscFree(ptrs_buffer)); 341 } 342 343 /* Compute B and B_delta (local actions) */ 344 PetscCall(PetscMalloc1(pcis->n_neigh, &aux_sums)); 345 PetscCall(PetscMalloc1(n_local_lambda, &l2g_indices)); 346 PetscCall(PetscMalloc1(n_local_lambda, &vals_B_delta)); 347 PetscCall(PetscMalloc1(n_local_lambda, &cols_B_delta)); 348 if (!pcbddc->use_deluxe_scaling) { 349 PetscCall(PetscMalloc1(n_local_lambda, &scaling_factors)); 350 } else { 351 scaling_factors = NULL; 352 all_factors = NULL; 353 } 354 PetscCall(ISGetIndices(subset_n, &aux_global_numbering)); 355 partial_sum = 0; 356 cum = 0; 357 for (i = 0; i < dual_size; i++) { 358 n_global_lambda = aux_global_numbering[cum]; 359 j = mat_graph->count[aux_local_numbering_1[i]]; 360 aux_sums[0] = 0; 361 for (s = 1; s < j; s++) aux_sums[s] = aux_sums[s - 1] + j - s + 1; 362 if (all_factors) array = all_factors[aux_local_numbering_1[i]]; 363 n_neg_values = 0; 364 while (n_neg_values < j && mat_graph->neighbours_set[aux_local_numbering_1[i]][n_neg_values] < rank) n_neg_values++; 365 n_pos_values = j - n_neg_values; 366 if (fully_redundant) { 367 for (s = 0; s < n_neg_values; s++) { 368 l2g_indices[partial_sum + s] = aux_sums[s] + n_neg_values - s - 1 + n_global_lambda; 369 cols_B_delta[partial_sum + s] = dual_dofs_boundary_indices[i]; 370 vals_B_delta[partial_sum + s] = -1.0; 371 if (!pcbddc->use_deluxe_scaling) scaling_factors[partial_sum + s] = array[s]; 372 } 373 for (s = 0; s < n_pos_values; s++) { 374 l2g_indices[partial_sum + s + n_neg_values] = aux_sums[n_neg_values] + s + n_global_lambda; 375 cols_B_delta[partial_sum + s + n_neg_values] = dual_dofs_boundary_indices[i]; 376 vals_B_delta[partial_sum + s + n_neg_values] = 1.0; 377 if (!pcbddc->use_deluxe_scaling) scaling_factors[partial_sum + s + n_neg_values] = array[s + n_neg_values]; 378 } 379 partial_sum += j; 380 } else { 381 /* l2g_indices and default cols and vals of B_delta */ 382 for (s = 0; s < j; s++) { 383 l2g_indices[partial_sum + s] = n_global_lambda + s; 384 cols_B_delta[partial_sum + s] = dual_dofs_boundary_indices[i]; 385 vals_B_delta[partial_sum + s] = 0.0; 386 } 387 /* B_delta */ 388 if (n_neg_values > 0) { /* there's a rank next to me to the left */ 389 vals_B_delta[partial_sum + n_neg_values - 1] = -1.0; 390 } 391 if (n_neg_values < j) { /* there's a rank next to me to the right */ 392 vals_B_delta[partial_sum + n_neg_values] = 1.0; 393 } 394 /* scaling as in Klawonn-Widlund 1999 */ 395 if (!pcbddc->use_deluxe_scaling) { 396 for (s = 0; s < n_neg_values; s++) { 397 scalar_value = 0.0; 398 for (k = 0; k < s + 1; k++) scalar_value += array[k]; 399 scaling_factors[partial_sum + s] = -scalar_value; 400 } 401 for (s = 0; s < n_pos_values; s++) { 402 scalar_value = 0.0; 403 for (k = s + n_neg_values; k < j; k++) scalar_value += array[k]; 404 scaling_factors[partial_sum + s + n_neg_values] = scalar_value; 405 } 406 } 407 partial_sum += j; 408 } 409 cum += aux_local_numbering_2[i]; 410 } 411 PetscCall(ISRestoreIndices(subset_n, &aux_global_numbering)); 412 PetscCall(ISDestroy(&subset_mult)); 413 PetscCall(ISDestroy(&subset_n)); 414 PetscCall(PetscFree(aux_sums)); 415 PetscCall(PetscFree(aux_local_numbering_1)); 416 PetscCall(PetscFree(dual_dofs_boundary_indices)); 417 if (all_factors) { 418 PetscCall(PetscFree(all_factors[0])); 419 PetscCall(PetscFree(all_factors)); 420 } 421 422 /* Create local part of B_delta */ 423 PetscCall(MatCreate(PETSC_COMM_SELF, &fetidpmat_ctx->B_delta)); 424 PetscCall(MatSetSizes(fetidpmat_ctx->B_delta, n_local_lambda, pcis->n_B, n_local_lambda, pcis->n_B)); 425 PetscCall(MatSetType(fetidpmat_ctx->B_delta, MATSEQAIJ)); 426 PetscCall(MatSeqAIJSetPreallocation(fetidpmat_ctx->B_delta, 1, NULL)); 427 PetscCall(MatSetOption(fetidpmat_ctx->B_delta, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE)); 428 for (i = 0; i < n_local_lambda; i++) PetscCall(MatSetValue(fetidpmat_ctx->B_delta, i, cols_B_delta[i], vals_B_delta[i], INSERT_VALUES)); 429 PetscCall(PetscFree(vals_B_delta)); 430 PetscCall(MatAssemblyBegin(fetidpmat_ctx->B_delta, MAT_FINAL_ASSEMBLY)); 431 PetscCall(MatAssemblyEnd(fetidpmat_ctx->B_delta, MAT_FINAL_ASSEMBLY)); 432 433 BD1 = NULL; 434 BD2 = NULL; 435 if (fully_redundant) { 436 PetscCheck(!pcbddc->use_deluxe_scaling, comm, PETSC_ERR_SUP, "Deluxe FETIDP with fully-redundant multipliers to be implemented"); 437 PetscCall(MatCreate(PETSC_COMM_SELF, &ScalingMat)); 438 PetscCall(MatSetSizes(ScalingMat, n_local_lambda, n_local_lambda, n_local_lambda, n_local_lambda)); 439 PetscCall(MatSetType(ScalingMat, MATSEQAIJ)); 440 PetscCall(MatSeqAIJSetPreallocation(ScalingMat, 1, NULL)); 441 for (i = 0; i < n_local_lambda; i++) PetscCall(MatSetValue(ScalingMat, i, i, scaling_factors[i], INSERT_VALUES)); 442 PetscCall(MatAssemblyBegin(ScalingMat, MAT_FINAL_ASSEMBLY)); 443 PetscCall(MatAssemblyEnd(ScalingMat, MAT_FINAL_ASSEMBLY)); 444 PetscCall(MatMatMult(ScalingMat, fetidpmat_ctx->B_delta, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &fetidpmat_ctx->B_Ddelta)); 445 PetscCall(MatDestroy(&ScalingMat)); 446 } else { 447 PetscCall(MatCreate(PETSC_COMM_SELF, &fetidpmat_ctx->B_Ddelta)); 448 PetscCall(MatSetSizes(fetidpmat_ctx->B_Ddelta, n_local_lambda, pcis->n_B, n_local_lambda, pcis->n_B)); 449 if (!pcbddc->use_deluxe_scaling || !pcbddc->sub_schurs) { 450 PetscCall(MatSetType(fetidpmat_ctx->B_Ddelta, MATSEQAIJ)); 451 PetscCall(MatSeqAIJSetPreallocation(fetidpmat_ctx->B_Ddelta, 1, NULL)); 452 for (i = 0; i < n_local_lambda; i++) PetscCall(MatSetValue(fetidpmat_ctx->B_Ddelta, i, cols_B_delta[i], scaling_factors[i], INSERT_VALUES)); 453 PetscCall(MatAssemblyBegin(fetidpmat_ctx->B_Ddelta, MAT_FINAL_ASSEMBLY)); 454 PetscCall(MatAssemblyEnd(fetidpmat_ctx->B_Ddelta, MAT_FINAL_ASSEMBLY)); 455 } else { 456 /* scaling as in Klawonn-Widlund 1999 */ 457 PCBDDCDeluxeScaling deluxe_ctx = pcbddc->deluxe_ctx; 458 PCBDDCSubSchurs sub_schurs = pcbddc->sub_schurs; 459 Mat T; 460 PetscScalar *W, lwork, *Bwork; 461 const PetscInt *idxs = NULL; 462 PetscInt cum, mss, *nnz; 463 PetscBLASInt *pivots, B_lwork, B_N, B_ierr; 464 465 PetscCheck(pcbddc->deluxe_singlemat, comm, PETSC_ERR_USER, "Cannot compute B_Ddelta! rerun with -pc_bddc_deluxe_singlemat"); 466 mss = 0; 467 PetscCall(PetscCalloc1(pcis->n_B, &nnz)); 468 if (sub_schurs->is_Ej_all) { 469 PetscCall(ISGetIndices(sub_schurs->is_Ej_all, &idxs)); 470 for (i = 0, cum = 0; i < sub_schurs->n_subs; i++) { 471 PetscInt subset_size; 472 473 PetscCall(ISGetLocalSize(sub_schurs->is_subs[i], &subset_size)); 474 for (j = 0; j < subset_size; j++) nnz[idxs[j + cum]] = subset_size; 475 mss = PetscMax(mss, subset_size); 476 cum += subset_size; 477 } 478 } 479 PetscCall(MatCreate(PETSC_COMM_SELF, &T)); 480 PetscCall(MatSetSizes(T, pcis->n_B, pcis->n_B, pcis->n_B, pcis->n_B)); 481 PetscCall(MatSetType(T, MATSEQAIJ)); 482 PetscCall(MatSeqAIJSetPreallocation(T, 0, nnz)); 483 PetscCall(PetscFree(nnz)); 484 485 /* workspace allocation */ 486 B_lwork = 0; 487 if (mss) { 488 PetscScalar dummy = 1; 489 490 B_lwork = -1; 491 PetscCall(PetscBLASIntCast(mss, &B_N)); 492 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 493 PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&B_N, &dummy, &B_N, &B_N, &lwork, &B_lwork, &B_ierr)); 494 PetscCall(PetscFPTrapPop()); 495 PetscCheck(!B_ierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in query to GETRI Lapack routine %d", (int)B_ierr); 496 PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(lwork), &B_lwork)); 497 } 498 PetscCall(PetscMalloc3(mss * mss, &W, mss, &pivots, B_lwork, &Bwork)); 499 500 for (i = 0, cum = 0; i < sub_schurs->n_subs; i++) { 501 const PetscScalar *M; 502 PetscInt subset_size; 503 504 PetscCall(ISGetLocalSize(sub_schurs->is_subs[i], &subset_size)); 505 PetscCall(PetscBLASIntCast(subset_size, &B_N)); 506 PetscCall(MatDenseGetArrayRead(deluxe_ctx->seq_mat[i], &M)); 507 PetscCall(PetscArraycpy(W, M, subset_size * subset_size)); 508 PetscCall(MatDenseRestoreArrayRead(deluxe_ctx->seq_mat[i], &M)); 509 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 510 PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&B_N, &B_N, W, &B_N, pivots, &B_ierr)); 511 PetscCheck(!B_ierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in GETRF Lapack routine %d", (int)B_ierr); 512 PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&B_N, W, &B_N, pivots, Bwork, &B_lwork, &B_ierr)); 513 PetscCheck(!B_ierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in GETRI Lapack routine %d", (int)B_ierr); 514 PetscCall(PetscFPTrapPop()); 515 /* silent static analyzer */ 516 PetscCheck(idxs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "IDXS not present"); 517 PetscCall(MatSetValues(T, subset_size, idxs + cum, subset_size, idxs + cum, W, INSERT_VALUES)); 518 cum += subset_size; 519 } 520 PetscCall(MatAssemblyBegin(T, MAT_FINAL_ASSEMBLY)); 521 PetscCall(MatAssemblyEnd(T, MAT_FINAL_ASSEMBLY)); 522 PetscCall(MatMatTransposeMult(T, fetidpmat_ctx->B_delta, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &BD1)); 523 PetscCall(MatMatMult(fetidpmat_ctx->B_delta, BD1, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &BD2)); 524 PetscCall(MatDestroy(&T)); 525 PetscCall(PetscFree3(W, pivots, Bwork)); 526 if (sub_schurs->is_Ej_all) PetscCall(ISRestoreIndices(sub_schurs->is_Ej_all, &idxs)); 527 } 528 } 529 PetscCall(PetscFree(scaling_factors)); 530 PetscCall(PetscFree(cols_B_delta)); 531 532 /* Layout of multipliers */ 533 PetscCall(PetscLayoutCreate(comm, &llay)); 534 PetscCall(PetscLayoutSetBlockSize(llay, 1)); 535 PetscCall(PetscLayoutSetSize(llay, fetidpmat_ctx->n_lambda)); 536 PetscCall(PetscLayoutSetUp(llay)); 537 PetscCall(PetscLayoutGetLocalSize(llay, &fetidpmat_ctx->n)); 538 539 /* Local work vector of multipliers */ 540 PetscCall(VecCreate(PETSC_COMM_SELF, &fetidpmat_ctx->lambda_local)); 541 PetscCall(VecSetSizes(fetidpmat_ctx->lambda_local, n_local_lambda, n_local_lambda)); 542 PetscCall(VecSetType(fetidpmat_ctx->lambda_local, VECSEQ)); 543 544 if (BD2) { 545 ISLocalToGlobalMapping l2g; 546 Mat T, TA, *pT; 547 IS is; 548 PetscInt nl, N; 549 BDdelta_DN ctx; 550 551 PetscCall(PetscLayoutGetLocalSize(llay, &nl)); 552 PetscCall(PetscLayoutGetSize(llay, &N)); 553 PetscCall(MatCreate(comm, &T)); 554 PetscCall(MatSetSizes(T, nl, nl, N, N)); 555 PetscCall(MatSetType(T, MATIS)); 556 PetscCall(ISLocalToGlobalMappingCreate(comm, 1, n_local_lambda, l2g_indices, PETSC_COPY_VALUES, &l2g)); 557 PetscCall(MatSetLocalToGlobalMapping(T, l2g, l2g)); 558 PetscCall(ISLocalToGlobalMappingDestroy(&l2g)); 559 PetscCall(MatISSetLocalMat(T, BD2)); 560 PetscCall(MatAssemblyBegin(T, MAT_FINAL_ASSEMBLY)); 561 PetscCall(MatAssemblyEnd(T, MAT_FINAL_ASSEMBLY)); 562 PetscCall(MatDestroy(&BD2)); 563 PetscCall(MatConvert(T, MATAIJ, MAT_INITIAL_MATRIX, &TA)); 564 PetscCall(MatDestroy(&T)); 565 PetscCall(ISCreateGeneral(comm, n_local_lambda, l2g_indices, PETSC_USE_POINTER, &is)); 566 PetscCall(MatCreateSubMatrices(TA, 1, &is, &is, MAT_INITIAL_MATRIX, &pT)); 567 PetscCall(MatDestroy(&TA)); 568 PetscCall(ISDestroy(&is)); 569 BD2 = pT[0]; 570 PetscCall(PetscFree(pT)); 571 572 /* B_Ddelta for non-redundant multipliers with deluxe scaling */ 573 PetscCall(PetscNew(&ctx)); 574 PetscCall(MatSetType(fetidpmat_ctx->B_Ddelta, MATSHELL)); 575 PetscCall(MatShellSetContext(fetidpmat_ctx->B_Ddelta, ctx)); 576 PetscCall(MatShellSetOperation(fetidpmat_ctx->B_Ddelta, MATOP_MULT, (void (*)(void))MatMult_BDdelta_deluxe_nonred)); 577 PetscCall(MatShellSetOperation(fetidpmat_ctx->B_Ddelta, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_BDdelta_deluxe_nonred)); 578 PetscCall(MatShellSetOperation(fetidpmat_ctx->B_Ddelta, MATOP_DESTROY, (void (*)(void))MatDestroy_BDdelta_deluxe_nonred)); 579 PetscCall(MatSetUp(fetidpmat_ctx->B_Ddelta)); 580 581 PetscCall(PetscObjectReference((PetscObject)BD1)); 582 ctx->BD = BD1; 583 PetscCall(KSPCreate(PETSC_COMM_SELF, &ctx->kBD)); 584 PetscCall(KSPSetOperators(ctx->kBD, BD2, BD2)); 585 PetscCall(VecDuplicate(fetidpmat_ctx->lambda_local, &ctx->work)); 586 fetidpmat_ctx->deluxe_nonred = PETSC_TRUE; 587 } 588 PetscCall(MatDestroy(&BD1)); 589 PetscCall(MatDestroy(&BD2)); 590 591 /* fetidpmat sizes */ 592 fetidpmat_ctx->n += nPgl; 593 fetidpmat_ctx->N = fetidpmat_ctx->n_lambda + nPg; 594 595 /* Global vector for FETI-DP linear system */ 596 PetscCall(VecCreate(comm, &fetidp_global)); 597 PetscCall(VecSetSizes(fetidp_global, fetidpmat_ctx->n, fetidpmat_ctx->N)); 598 PetscCall(VecSetType(fetidp_global, VECMPI)); 599 PetscCall(VecSetUp(fetidp_global)); 600 601 /* Decide layout for fetidp dofs: if it is a saddle point problem 602 pressure is ordered first in the local part of the global vector 603 of the FETI-DP linear system */ 604 if (nPg) { 605 Vec v; 606 IS IS_l2g_p, ais; 607 PetscLayout alay; 608 const PetscInt *idxs, *pranges, *aranges, *lranges; 609 PetscInt *l2g_indices_p, rst; 610 PetscMPIInt rank; 611 612 PetscCall(PetscMalloc1(nPl, &l2g_indices_p)); 613 PetscCall(VecGetLayout(fetidp_global, &alay)); 614 PetscCall(PetscLayoutGetRanges(alay, &aranges)); 615 PetscCall(PetscLayoutGetRanges(play, &pranges)); 616 PetscCall(PetscLayoutGetRanges(llay, &lranges)); 617 618 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)fetidp_global), &rank)); 619 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)fetidp_global), pranges[rank + 1] - pranges[rank], aranges[rank], 1, &fetidpmat_ctx->pressure)); 620 PetscCall(PetscObjectSetName((PetscObject)fetidpmat_ctx->pressure, "F_P")); 621 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)fetidp_global), lranges[rank + 1] - lranges[rank], aranges[rank] + pranges[rank + 1] - pranges[rank], 1, &fetidpmat_ctx->lagrange)); 622 PetscCall(PetscObjectSetName((PetscObject)fetidpmat_ctx->lagrange, "F_L")); 623 PetscCall(ISLocalToGlobalMappingGetIndices(l2gmap_p, &idxs)); 624 /* shift local to global indices for pressure */ 625 for (i = 0; i < nPl; i++) { 626 PetscMPIInt owner; 627 628 PetscCall(PetscLayoutFindOwner(play, idxs[i], &owner)); 629 l2g_indices_p[i] = idxs[i] - pranges[owner] + aranges[owner]; 630 } 631 PetscCall(ISLocalToGlobalMappingRestoreIndices(l2gmap_p, &idxs)); 632 PetscCall(ISCreateGeneral(comm, nPl, l2g_indices_p, PETSC_OWN_POINTER, &IS_l2g_p)); 633 634 /* local to global scatter for pressure */ 635 PetscCall(VecScatterCreate(fetidpmat_ctx->vP, NULL, fetidp_global, IS_l2g_p, &fetidpmat_ctx->l2g_p)); 636 PetscCall(ISDestroy(&IS_l2g_p)); 637 638 /* scatter for lagrange multipliers only */ 639 PetscCall(VecCreate(comm, &v)); 640 PetscCall(VecSetType(v, VECSTANDARD)); 641 PetscCall(VecSetLayout(v, llay)); 642 PetscCall(VecSetUp(v)); 643 PetscCall(ISCreateGeneral(comm, n_local_lambda, l2g_indices, PETSC_COPY_VALUES, &ais)); 644 PetscCall(VecScatterCreate(fetidpmat_ctx->lambda_local, NULL, v, ais, &fetidpmat_ctx->l2g_lambda_only)); 645 PetscCall(ISDestroy(&ais)); 646 PetscCall(VecDestroy(&v)); 647 648 /* shift local to global indices for multipliers */ 649 for (i = 0; i < n_local_lambda; i++) { 650 PetscInt ps; 651 PetscMPIInt owner; 652 653 PetscCall(PetscLayoutFindOwner(llay, l2g_indices[i], &owner)); 654 ps = pranges[owner + 1] - pranges[owner]; 655 l2g_indices[i] = l2g_indices[i] - lranges[owner] + aranges[owner] + ps; 656 } 657 658 /* scatter from alldofs to pressures global fetidp vector */ 659 PetscCall(PetscLayoutGetRange(alay, &rst, NULL)); 660 PetscCall(ISCreateStride(comm, nPgl, rst, 1, &ais)); 661 PetscCall(VecScatterCreate(pcis->vec1_global, pP, fetidp_global, ais, &fetidpmat_ctx->g2g_p)); 662 PetscCall(ISDestroy(&ais)); 663 } 664 PetscCall(PetscLayoutDestroy(&llay)); 665 PetscCall(PetscLayoutDestroy(&play)); 666 PetscCall(ISCreateGeneral(comm, n_local_lambda, l2g_indices, PETSC_OWN_POINTER, &IS_l2g_lambda)); 667 668 /* scatter from local to global multipliers */ 669 PetscCall(VecScatterCreate(fetidpmat_ctx->lambda_local, NULL, fetidp_global, IS_l2g_lambda, &fetidpmat_ctx->l2g_lambda)); 670 PetscCall(ISDestroy(&IS_l2g_lambda)); 671 PetscCall(ISLocalToGlobalMappingDestroy(&l2gmap_p)); 672 PetscCall(VecDestroy(&fetidp_global)); 673 674 /* Create some work vectors needed by fetidp */ 675 PetscCall(VecDuplicate(pcis->vec1_B, &fetidpmat_ctx->temp_solution_B)); 676 PetscCall(VecDuplicate(pcis->vec1_D, &fetidpmat_ctx->temp_solution_D)); 677 PetscFunctionReturn(PETSC_SUCCESS); 678 } 679 680 PetscErrorCode PCBDDCSetupFETIDPPCContext(Mat fetimat, FETIDPPC_ctx fetidppc_ctx) 681 { 682 FETIDPMat_ctx mat_ctx; 683 PC_BDDC *pcbddc = (PC_BDDC *)fetidppc_ctx->pc->data; 684 PC_IS *pcis = (PC_IS *)fetidppc_ctx->pc->data; 685 PetscBool lumped = PETSC_FALSE; 686 687 PetscFunctionBegin; 688 PetscCall(MatShellGetContext(fetimat, &mat_ctx)); 689 /* get references from objects created when setting up feti mat context */ 690 PetscCall(PetscObjectReference((PetscObject)mat_ctx->lambda_local)); 691 fetidppc_ctx->lambda_local = mat_ctx->lambda_local; 692 PetscCall(PetscObjectReference((PetscObject)mat_ctx->B_Ddelta)); 693 fetidppc_ctx->B_Ddelta = mat_ctx->B_Ddelta; 694 if (mat_ctx->deluxe_nonred) { 695 PC pc, mpc; 696 BDdelta_DN ctx; 697 MatSolverType solver; 698 const char *prefix; 699 700 PetscCall(MatShellGetContext(mat_ctx->B_Ddelta, &ctx)); 701 PetscCall(KSPSetType(ctx->kBD, KSPPREONLY)); 702 PetscCall(KSPGetPC(ctx->kBD, &mpc)); 703 PetscCall(KSPGetPC(pcbddc->ksp_D, &pc)); 704 PetscCall(PCSetType(mpc, PCLU)); 705 PetscCall(PCFactorGetMatSolverType(pc, (MatSolverType *)&solver)); 706 if (solver) PetscCall(PCFactorSetMatSolverType(mpc, solver)); 707 PetscCall(MatGetOptionsPrefix(fetimat, &prefix)); 708 PetscCall(KSPSetOptionsPrefix(ctx->kBD, prefix)); 709 PetscCall(KSPAppendOptionsPrefix(ctx->kBD, "bddelta_")); 710 PetscCall(KSPSetFromOptions(ctx->kBD)); 711 } 712 713 if (mat_ctx->l2g_lambda_only) { 714 PetscCall(PetscObjectReference((PetscObject)mat_ctx->l2g_lambda_only)); 715 fetidppc_ctx->l2g_lambda = mat_ctx->l2g_lambda_only; 716 } else { 717 PetscCall(PetscObjectReference((PetscObject)mat_ctx->l2g_lambda)); 718 fetidppc_ctx->l2g_lambda = mat_ctx->l2g_lambda; 719 } 720 /* Dirichlet preconditioner */ 721 PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)fetimat)->prefix, "-pc_lumped", &lumped, NULL)); 722 if (!lumped) { 723 IS iV; 724 PetscBool discrete_harmonic = PETSC_FALSE; 725 726 PetscCall(PetscObjectQuery((PetscObject)fetidppc_ctx->pc, "__KSPFETIDP_iV", (PetscObject *)&iV)); 727 if (iV) PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)fetimat)->prefix, "-pc_discrete_harmonic", &discrete_harmonic, NULL)); 728 if (discrete_harmonic) { 729 KSP sksp; 730 PC pc; 731 PCBDDCSubSchurs sub_schurs = pcbddc->sub_schurs; 732 Mat A_II, A_IB, A_BI; 733 IS iP = NULL; 734 PetscBool isshell, reuse = PETSC_FALSE; 735 KSPType ksptype; 736 const char *prefix; 737 738 /* 739 We constructs a Schur complement for 740 741 | A_II A_ID | 742 | A_DI A_DD | 743 744 instead of 745 746 | A_II B^t_II A_ID | 747 | B_II -C_II B_ID | 748 | A_DI B^t_ID A_DD | 749 750 */ 751 if (sub_schurs && sub_schurs->reuse_solver) { 752 PetscCall(PetscObjectQuery((PetscObject)sub_schurs->A, "__KSPFETIDP_iP", (PetscObject *)&iP)); 753 if (iP) reuse = PETSC_TRUE; 754 } 755 if (!reuse) { 756 IS aB; 757 PetscInt nb; 758 PetscCall(ISGetLocalSize(pcis->is_B_local, &nb)); 759 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pcis->A_II), nb, 0, 1, &aB)); 760 PetscCall(MatCreateSubMatrix(pcis->A_II, iV, iV, MAT_INITIAL_MATRIX, &A_II)); 761 PetscCall(MatCreateSubMatrix(pcis->A_IB, iV, aB, MAT_INITIAL_MATRIX, &A_IB)); 762 PetscCall(MatCreateSubMatrix(pcis->A_BI, aB, iV, MAT_INITIAL_MATRIX, &A_BI)); 763 PetscCall(ISDestroy(&aB)); 764 } else { 765 PetscCall(MatCreateSubMatrix(sub_schurs->A, pcis->is_I_local, pcis->is_B_local, MAT_INITIAL_MATRIX, &A_IB)); 766 PetscCall(MatCreateSubMatrix(sub_schurs->A, pcis->is_B_local, pcis->is_I_local, MAT_INITIAL_MATRIX, &A_BI)); 767 PetscCall(PetscObjectReference((PetscObject)pcis->A_II)); 768 A_II = pcis->A_II; 769 } 770 PetscCall(MatCreateSchurComplement(A_II, A_II, A_IB, A_BI, pcis->A_BB, &fetidppc_ctx->S_j)); 771 772 /* propagate settings of solver */ 773 PetscCall(MatSchurComplementGetKSP(fetidppc_ctx->S_j, &sksp)); 774 PetscCall(KSPGetType(pcis->ksp_D, &ksptype)); 775 PetscCall(KSPSetType(sksp, ksptype)); 776 PetscCall(KSPGetPC(pcis->ksp_D, &pc)); 777 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCSHELL, &isshell)); 778 if (!isshell) { 779 MatSolverType solver; 780 PCType pctype; 781 782 PetscCall(PCGetType(pc, &pctype)); 783 PetscCall(PCFactorGetMatSolverType(pc, (MatSolverType *)&solver)); 784 PetscCall(KSPGetPC(sksp, &pc)); 785 PetscCall(PCSetType(pc, pctype)); 786 if (solver) PetscCall(PCFactorSetMatSolverType(pc, solver)); 787 } else { 788 PetscCall(KSPGetPC(sksp, &pc)); 789 PetscCall(PCSetType(pc, PCLU)); 790 } 791 PetscCall(MatDestroy(&A_II)); 792 PetscCall(MatDestroy(&A_IB)); 793 PetscCall(MatDestroy(&A_BI)); 794 PetscCall(MatGetOptionsPrefix(fetimat, &prefix)); 795 PetscCall(KSPSetOptionsPrefix(sksp, prefix)); 796 PetscCall(KSPAppendOptionsPrefix(sksp, "harmonic_")); 797 PetscCall(KSPSetFromOptions(sksp)); 798 if (reuse) { 799 PetscCall(KSPSetPC(sksp, sub_schurs->reuse_solver->interior_solver)); 800 PetscCall(PetscObjectIncrementTabLevel((PetscObject)sub_schurs->reuse_solver->interior_solver, (PetscObject)sksp, 0)); 801 } 802 } else { /* default Dirichlet preconditioner is pde-harmonic */ 803 PetscCall(MatCreateSchurComplement(pcis->A_II, pcis->A_II, pcis->A_IB, pcis->A_BI, pcis->A_BB, &fetidppc_ctx->S_j)); 804 PetscCall(MatSchurComplementSetKSP(fetidppc_ctx->S_j, pcis->ksp_D)); 805 } 806 } else { 807 PetscCall(PetscObjectReference((PetscObject)pcis->A_BB)); 808 fetidppc_ctx->S_j = pcis->A_BB; 809 } 810 /* saddle-point */ 811 if (mat_ctx->xPg) { 812 PetscCall(PetscObjectReference((PetscObject)mat_ctx->xPg)); 813 fetidppc_ctx->xPg = mat_ctx->xPg; 814 PetscCall(PetscObjectReference((PetscObject)mat_ctx->yPg)); 815 fetidppc_ctx->yPg = mat_ctx->yPg; 816 } 817 PetscFunctionReturn(PETSC_SUCCESS); 818 } 819 820 PetscErrorCode FETIDPMatMult_Kernel(Mat fetimat, Vec x, Vec y, PetscBool trans) 821 { 822 FETIDPMat_ctx mat_ctx; 823 PC_BDDC *pcbddc; 824 PC_IS *pcis; 825 826 PetscFunctionBegin; 827 PetscCall(MatShellGetContext(fetimat, &mat_ctx)); 828 pcis = (PC_IS *)mat_ctx->pc->data; 829 pcbddc = (PC_BDDC *)mat_ctx->pc->data; 830 /* Application of B_delta^T */ 831 PetscCall(VecSet(pcis->vec1_B, 0.)); 832 PetscCall(VecScatterBegin(mat_ctx->l2g_lambda, x, mat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE)); 833 PetscCall(VecScatterEnd(mat_ctx->l2g_lambda, x, mat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE)); 834 PetscCall(MatMultTranspose(mat_ctx->B_delta, mat_ctx->lambda_local, pcis->vec1_B)); 835 836 /* Add contribution from saddle point */ 837 if (mat_ctx->l2g_p) { 838 PetscCall(VecScatterBegin(mat_ctx->l2g_p, x, mat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE)); 839 PetscCall(VecScatterEnd(mat_ctx->l2g_p, x, mat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE)); 840 if (pcbddc->switch_static) { 841 if (trans) { 842 PetscCall(MatMultTranspose(mat_ctx->B_BI, mat_ctx->vP, pcis->vec1_D)); 843 } else { 844 PetscCall(MatMult(mat_ctx->Bt_BI, mat_ctx->vP, pcis->vec1_D)); 845 } 846 } 847 if (trans) { 848 PetscCall(MatMultTransposeAdd(mat_ctx->B_BB, mat_ctx->vP, pcis->vec1_B, pcis->vec1_B)); 849 } else { 850 PetscCall(MatMultAdd(mat_ctx->Bt_BB, mat_ctx->vP, pcis->vec1_B, pcis->vec1_B)); 851 } 852 } else { 853 if (pcbddc->switch_static) PetscCall(VecSet(pcis->vec1_D, 0.0)); 854 } 855 /* Application of \widetilde{S}^-1 */ 856 PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n)); 857 PetscCall(PCBDDCApplyInterfacePreconditioner(mat_ctx->pc, trans)); 858 PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n)); 859 PetscCall(VecSet(y, 0.0)); 860 /* Application of B_delta */ 861 PetscCall(MatMult(mat_ctx->B_delta, pcis->vec1_B, mat_ctx->lambda_local)); 862 /* Contribution from boundary pressures */ 863 if (mat_ctx->C) { 864 const PetscScalar *lx; 865 PetscScalar *ly; 866 867 /* pressure ordered first in the local part of x and y */ 868 PetscCall(VecGetArrayRead(x, &lx)); 869 PetscCall(VecGetArray(y, &ly)); 870 PetscCall(VecPlaceArray(mat_ctx->xPg, lx)); 871 PetscCall(VecPlaceArray(mat_ctx->yPg, ly)); 872 if (trans) { 873 PetscCall(MatMultTranspose(mat_ctx->C, mat_ctx->xPg, mat_ctx->yPg)); 874 } else { 875 PetscCall(MatMult(mat_ctx->C, mat_ctx->xPg, mat_ctx->yPg)); 876 } 877 PetscCall(VecResetArray(mat_ctx->xPg)); 878 PetscCall(VecResetArray(mat_ctx->yPg)); 879 PetscCall(VecRestoreArrayRead(x, &lx)); 880 PetscCall(VecRestoreArray(y, &ly)); 881 } 882 /* Add contribution from saddle point */ 883 if (mat_ctx->l2g_p) { 884 if (trans) { 885 PetscCall(MatMultTranspose(mat_ctx->Bt_BB, pcis->vec1_B, mat_ctx->vP)); 886 } else { 887 PetscCall(MatMult(mat_ctx->B_BB, pcis->vec1_B, mat_ctx->vP)); 888 } 889 if (pcbddc->switch_static) { 890 if (trans) { 891 PetscCall(MatMultTransposeAdd(mat_ctx->Bt_BI, pcis->vec1_D, mat_ctx->vP, mat_ctx->vP)); 892 } else { 893 PetscCall(MatMultAdd(mat_ctx->B_BI, pcis->vec1_D, mat_ctx->vP, mat_ctx->vP)); 894 } 895 } 896 PetscCall(VecScatterBegin(mat_ctx->l2g_p, mat_ctx->vP, y, ADD_VALUES, SCATTER_FORWARD)); 897 PetscCall(VecScatterEnd(mat_ctx->l2g_p, mat_ctx->vP, y, ADD_VALUES, SCATTER_FORWARD)); 898 } 899 PetscCall(VecScatterBegin(mat_ctx->l2g_lambda, mat_ctx->lambda_local, y, ADD_VALUES, SCATTER_FORWARD)); 900 PetscCall(VecScatterEnd(mat_ctx->l2g_lambda, mat_ctx->lambda_local, y, ADD_VALUES, SCATTER_FORWARD)); 901 PetscFunctionReturn(PETSC_SUCCESS); 902 } 903 904 PetscErrorCode FETIDPMatMult(Mat fetimat, Vec x, Vec y) 905 { 906 PetscFunctionBegin; 907 PetscCall(FETIDPMatMult_Kernel(fetimat, x, y, PETSC_FALSE)); 908 PetscFunctionReturn(PETSC_SUCCESS); 909 } 910 911 PetscErrorCode FETIDPMatMultTranspose(Mat fetimat, Vec x, Vec y) 912 { 913 PetscFunctionBegin; 914 PetscCall(FETIDPMatMult_Kernel(fetimat, x, y, PETSC_TRUE)); 915 PetscFunctionReturn(PETSC_SUCCESS); 916 } 917 918 PetscErrorCode FETIDPPCApply_Kernel(PC fetipc, Vec x, Vec y, PetscBool trans) 919 { 920 FETIDPPC_ctx pc_ctx; 921 PC_IS *pcis; 922 923 PetscFunctionBegin; 924 PetscCall(PCShellGetContext(fetipc, &pc_ctx)); 925 pcis = (PC_IS *)pc_ctx->pc->data; 926 /* Application of B_Ddelta^T */ 927 PetscCall(VecScatterBegin(pc_ctx->l2g_lambda, x, pc_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE)); 928 PetscCall(VecScatterEnd(pc_ctx->l2g_lambda, x, pc_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE)); 929 PetscCall(VecSet(pcis->vec2_B, 0.0)); 930 PetscCall(MatMultTranspose(pc_ctx->B_Ddelta, pc_ctx->lambda_local, pcis->vec2_B)); 931 /* Application of local Schur complement */ 932 if (trans) { 933 PetscCall(MatMultTranspose(pc_ctx->S_j, pcis->vec2_B, pcis->vec1_B)); 934 } else { 935 PetscCall(MatMult(pc_ctx->S_j, pcis->vec2_B, pcis->vec1_B)); 936 } 937 /* Application of B_Ddelta */ 938 PetscCall(MatMult(pc_ctx->B_Ddelta, pcis->vec1_B, pc_ctx->lambda_local)); 939 PetscCall(VecSet(y, 0.0)); 940 PetscCall(VecScatterBegin(pc_ctx->l2g_lambda, pc_ctx->lambda_local, y, ADD_VALUES, SCATTER_FORWARD)); 941 PetscCall(VecScatterEnd(pc_ctx->l2g_lambda, pc_ctx->lambda_local, y, ADD_VALUES, SCATTER_FORWARD)); 942 PetscFunctionReturn(PETSC_SUCCESS); 943 } 944 945 PetscErrorCode FETIDPPCApply(PC pc, Vec x, Vec y) 946 { 947 PetscFunctionBegin; 948 PetscCall(FETIDPPCApply_Kernel(pc, x, y, PETSC_FALSE)); 949 PetscFunctionReturn(PETSC_SUCCESS); 950 } 951 952 PetscErrorCode FETIDPPCApplyTranspose(PC pc, Vec x, Vec y) 953 { 954 PetscFunctionBegin; 955 PetscCall(FETIDPPCApply_Kernel(pc, x, y, PETSC_TRUE)); 956 PetscFunctionReturn(PETSC_SUCCESS); 957 } 958 959 PetscErrorCode FETIDPPCView(PC pc, PetscViewer viewer) 960 { 961 FETIDPPC_ctx pc_ctx; 962 PetscBool iascii; 963 PetscViewer sviewer; 964 965 PetscFunctionBegin; 966 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 967 if (iascii) { 968 PetscMPIInt rank; 969 PetscBool isschur, isshell; 970 971 PetscCall(PCShellGetContext(pc, &pc_ctx)); 972 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank)); 973 PetscCall(PetscObjectTypeCompare((PetscObject)pc_ctx->S_j, MATSCHURCOMPLEMENT, &isschur)); 974 if (isschur) { 975 PetscCall(PetscViewerASCIIPrintf(viewer, " Dirichlet preconditioner (just from rank 0)\n")); 976 } else { 977 PetscCall(PetscViewerASCIIPrintf(viewer, " Lumped preconditioner (just from rank 0)\n")); 978 } 979 PetscCall(PetscViewerGetSubViewer(viewer, PetscObjectComm((PetscObject)pc_ctx->S_j), &sviewer)); 980 if (rank == 0) { 981 PetscCall(PetscViewerPushFormat(sviewer, PETSC_VIEWER_ASCII_INFO)); 982 PetscCall(PetscViewerASCIIPushTab(sviewer)); 983 PetscCall(MatView(pc_ctx->S_j, sviewer)); 984 PetscCall(PetscViewerASCIIPopTab(sviewer)); 985 PetscCall(PetscViewerPopFormat(sviewer)); 986 } 987 PetscCall(PetscViewerRestoreSubViewer(viewer, PetscObjectComm((PetscObject)pc_ctx->S_j), &sviewer)); 988 PetscCall(PetscObjectTypeCompare((PetscObject)pc_ctx->B_Ddelta, MATSHELL, &isshell)); 989 if (isshell) { 990 BDdelta_DN ctx; 991 PetscCall(PetscViewerASCIIPrintf(viewer, " FETI-DP BDdelta: DB^t * (B D^-1 B^t)^-1 for deluxe scaling (just from rank 0)\n")); 992 PetscCall(MatShellGetContext(pc_ctx->B_Ddelta, &ctx)); 993 PetscCall(PetscViewerGetSubViewer(viewer, PetscObjectComm((PetscObject)pc_ctx->S_j), &sviewer)); 994 if (rank == 0) { 995 PetscInt tl; 996 997 PetscCall(PetscViewerASCIIGetTab(sviewer, &tl)); 998 PetscCall(PetscObjectSetTabLevel((PetscObject)ctx->kBD, tl)); 999 PetscCall(KSPView(ctx->kBD, sviewer)); 1000 PetscCall(PetscViewerPushFormat(sviewer, PETSC_VIEWER_ASCII_INFO)); 1001 PetscCall(MatView(ctx->BD, sviewer)); 1002 PetscCall(PetscViewerPopFormat(sviewer)); 1003 } 1004 PetscCall(PetscViewerRestoreSubViewer(viewer, PetscObjectComm((PetscObject)pc_ctx->S_j), &sviewer)); 1005 } 1006 PetscCall(PetscViewerFlush(viewer)); 1007 } 1008 PetscFunctionReturn(PETSC_SUCCESS); 1009 } 1010