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