Actual source code: bddcfetidp.c
1: #include <petsc/private/pcbddcimpl.h>
2: #include <petsc/private/pcbddcprivateimpl.h>
3: #include <petscblaslapack.h>
5: static PetscErrorCode MatMult_BDdelta_deluxe_nonred(Mat A, Vec x, Vec y)
6: {
7: BDdelta_DN ctx;
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: }
17: static PetscErrorCode MatMultTranspose_BDdelta_deluxe_nonred(Mat A, Vec x, Vec y)
18: {
19: BDdelta_DN ctx;
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: }
29: static PetscErrorCode MatDestroy_BDdelta_deluxe_nonred(Mat A)
30: {
31: BDdelta_DN ctx;
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: }
42: PetscErrorCode PCBDDCCreateFETIDPMatContext(PC pc, FETIDPMat_ctx *fetidpmat_ctx)
43: {
44: FETIDPMat_ctx newctx;
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: }
55: PetscErrorCode PCBDDCCreateFETIDPPCContext(PC pc, FETIDPPC_ctx *fetidppc_ctx)
56: {
57: FETIDPPC_ctx newctx;
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: }
68: PetscErrorCode PCBDDCDestroyFETIDPMat(Mat A)
69: {
70: FETIDPMat_ctx mat_ctx;
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: }
99: PetscErrorCode PCBDDCDestroyFETIDPPC(PC pc)
100: {
101: FETIDPPC_ctx pc_ctx;
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: }
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;
141: /* saddlepoint */
142: ISLocalToGlobalMapping l2gmap_p;
143: PetscLayout play;
144: IS gP, pP;
145: PetscInt nPl, nPg, nPgl;
147: PetscFunctionBegin;
148: PetscCall(PetscObjectGetComm((PetscObject)(fetidpmat_ctx->pc), &comm));
149: PetscCallMPI(MPI_Comm_rank(comm, &rank));
150: PetscCallMPI(MPI_Comm_size(comm, &size));
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));
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;
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));
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));
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));
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));
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));
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: }
218: /* Default type of lagrange multipliers is non-redundant */
219: fully_redundant = fetidpmat_ctx->fully_redundant;
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;
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));
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));
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;
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));
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: }
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;
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: }
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: }
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: }
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));
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;
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;
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));
485: /* workspace allocation */
486: B_lwork = 0;
487: if (mss) {
488: PetscScalar dummy = 1;
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));
500: for (i = 0, cum = 0; i < sub_schurs->n_subs; i++) {
501: const PetscScalar *M;
502: PetscInt subset_size;
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));
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));
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));
544: if (BD2) {
545: ISLocalToGlobalMapping l2g;
546: Mat T, TA, *pT;
547: IS is;
548: PetscInt nl, N;
549: BDdelta_DN ctx;
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));
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));
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));
592: /* fetidpmat sizes */
593: fetidpmat_ctx->n += nPgl;
594: fetidpmat_ctx->N = fetidpmat_ctx->n_lambda + nPg;
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));
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;
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));
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;
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));
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));
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));
649: /* shift local to global indices for multipliers */
650: for (i = 0; i < n_local_lambda; i++) {
651: PetscInt ps;
652: PetscMPIInt owner;
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: }
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));
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));
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: }
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;
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;
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: }
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;
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;
739: /*
740: We constructs a Schur complement for
742: | A_II A_ID |
743: | A_DI A_DD |
745: instead of
747: | A_II B^t_II A_ID |
748: | B_II -C_II B_ID |
749: | A_DI B^t_ID A_DD |
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));
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;
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: }
821: static 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;
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));
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;
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: }
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: }
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: }
919: static PetscErrorCode FETIDPPCApply_Kernel(PC fetipc, Vec x, Vec y, PetscBool trans)
920: {
921: FETIDPPC_ctx pc_ctx;
922: PC_IS *pcis;
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: }
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: }
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: }
960: PetscErrorCode FETIDPPCView(PC pc, PetscViewer viewer)
961: {
962: FETIDPPC_ctx pc_ctx;
963: PetscBool iascii;
964: PetscViewer sviewer;
966: PetscFunctionBegin;
967: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
968: if (iascii) {
969: PetscMPIInt rank;
970: PetscBool isschur, isshell;
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;
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: }