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(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: }
101: PetscErrorCode PCBDDCDestroyFETIDPPC(PC pc)
102: {
103: FETIDPPC_ctx pc_ctx;
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: }
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;
144: /* saddlepoint */
145: ISLocalToGlobalMapping l2gmap_p;
146: PetscLayout play;
147: IS gP, pP;
148: PetscInt nPl, nPg, nPgl;
150: PetscFunctionBegin;
151: PetscCall(PetscObjectGetComm((PetscObject)fetidpmat_ctx->pc, &comm));
152: PetscCallMPI(MPI_Comm_rank(comm, &rank));
153: PetscCallMPI(MPI_Comm_size(comm, &size));
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));
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;
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));
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));
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));
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));
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));
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));
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));
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: }
229: /* Default type of lagrange multipliers is non-redundant */
230: fully_redundant = fetidpmat_ctx->fully_redundant;
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;
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));
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));
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];
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: }
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;
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));
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: }
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;
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: }
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: }
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));
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));
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;
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;
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));
520: /* workspace allocation */
521: B_lwork = 0;
522: if (mss) {
523: PetscScalar dummy = 1;
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));
535: for (i = 0, cum = 0; i < sub_schurs->n_subs; i++) {
536: const PetscScalar *M;
537: PetscInt subset_size;
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));
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));
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));
579: if (BD2) {
580: ISLocalToGlobalMapping l2g;
581: Mat T, TA, *pT;
582: IS is;
583: PetscInt nl, N;
584: BDdelta_DN ctx;
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));
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, (void (*)(void))MatMult_BDdelta_deluxe_nonred));
612: PetscCall(MatShellSetOperation(fetidpmat_ctx->B_Ddelta, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_BDdelta_deluxe_nonred));
613: PetscCall(MatShellSetOperation(fetidpmat_ctx->B_Ddelta, MATOP_DESTROY, (void (*)(void))MatDestroy_BDdelta_deluxe_nonred));
614: PetscCall(MatSetUp(fetidpmat_ctx->B_Ddelta));
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));
627: /* fetidpmat sizes */
628: fetidpmat_ctx->n += nPgl;
629: fetidpmat_ctx->N = fetidpmat_ctx->n_lambda + nPg;
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));
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;
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));
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;
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));
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));
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));
684: /* shift local to global indices for multipliers */
685: for (i = 0; i < n_local_lambda; i++) {
686: PetscInt ps;
687: PetscMPIInt owner;
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: }
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));
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));
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: }
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;
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;
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: }
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;
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;
774: /*
775: We constructs a Schur complement for
777: | A_II A_ID |
778: | A_DI A_DD |
780: instead of
782: | A_II B^t_II A_ID |
783: | B_II -C_II B_ID |
784: | A_DI B^t_ID A_DD |
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));
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;
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: }
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;
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));
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;
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: }
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: }
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: }
956: static PetscErrorCode FETIDPPCApply_Kernel(PC fetipc, Vec x, Vec y, PetscBool trans)
957: {
958: FETIDPPC_ctx pc_ctx;
959: PC_IS *pcis;
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: }
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: }
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: }
997: PetscErrorCode FETIDPPCView(PC pc, PetscViewer viewer)
998: {
999: FETIDPPC_ctx pc_ctx;
1000: PetscBool iascii;
1001: PetscViewer sviewer;
1003: PetscFunctionBegin;
1004: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1005: if (iascii) {
1006: PetscMPIInt rank;
1007: PetscBool isschur, isshell;
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;
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: }