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