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(KSPSetNestLevel(ctx->kBD, fetidpmat_ctx->pc->kspnestlevel)); 585 PetscCall(KSPSetOperators(ctx->kBD, BD2, BD2)); 586 PetscCall(VecDuplicate(fetidpmat_ctx->lambda_local, &ctx->work)); 587 fetidpmat_ctx->deluxe_nonred = PETSC_TRUE; 588 } 589 PetscCall(MatDestroy(&BD1)); 590 PetscCall(MatDestroy(&BD2)); 591 592 /* fetidpmat sizes */ 593 fetidpmat_ctx->n += nPgl; 594 fetidpmat_ctx->N = fetidpmat_ctx->n_lambda + nPg; 595 596 /* Global vector for FETI-DP linear system */ 597 PetscCall(VecCreate(comm, &fetidp_global)); 598 PetscCall(VecSetSizes(fetidp_global, fetidpmat_ctx->n, fetidpmat_ctx->N)); 599 PetscCall(VecSetType(fetidp_global, VECMPI)); 600 PetscCall(VecSetUp(fetidp_global)); 601 602 /* Decide layout for fetidp dofs: if it is a saddle point problem 603 pressure is ordered first in the local part of the global vector 604 of the FETI-DP linear system */ 605 if (nPg) { 606 Vec v; 607 IS IS_l2g_p, ais; 608 PetscLayout alay; 609 const PetscInt *idxs, *pranges, *aranges, *lranges; 610 PetscInt *l2g_indices_p, rst; 611 PetscMPIInt rank; 612 613 PetscCall(PetscMalloc1(nPl, &l2g_indices_p)); 614 PetscCall(VecGetLayout(fetidp_global, &alay)); 615 PetscCall(PetscLayoutGetRanges(alay, &aranges)); 616 PetscCall(PetscLayoutGetRanges(play, &pranges)); 617 PetscCall(PetscLayoutGetRanges(llay, &lranges)); 618 619 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)fetidp_global), &rank)); 620 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)fetidp_global), pranges[rank + 1] - pranges[rank], aranges[rank], 1, &fetidpmat_ctx->pressure)); 621 PetscCall(PetscObjectSetName((PetscObject)fetidpmat_ctx->pressure, "F_P")); 622 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)fetidp_global), lranges[rank + 1] - lranges[rank], aranges[rank] + pranges[rank + 1] - pranges[rank], 1, &fetidpmat_ctx->lagrange)); 623 PetscCall(PetscObjectSetName((PetscObject)fetidpmat_ctx->lagrange, "F_L")); 624 PetscCall(ISLocalToGlobalMappingGetIndices(l2gmap_p, &idxs)); 625 /* shift local to global indices for pressure */ 626 for (i = 0; i < nPl; i++) { 627 PetscMPIInt owner; 628 629 PetscCall(PetscLayoutFindOwner(play, idxs[i], &owner)); 630 l2g_indices_p[i] = idxs[i] - pranges[owner] + aranges[owner]; 631 } 632 PetscCall(ISLocalToGlobalMappingRestoreIndices(l2gmap_p, &idxs)); 633 PetscCall(ISCreateGeneral(comm, nPl, l2g_indices_p, PETSC_OWN_POINTER, &IS_l2g_p)); 634 635 /* local to global scatter for pressure */ 636 PetscCall(VecScatterCreate(fetidpmat_ctx->vP, NULL, fetidp_global, IS_l2g_p, &fetidpmat_ctx->l2g_p)); 637 PetscCall(ISDestroy(&IS_l2g_p)); 638 639 /* scatter for lagrange multipliers only */ 640 PetscCall(VecCreate(comm, &v)); 641 PetscCall(VecSetType(v, VECSTANDARD)); 642 PetscCall(VecSetLayout(v, llay)); 643 PetscCall(VecSetUp(v)); 644 PetscCall(ISCreateGeneral(comm, n_local_lambda, l2g_indices, PETSC_COPY_VALUES, &ais)); 645 PetscCall(VecScatterCreate(fetidpmat_ctx->lambda_local, NULL, v, ais, &fetidpmat_ctx->l2g_lambda_only)); 646 PetscCall(ISDestroy(&ais)); 647 PetscCall(VecDestroy(&v)); 648 649 /* shift local to global indices for multipliers */ 650 for (i = 0; i < n_local_lambda; i++) { 651 PetscInt ps; 652 PetscMPIInt owner; 653 654 PetscCall(PetscLayoutFindOwner(llay, l2g_indices[i], &owner)); 655 ps = pranges[owner + 1] - pranges[owner]; 656 l2g_indices[i] = l2g_indices[i] - lranges[owner] + aranges[owner] + ps; 657 } 658 659 /* scatter from alldofs to pressures global fetidp vector */ 660 PetscCall(PetscLayoutGetRange(alay, &rst, NULL)); 661 PetscCall(ISCreateStride(comm, nPgl, rst, 1, &ais)); 662 PetscCall(VecScatterCreate(pcis->vec1_global, pP, fetidp_global, ais, &fetidpmat_ctx->g2g_p)); 663 PetscCall(ISDestroy(&ais)); 664 } 665 PetscCall(PetscLayoutDestroy(&llay)); 666 PetscCall(PetscLayoutDestroy(&play)); 667 PetscCall(ISCreateGeneral(comm, n_local_lambda, l2g_indices, PETSC_OWN_POINTER, &IS_l2g_lambda)); 668 669 /* scatter from local to global multipliers */ 670 PetscCall(VecScatterCreate(fetidpmat_ctx->lambda_local, NULL, fetidp_global, IS_l2g_lambda, &fetidpmat_ctx->l2g_lambda)); 671 PetscCall(ISDestroy(&IS_l2g_lambda)); 672 PetscCall(ISLocalToGlobalMappingDestroy(&l2gmap_p)); 673 PetscCall(VecDestroy(&fetidp_global)); 674 675 /* Create some work vectors needed by fetidp */ 676 PetscCall(VecDuplicate(pcis->vec1_B, &fetidpmat_ctx->temp_solution_B)); 677 PetscCall(VecDuplicate(pcis->vec1_D, &fetidpmat_ctx->temp_solution_D)); 678 PetscFunctionReturn(PETSC_SUCCESS); 679 } 680 681 PetscErrorCode PCBDDCSetupFETIDPPCContext(Mat fetimat, FETIDPPC_ctx fetidppc_ctx) 682 { 683 FETIDPMat_ctx mat_ctx; 684 PC_BDDC *pcbddc = (PC_BDDC *)fetidppc_ctx->pc->data; 685 PC_IS *pcis = (PC_IS *)fetidppc_ctx->pc->data; 686 PetscBool lumped = PETSC_FALSE; 687 688 PetscFunctionBegin; 689 PetscCall(MatShellGetContext(fetimat, &mat_ctx)); 690 /* get references from objects created when setting up feti mat context */ 691 PetscCall(PetscObjectReference((PetscObject)mat_ctx->lambda_local)); 692 fetidppc_ctx->lambda_local = mat_ctx->lambda_local; 693 PetscCall(PetscObjectReference((PetscObject)mat_ctx->B_Ddelta)); 694 fetidppc_ctx->B_Ddelta = mat_ctx->B_Ddelta; 695 if (mat_ctx->deluxe_nonred) { 696 PC pc, mpc; 697 BDdelta_DN ctx; 698 MatSolverType solver; 699 const char *prefix; 700 701 PetscCall(MatShellGetContext(mat_ctx->B_Ddelta, &ctx)); 702 PetscCall(KSPSetType(ctx->kBD, KSPPREONLY)); 703 PetscCall(KSPGetPC(ctx->kBD, &mpc)); 704 PetscCall(KSPGetPC(pcbddc->ksp_D, &pc)); 705 PetscCall(PCSetType(mpc, PCLU)); 706 PetscCall(PCFactorGetMatSolverType(pc, (MatSolverType *)&solver)); 707 if (solver) PetscCall(PCFactorSetMatSolverType(mpc, solver)); 708 PetscCall(MatGetOptionsPrefix(fetimat, &prefix)); 709 PetscCall(KSPSetOptionsPrefix(ctx->kBD, prefix)); 710 PetscCall(KSPAppendOptionsPrefix(ctx->kBD, "bddelta_")); 711 PetscCall(KSPSetFromOptions(ctx->kBD)); 712 } 713 714 if (mat_ctx->l2g_lambda_only) { 715 PetscCall(PetscObjectReference((PetscObject)mat_ctx->l2g_lambda_only)); 716 fetidppc_ctx->l2g_lambda = mat_ctx->l2g_lambda_only; 717 } else { 718 PetscCall(PetscObjectReference((PetscObject)mat_ctx->l2g_lambda)); 719 fetidppc_ctx->l2g_lambda = mat_ctx->l2g_lambda; 720 } 721 /* Dirichlet preconditioner */ 722 PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)fetimat)->prefix, "-pc_lumped", &lumped, NULL)); 723 if (!lumped) { 724 IS iV; 725 PetscBool discrete_harmonic = PETSC_FALSE; 726 727 PetscCall(PetscObjectQuery((PetscObject)fetidppc_ctx->pc, "__KSPFETIDP_iV", (PetscObject *)&iV)); 728 if (iV) PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)fetimat)->prefix, "-pc_discrete_harmonic", &discrete_harmonic, NULL)); 729 if (discrete_harmonic) { 730 KSP sksp; 731 PC pc; 732 PCBDDCSubSchurs sub_schurs = pcbddc->sub_schurs; 733 Mat A_II, A_IB, A_BI; 734 IS iP = NULL; 735 PetscBool isshell, reuse = PETSC_FALSE; 736 KSPType ksptype; 737 const char *prefix; 738 739 /* 740 We constructs a Schur complement for 741 742 | A_II A_ID | 743 | A_DI A_DD | 744 745 instead of 746 747 | A_II B^t_II A_ID | 748 | B_II -C_II B_ID | 749 | A_DI B^t_ID A_DD | 750 751 */ 752 if (sub_schurs && sub_schurs->reuse_solver) { 753 PetscCall(PetscObjectQuery((PetscObject)sub_schurs->A, "__KSPFETIDP_iP", (PetscObject *)&iP)); 754 if (iP) reuse = PETSC_TRUE; 755 } 756 if (!reuse) { 757 IS aB; 758 PetscInt nb; 759 PetscCall(ISGetLocalSize(pcis->is_B_local, &nb)); 760 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pcis->A_II), nb, 0, 1, &aB)); 761 PetscCall(MatCreateSubMatrix(pcis->A_II, iV, iV, MAT_INITIAL_MATRIX, &A_II)); 762 PetscCall(MatCreateSubMatrix(pcis->A_IB, iV, aB, MAT_INITIAL_MATRIX, &A_IB)); 763 PetscCall(MatCreateSubMatrix(pcis->A_BI, aB, iV, MAT_INITIAL_MATRIX, &A_BI)); 764 PetscCall(ISDestroy(&aB)); 765 } else { 766 PetscCall(MatCreateSubMatrix(sub_schurs->A, pcis->is_I_local, pcis->is_B_local, MAT_INITIAL_MATRIX, &A_IB)); 767 PetscCall(MatCreateSubMatrix(sub_schurs->A, pcis->is_B_local, pcis->is_I_local, MAT_INITIAL_MATRIX, &A_BI)); 768 PetscCall(PetscObjectReference((PetscObject)pcis->A_II)); 769 A_II = pcis->A_II; 770 } 771 PetscCall(MatCreateSchurComplement(A_II, A_II, A_IB, A_BI, pcis->A_BB, &fetidppc_ctx->S_j)); 772 773 /* propagate settings of solver */ 774 PetscCall(MatSchurComplementGetKSP(fetidppc_ctx->S_j, &sksp)); 775 PetscCall(KSPGetType(pcis->ksp_D, &ksptype)); 776 PetscCall(KSPSetType(sksp, ksptype)); 777 PetscCall(KSPGetPC(pcis->ksp_D, &pc)); 778 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCSHELL, &isshell)); 779 if (!isshell) { 780 MatSolverType solver; 781 PCType pctype; 782 783 PetscCall(PCGetType(pc, &pctype)); 784 PetscCall(PCFactorGetMatSolverType(pc, (MatSolverType *)&solver)); 785 PetscCall(KSPGetPC(sksp, &pc)); 786 PetscCall(PCSetType(pc, pctype)); 787 if (solver) PetscCall(PCFactorSetMatSolverType(pc, solver)); 788 } else { 789 PetscCall(KSPGetPC(sksp, &pc)); 790 PetscCall(PCSetType(pc, PCLU)); 791 } 792 PetscCall(MatDestroy(&A_II)); 793 PetscCall(MatDestroy(&A_IB)); 794 PetscCall(MatDestroy(&A_BI)); 795 PetscCall(MatGetOptionsPrefix(fetimat, &prefix)); 796 PetscCall(KSPSetOptionsPrefix(sksp, prefix)); 797 PetscCall(KSPAppendOptionsPrefix(sksp, "harmonic_")); 798 PetscCall(KSPSetFromOptions(sksp)); 799 if (reuse) { 800 PetscCall(KSPSetPC(sksp, sub_schurs->reuse_solver->interior_solver)); 801 PetscCall(PetscObjectIncrementTabLevel((PetscObject)sub_schurs->reuse_solver->interior_solver, (PetscObject)sksp, 0)); 802 } 803 } else { /* default Dirichlet preconditioner is pde-harmonic */ 804 PetscCall(MatCreateSchurComplement(pcis->A_II, pcis->A_II, pcis->A_IB, pcis->A_BI, pcis->A_BB, &fetidppc_ctx->S_j)); 805 PetscCall(MatSchurComplementSetKSP(fetidppc_ctx->S_j, pcis->ksp_D)); 806 } 807 } else { 808 PetscCall(PetscObjectReference((PetscObject)pcis->A_BB)); 809 fetidppc_ctx->S_j = pcis->A_BB; 810 } 811 /* saddle-point */ 812 if (mat_ctx->xPg) { 813 PetscCall(PetscObjectReference((PetscObject)mat_ctx->xPg)); 814 fetidppc_ctx->xPg = mat_ctx->xPg; 815 PetscCall(PetscObjectReference((PetscObject)mat_ctx->yPg)); 816 fetidppc_ctx->yPg = mat_ctx->yPg; 817 } 818 PetscFunctionReturn(PETSC_SUCCESS); 819 } 820 821 PetscErrorCode FETIDPMatMult_Kernel(Mat fetimat, Vec x, Vec y, PetscBool trans) 822 { 823 FETIDPMat_ctx mat_ctx; 824 PC_BDDC *pcbddc; 825 PC_IS *pcis; 826 827 PetscFunctionBegin; 828 PetscCall(MatShellGetContext(fetimat, &mat_ctx)); 829 pcis = (PC_IS *)mat_ctx->pc->data; 830 pcbddc = (PC_BDDC *)mat_ctx->pc->data; 831 /* Application of B_delta^T */ 832 PetscCall(VecSet(pcis->vec1_B, 0.)); 833 PetscCall(VecScatterBegin(mat_ctx->l2g_lambda, x, mat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE)); 834 PetscCall(VecScatterEnd(mat_ctx->l2g_lambda, x, mat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE)); 835 PetscCall(MatMultTranspose(mat_ctx->B_delta, mat_ctx->lambda_local, pcis->vec1_B)); 836 837 /* Add contribution from saddle point */ 838 if (mat_ctx->l2g_p) { 839 PetscCall(VecScatterBegin(mat_ctx->l2g_p, x, mat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE)); 840 PetscCall(VecScatterEnd(mat_ctx->l2g_p, x, mat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE)); 841 if (pcbddc->switch_static) { 842 if (trans) { 843 PetscCall(MatMultTranspose(mat_ctx->B_BI, mat_ctx->vP, pcis->vec1_D)); 844 } else { 845 PetscCall(MatMult(mat_ctx->Bt_BI, mat_ctx->vP, pcis->vec1_D)); 846 } 847 } 848 if (trans) { 849 PetscCall(MatMultTransposeAdd(mat_ctx->B_BB, mat_ctx->vP, pcis->vec1_B, pcis->vec1_B)); 850 } else { 851 PetscCall(MatMultAdd(mat_ctx->Bt_BB, mat_ctx->vP, pcis->vec1_B, pcis->vec1_B)); 852 } 853 } else { 854 if (pcbddc->switch_static) PetscCall(VecSet(pcis->vec1_D, 0.0)); 855 } 856 /* Application of \widetilde{S}^-1 */ 857 PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n)); 858 PetscCall(PCBDDCApplyInterfacePreconditioner(mat_ctx->pc, trans)); 859 PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n)); 860 PetscCall(VecSet(y, 0.0)); 861 /* Application of B_delta */ 862 PetscCall(MatMult(mat_ctx->B_delta, pcis->vec1_B, mat_ctx->lambda_local)); 863 /* Contribution from boundary pressures */ 864 if (mat_ctx->C) { 865 const PetscScalar *lx; 866 PetscScalar *ly; 867 868 /* pressure ordered first in the local part of x and y */ 869 PetscCall(VecGetArrayRead(x, &lx)); 870 PetscCall(VecGetArray(y, &ly)); 871 PetscCall(VecPlaceArray(mat_ctx->xPg, lx)); 872 PetscCall(VecPlaceArray(mat_ctx->yPg, ly)); 873 if (trans) { 874 PetscCall(MatMultTranspose(mat_ctx->C, mat_ctx->xPg, mat_ctx->yPg)); 875 } else { 876 PetscCall(MatMult(mat_ctx->C, mat_ctx->xPg, mat_ctx->yPg)); 877 } 878 PetscCall(VecResetArray(mat_ctx->xPg)); 879 PetscCall(VecResetArray(mat_ctx->yPg)); 880 PetscCall(VecRestoreArrayRead(x, &lx)); 881 PetscCall(VecRestoreArray(y, &ly)); 882 } 883 /* Add contribution from saddle point */ 884 if (mat_ctx->l2g_p) { 885 if (trans) { 886 PetscCall(MatMultTranspose(mat_ctx->Bt_BB, pcis->vec1_B, mat_ctx->vP)); 887 } else { 888 PetscCall(MatMult(mat_ctx->B_BB, pcis->vec1_B, mat_ctx->vP)); 889 } 890 if (pcbddc->switch_static) { 891 if (trans) { 892 PetscCall(MatMultTransposeAdd(mat_ctx->Bt_BI, pcis->vec1_D, mat_ctx->vP, mat_ctx->vP)); 893 } else { 894 PetscCall(MatMultAdd(mat_ctx->B_BI, pcis->vec1_D, mat_ctx->vP, mat_ctx->vP)); 895 } 896 } 897 PetscCall(VecScatterBegin(mat_ctx->l2g_p, mat_ctx->vP, y, ADD_VALUES, SCATTER_FORWARD)); 898 PetscCall(VecScatterEnd(mat_ctx->l2g_p, mat_ctx->vP, y, ADD_VALUES, SCATTER_FORWARD)); 899 } 900 PetscCall(VecScatterBegin(mat_ctx->l2g_lambda, mat_ctx->lambda_local, y, ADD_VALUES, SCATTER_FORWARD)); 901 PetscCall(VecScatterEnd(mat_ctx->l2g_lambda, mat_ctx->lambda_local, y, ADD_VALUES, SCATTER_FORWARD)); 902 PetscFunctionReturn(PETSC_SUCCESS); 903 } 904 905 PetscErrorCode FETIDPMatMult(Mat fetimat, Vec x, Vec y) 906 { 907 PetscFunctionBegin; 908 PetscCall(FETIDPMatMult_Kernel(fetimat, x, y, PETSC_FALSE)); 909 PetscFunctionReturn(PETSC_SUCCESS); 910 } 911 912 PetscErrorCode FETIDPMatMultTranspose(Mat fetimat, Vec x, Vec y) 913 { 914 PetscFunctionBegin; 915 PetscCall(FETIDPMatMult_Kernel(fetimat, x, y, PETSC_TRUE)); 916 PetscFunctionReturn(PETSC_SUCCESS); 917 } 918 919 PetscErrorCode FETIDPPCApply_Kernel(PC fetipc, Vec x, Vec y, PetscBool trans) 920 { 921 FETIDPPC_ctx pc_ctx; 922 PC_IS *pcis; 923 924 PetscFunctionBegin; 925 PetscCall(PCShellGetContext(fetipc, &pc_ctx)); 926 pcis = (PC_IS *)pc_ctx->pc->data; 927 /* Application of B_Ddelta^T */ 928 PetscCall(VecScatterBegin(pc_ctx->l2g_lambda, x, pc_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE)); 929 PetscCall(VecScatterEnd(pc_ctx->l2g_lambda, x, pc_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE)); 930 PetscCall(VecSet(pcis->vec2_B, 0.0)); 931 PetscCall(MatMultTranspose(pc_ctx->B_Ddelta, pc_ctx->lambda_local, pcis->vec2_B)); 932 /* Application of local Schur complement */ 933 if (trans) { 934 PetscCall(MatMultTranspose(pc_ctx->S_j, pcis->vec2_B, pcis->vec1_B)); 935 } else { 936 PetscCall(MatMult(pc_ctx->S_j, pcis->vec2_B, pcis->vec1_B)); 937 } 938 /* Application of B_Ddelta */ 939 PetscCall(MatMult(pc_ctx->B_Ddelta, pcis->vec1_B, pc_ctx->lambda_local)); 940 PetscCall(VecSet(y, 0.0)); 941 PetscCall(VecScatterBegin(pc_ctx->l2g_lambda, pc_ctx->lambda_local, y, ADD_VALUES, SCATTER_FORWARD)); 942 PetscCall(VecScatterEnd(pc_ctx->l2g_lambda, pc_ctx->lambda_local, y, ADD_VALUES, SCATTER_FORWARD)); 943 PetscFunctionReturn(PETSC_SUCCESS); 944 } 945 946 PetscErrorCode FETIDPPCApply(PC pc, Vec x, Vec y) 947 { 948 PetscFunctionBegin; 949 PetscCall(FETIDPPCApply_Kernel(pc, x, y, PETSC_FALSE)); 950 PetscFunctionReturn(PETSC_SUCCESS); 951 } 952 953 PetscErrorCode FETIDPPCApplyTranspose(PC pc, Vec x, Vec y) 954 { 955 PetscFunctionBegin; 956 PetscCall(FETIDPPCApply_Kernel(pc, x, y, PETSC_TRUE)); 957 PetscFunctionReturn(PETSC_SUCCESS); 958 } 959 960 PetscErrorCode FETIDPPCView(PC pc, PetscViewer viewer) 961 { 962 FETIDPPC_ctx pc_ctx; 963 PetscBool iascii; 964 PetscViewer sviewer; 965 966 PetscFunctionBegin; 967 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 968 if (iascii) { 969 PetscMPIInt rank; 970 PetscBool isschur, isshell; 971 972 PetscCall(PCShellGetContext(pc, &pc_ctx)); 973 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank)); 974 PetscCall(PetscObjectTypeCompare((PetscObject)pc_ctx->S_j, MATSCHURCOMPLEMENT, &isschur)); 975 if (isschur) { 976 PetscCall(PetscViewerASCIIPrintf(viewer, " Dirichlet preconditioner (just from rank 0)\n")); 977 } else { 978 PetscCall(PetscViewerASCIIPrintf(viewer, " Lumped preconditioner (just from rank 0)\n")); 979 } 980 PetscCall(PetscViewerGetSubViewer(viewer, PetscObjectComm((PetscObject)pc_ctx->S_j), &sviewer)); 981 if (rank == 0) { 982 PetscCall(PetscViewerPushFormat(sviewer, PETSC_VIEWER_ASCII_INFO)); 983 PetscCall(PetscViewerASCIIPushTab(sviewer)); 984 PetscCall(MatView(pc_ctx->S_j, sviewer)); 985 PetscCall(PetscViewerASCIIPopTab(sviewer)); 986 PetscCall(PetscViewerPopFormat(sviewer)); 987 } 988 PetscCall(PetscViewerRestoreSubViewer(viewer, PetscObjectComm((PetscObject)pc_ctx->S_j), &sviewer)); 989 PetscCall(PetscObjectTypeCompare((PetscObject)pc_ctx->B_Ddelta, MATSHELL, &isshell)); 990 if (isshell) { 991 BDdelta_DN ctx; 992 PetscCall(PetscViewerASCIIPrintf(viewer, " FETI-DP BDdelta: DB^t * (B D^-1 B^t)^-1 for deluxe scaling (just from rank 0)\n")); 993 PetscCall(MatShellGetContext(pc_ctx->B_Ddelta, &ctx)); 994 PetscCall(PetscViewerGetSubViewer(viewer, PetscObjectComm((PetscObject)pc_ctx->S_j), &sviewer)); 995 if (rank == 0) { 996 PetscInt tl; 997 998 PetscCall(PetscViewerASCIIGetTab(sviewer, &tl)); 999 PetscCall(PetscObjectSetTabLevel((PetscObject)ctx->kBD, tl)); 1000 PetscCall(KSPView(ctx->kBD, sviewer)); 1001 PetscCall(PetscViewerPushFormat(sviewer, PETSC_VIEWER_ASCII_INFO)); 1002 PetscCall(MatView(ctx->BD, sviewer)); 1003 PetscCall(PetscViewerPopFormat(sviewer)); 1004 } 1005 PetscCall(PetscViewerRestoreSubViewer(viewer, PetscObjectComm((PetscObject)pc_ctx->S_j), &sviewer)); 1006 } 1007 PetscCall(PetscViewerFlush(viewer)); 1008 } 1009 PetscFunctionReturn(PETSC_SUCCESS); 1010 } 1011