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