Actual source code: bddcscalingbasic.c
1: #include <petsc/private/pcbddcimpl.h>
2: #include <petsc/private/pcbddcprivateimpl.h>
3: #include <petscblaslapack.h>
4: #include <../src/mat/impls/dense/seq/dense.h>
6: /* prototypes for deluxe functions */
7: static PetscErrorCode PCBDDCScalingCreate_Deluxe(PC);
8: static PetscErrorCode PCBDDCScalingDestroy_Deluxe(PC);
9: static PetscErrorCode PCBDDCScalingSetUp_Deluxe(PC);
10: static PetscErrorCode PCBDDCScalingSetUp_Deluxe_Private(PC);
11: static PetscErrorCode PCBDDCScalingReset_Deluxe_Solvers(PCBDDCDeluxeScaling);
13: static PetscErrorCode PCBDDCMatTransposeMatSolve_SeqDense(Mat A, Mat B, Mat X)
14: {
15: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
16: const PetscScalar *b;
17: PetscScalar *x;
18: PetscInt n;
19: PetscBLASInt nrhs, m;
20: PetscBool flg;
22: PetscFunctionBegin;
23: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
24: PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATSEQDENSE, MATMPIDENSE, NULL));
25: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix B must be MATDENSE matrix");
26: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &flg, MATSEQDENSE, MATMPIDENSE, NULL));
27: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix X must be MATDENSE matrix");
29: PetscCall(MatGetSize(B, NULL, &n));
30: PetscCall(PetscBLASIntCast(n, &nrhs));
31: PetscCall(MatDenseGetArrayRead(B, &b));
32: PetscCall(MatDenseGetArray(X, &x));
33: PetscCall(PetscArraycpy(x, b, m * nrhs));
34: PetscCall(MatDenseRestoreArrayRead(B, &b));
36: PetscCheck(A->factortype == MAT_FACTOR_LU, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only LU factor supported");
37: PetscCallLAPACKInfo("LAPACKgetrs", LAPACKgetrs_("T", &m, &nrhs, mat->v, &mat->lda, mat->pivots, x, &m, &info));
39: PetscCall(MatDenseRestoreArray(X, &x));
40: PetscCall(PetscLogFlops(nrhs * (2.0 * m * m - m)));
41: PetscFunctionReturn(PETSC_SUCCESS);
42: }
44: static PetscErrorCode PCBDDCScalingExtension_Basic(PC pc, Vec local_interface_vector, Vec global_vector)
45: {
46: PC_IS *pcis = (PC_IS *)pc->data;
47: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
49: PetscFunctionBegin;
50: /* Apply partition of unity */
51: PetscCall(VecPointwiseMult(pcbddc->work_scaling, pcis->D, local_interface_vector));
52: PetscCall(VecSet(global_vector, 0.0));
53: PetscCall(VecScatterBegin(pcis->global_to_B, pcbddc->work_scaling, global_vector, ADD_VALUES, SCATTER_REVERSE));
54: PetscCall(VecScatterEnd(pcis->global_to_B, pcbddc->work_scaling, global_vector, ADD_VALUES, SCATTER_REVERSE));
55: PetscFunctionReturn(PETSC_SUCCESS);
56: }
58: static PetscErrorCode PCBDDCScalingExtension_Deluxe(PC pc, Vec x, Vec y)
59: {
60: PC_IS *pcis = (PC_IS *)pc->data;
61: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
62: PCBDDCDeluxeScaling deluxe_ctx = pcbddc->deluxe_ctx;
64: PetscFunctionBegin;
65: PetscCall(VecSet(pcbddc->work_scaling, 0.0));
66: PetscCall(VecSet(y, 0.0));
67: if (deluxe_ctx->n_simple) { /* scale deluxe vertices using diagonal scaling */
68: PetscInt i;
69: const PetscScalar *array_x, *array_D;
70: PetscScalar *array;
71: PetscCall(VecGetArrayRead(x, &array_x));
72: PetscCall(VecGetArrayRead(pcis->D, &array_D));
73: PetscCall(VecGetArray(pcbddc->work_scaling, &array));
74: for (i = 0; i < deluxe_ctx->n_simple; i++) array[deluxe_ctx->idx_simple_B[i]] = array_x[deluxe_ctx->idx_simple_B[i]] * array_D[deluxe_ctx->idx_simple_B[i]];
75: PetscCall(VecRestoreArray(pcbddc->work_scaling, &array));
76: PetscCall(VecRestoreArrayRead(pcis->D, &array_D));
77: PetscCall(VecRestoreArrayRead(x, &array_x));
78: }
79: /* sequential part : all problems and Schur applications collapsed into a single matrix-vector multiplication or a matvec and a solve */
80: if (deluxe_ctx->seq_mat) {
81: PetscInt i;
82: for (i = 0; i < deluxe_ctx->seq_n; i++) {
83: if (deluxe_ctx->change) {
84: PetscCall(VecScatterBegin(deluxe_ctx->seq_scctx[i], x, deluxe_ctx->seq_work2[i], INSERT_VALUES, SCATTER_FORWARD));
85: PetscCall(VecScatterEnd(deluxe_ctx->seq_scctx[i], x, deluxe_ctx->seq_work2[i], INSERT_VALUES, SCATTER_FORWARD));
86: if (deluxe_ctx->change_with_qr) {
87: Mat change;
89: PetscCall(KSPGetOperators(deluxe_ctx->change[i], &change, NULL));
90: PetscCall(MatMultTranspose(change, deluxe_ctx->seq_work2[i], deluxe_ctx->seq_work1[i]));
91: } else {
92: PetscCall(KSPSolve(deluxe_ctx->change[i], deluxe_ctx->seq_work2[i], deluxe_ctx->seq_work1[i]));
93: }
94: } else {
95: PetscCall(VecScatterBegin(deluxe_ctx->seq_scctx[i], x, deluxe_ctx->seq_work1[i], INSERT_VALUES, SCATTER_FORWARD));
96: PetscCall(VecScatterEnd(deluxe_ctx->seq_scctx[i], x, deluxe_ctx->seq_work1[i], INSERT_VALUES, SCATTER_FORWARD));
97: }
98: PetscCall(MatMultTranspose(deluxe_ctx->seq_mat[i], deluxe_ctx->seq_work1[i], deluxe_ctx->seq_work2[i]));
99: if (deluxe_ctx->seq_mat_inv_sum[i]) {
100: PetscScalar *x;
102: PetscCall(VecGetArray(deluxe_ctx->seq_work2[i], &x));
103: PetscCall(VecPlaceArray(deluxe_ctx->seq_work1[i], x));
104: PetscCall(VecRestoreArray(deluxe_ctx->seq_work2[i], &x));
105: PetscCall(MatSolveTranspose(deluxe_ctx->seq_mat_inv_sum[i], deluxe_ctx->seq_work1[i], deluxe_ctx->seq_work2[i]));
106: PetscCall(VecResetArray(deluxe_ctx->seq_work1[i]));
107: }
108: if (deluxe_ctx->change) {
109: Mat change;
111: PetscCall(KSPGetOperators(deluxe_ctx->change[i], &change, NULL));
112: PetscCall(MatMult(change, deluxe_ctx->seq_work2[i], deluxe_ctx->seq_work1[i]));
113: PetscCall(VecScatterBegin(deluxe_ctx->seq_scctx[i], deluxe_ctx->seq_work1[i], pcbddc->work_scaling, INSERT_VALUES, SCATTER_REVERSE));
114: PetscCall(VecScatterEnd(deluxe_ctx->seq_scctx[i], deluxe_ctx->seq_work1[i], pcbddc->work_scaling, INSERT_VALUES, SCATTER_REVERSE));
115: } else {
116: PetscCall(VecScatterBegin(deluxe_ctx->seq_scctx[i], deluxe_ctx->seq_work2[i], pcbddc->work_scaling, INSERT_VALUES, SCATTER_REVERSE));
117: PetscCall(VecScatterEnd(deluxe_ctx->seq_scctx[i], deluxe_ctx->seq_work2[i], pcbddc->work_scaling, INSERT_VALUES, SCATTER_REVERSE));
118: }
119: }
120: }
121: /* put local boundary part in global vector */
122: PetscCall(VecScatterBegin(pcis->global_to_B, pcbddc->work_scaling, y, ADD_VALUES, SCATTER_REVERSE));
123: PetscCall(VecScatterEnd(pcis->global_to_B, pcbddc->work_scaling, y, ADD_VALUES, SCATTER_REVERSE));
124: PetscFunctionReturn(PETSC_SUCCESS);
125: }
127: PetscErrorCode PCBDDCScalingExtension(PC pc, Vec local_interface_vector, Vec global_vector)
128: {
129: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
131: PetscFunctionBegin;
135: PetscCheck(local_interface_vector != pcbddc->work_scaling, PETSC_COMM_SELF, PETSC_ERR_SUP, "Local vector cannot be pcbddc->work_scaling!");
136: PetscUseMethod(pc, "PCBDDCScalingExtension_C", (PC, Vec, Vec), (pc, local_interface_vector, global_vector));
137: PetscFunctionReturn(PETSC_SUCCESS);
138: }
140: static PetscErrorCode PCBDDCScalingRestriction_Basic(PC pc, Vec global_vector, Vec local_interface_vector)
141: {
142: PC_IS *pcis = (PC_IS *)pc->data;
144: PetscFunctionBegin;
145: PetscCall(VecScatterBegin(pcis->global_to_B, global_vector, local_interface_vector, INSERT_VALUES, SCATTER_FORWARD));
146: PetscCall(VecScatterEnd(pcis->global_to_B, global_vector, local_interface_vector, INSERT_VALUES, SCATTER_FORWARD));
147: /* Apply partition of unity */
148: PetscCall(VecPointwiseMult(local_interface_vector, pcis->D, local_interface_vector));
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
152: static PetscErrorCode PCBDDCScalingRestriction_Deluxe(PC pc, Vec x, Vec y)
153: {
154: PC_IS *pcis = (PC_IS *)pc->data;
155: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
156: PCBDDCDeluxeScaling deluxe_ctx = pcbddc->deluxe_ctx;
158: PetscFunctionBegin;
159: /* get local boundary part of global vector */
160: PetscCall(VecScatterBegin(pcis->global_to_B, x, y, INSERT_VALUES, SCATTER_FORWARD));
161: PetscCall(VecScatterEnd(pcis->global_to_B, x, y, INSERT_VALUES, SCATTER_FORWARD));
162: if (deluxe_ctx->n_simple) { /* scale deluxe vertices using diagonal scaling */
163: PetscInt i;
164: PetscScalar *array_y;
165: const PetscScalar *array_D;
166: PetscCall(VecGetArray(y, &array_y));
167: PetscCall(VecGetArrayRead(pcis->D, &array_D));
168: for (i = 0; i < deluxe_ctx->n_simple; i++) array_y[deluxe_ctx->idx_simple_B[i]] *= array_D[deluxe_ctx->idx_simple_B[i]];
169: PetscCall(VecRestoreArrayRead(pcis->D, &array_D));
170: PetscCall(VecRestoreArray(y, &array_y));
171: }
172: /* sequential part : all problems and Schur applications collapsed into a single matrix-vector multiplication or a matvec and a solve */
173: if (deluxe_ctx->seq_mat) {
174: PetscInt i;
175: for (i = 0; i < deluxe_ctx->seq_n; i++) {
176: if (deluxe_ctx->change) {
177: Mat change;
179: PetscCall(VecScatterBegin(deluxe_ctx->seq_scctx[i], y, deluxe_ctx->seq_work2[i], INSERT_VALUES, SCATTER_FORWARD));
180: PetscCall(VecScatterEnd(deluxe_ctx->seq_scctx[i], y, deluxe_ctx->seq_work2[i], INSERT_VALUES, SCATTER_FORWARD));
181: PetscCall(KSPGetOperators(deluxe_ctx->change[i], &change, NULL));
182: PetscCall(MatMultTranspose(change, deluxe_ctx->seq_work2[i], deluxe_ctx->seq_work1[i]));
183: } else {
184: PetscCall(VecScatterBegin(deluxe_ctx->seq_scctx[i], y, deluxe_ctx->seq_work1[i], INSERT_VALUES, SCATTER_FORWARD));
185: PetscCall(VecScatterEnd(deluxe_ctx->seq_scctx[i], y, deluxe_ctx->seq_work1[i], INSERT_VALUES, SCATTER_FORWARD));
186: }
187: if (deluxe_ctx->seq_mat_inv_sum[i]) {
188: PetscScalar *x;
190: PetscCall(VecGetArray(deluxe_ctx->seq_work1[i], &x));
191: PetscCall(VecPlaceArray(deluxe_ctx->seq_work2[i], x));
192: PetscCall(VecRestoreArray(deluxe_ctx->seq_work1[i], &x));
193: PetscCall(MatSolve(deluxe_ctx->seq_mat_inv_sum[i], deluxe_ctx->seq_work1[i], deluxe_ctx->seq_work2[i]));
194: PetscCall(VecResetArray(deluxe_ctx->seq_work2[i]));
195: }
196: PetscCall(MatMult(deluxe_ctx->seq_mat[i], deluxe_ctx->seq_work1[i], deluxe_ctx->seq_work2[i]));
197: if (deluxe_ctx->change) {
198: if (deluxe_ctx->change_with_qr) {
199: Mat change;
201: PetscCall(KSPGetOperators(deluxe_ctx->change[i], &change, NULL));
202: PetscCall(MatMult(change, deluxe_ctx->seq_work2[i], deluxe_ctx->seq_work1[i]));
203: } else {
204: PetscCall(KSPSolveTranspose(deluxe_ctx->change[i], deluxe_ctx->seq_work2[i], deluxe_ctx->seq_work1[i]));
205: }
206: PetscCall(VecScatterBegin(deluxe_ctx->seq_scctx[i], deluxe_ctx->seq_work1[i], y, INSERT_VALUES, SCATTER_REVERSE));
207: PetscCall(VecScatterEnd(deluxe_ctx->seq_scctx[i], deluxe_ctx->seq_work1[i], y, INSERT_VALUES, SCATTER_REVERSE));
208: } else {
209: PetscCall(VecScatterBegin(deluxe_ctx->seq_scctx[i], deluxe_ctx->seq_work2[i], y, INSERT_VALUES, SCATTER_REVERSE));
210: PetscCall(VecScatterEnd(deluxe_ctx->seq_scctx[i], deluxe_ctx->seq_work2[i], y, INSERT_VALUES, SCATTER_REVERSE));
211: }
212: }
213: }
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
217: PetscErrorCode PCBDDCScalingRestriction(PC pc, Vec global_vector, Vec local_interface_vector)
218: {
219: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
221: PetscFunctionBegin;
225: PetscCheck(local_interface_vector != pcbddc->work_scaling, PETSC_COMM_SELF, PETSC_ERR_SUP, "Local vector cannot be pcbddc->work_scaling!");
226: PetscUseMethod(pc, "PCBDDCScalingRestriction_C", (PC, Vec, Vec), (pc, global_vector, local_interface_vector));
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
230: PetscErrorCode PCBDDCScalingSetUp(PC pc)
231: {
232: PC_IS *pcis = (PC_IS *)pc->data;
233: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
235: PetscFunctionBegin;
237: PetscCall(PetscLogEventBegin(PC_BDDC_Scaling[pcbddc->current_level], pc, 0, 0, 0));
238: /* create work vector for the operator */
239: PetscCall(VecDestroy(&pcbddc->work_scaling));
240: PetscCall(VecDuplicate(pcis->vec1_B, &pcbddc->work_scaling));
241: /* always rebuild pcis->D */
242: if (pcis->use_stiffness_scaling) {
243: PetscScalar *a;
244: PetscInt i, n;
246: PetscCall(MatGetDiagonal(pcbddc->local_mat, pcis->vec1_N));
247: PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec1_N, pcis->D, INSERT_VALUES, SCATTER_FORWARD));
248: PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec1_N, pcis->D, INSERT_VALUES, SCATTER_FORWARD));
249: PetscCall(VecAbs(pcis->D));
250: PetscCall(VecGetLocalSize(pcis->D, &n));
251: PetscCall(VecGetArray(pcis->D, &a));
252: for (i = 0; i < n; i++)
253: if (PetscAbsScalar(a[i]) < PETSC_SMALL) a[i] = 1.0;
254: PetscCall(VecRestoreArray(pcis->D, &a));
255: }
256: PetscCall(VecSet(pcis->vec1_global, 0.0));
257: PetscCall(VecScatterBegin(pcis->global_to_B, pcis->D, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
258: PetscCall(VecScatterEnd(pcis->global_to_B, pcis->D, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
259: PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_global, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
260: PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_global, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
261: PetscCall(VecPointwiseDivide(pcis->D, pcis->D, pcis->vec1_B));
262: /* now setup */
263: if (pcbddc->use_deluxe_scaling) {
264: if (!pcbddc->deluxe_ctx) PetscCall(PCBDDCScalingCreate_Deluxe(pc));
265: PetscCall(PCBDDCScalingSetUp_Deluxe(pc));
266: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCScalingRestriction_C", PCBDDCScalingRestriction_Deluxe));
267: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCScalingExtension_C", PCBDDCScalingExtension_Deluxe));
268: } else {
269: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCScalingRestriction_C", PCBDDCScalingRestriction_Basic));
270: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCScalingExtension_C", PCBDDCScalingExtension_Basic));
271: }
272: PetscCall(PetscLogEventEnd(PC_BDDC_Scaling[pcbddc->current_level], pc, 0, 0, 0));
274: /* test */
275: if (pcbddc->dbg_flag) {
276: Mat B0_B = NULL;
277: Vec B0_Bv = NULL, B0_Bv2 = NULL;
278: Vec vec2_global;
279: PetscViewer viewer = pcbddc->dbg_viewer;
280: PetscReal error;
282: /* extension -> from local to parallel */
283: PetscCall(VecSet(pcis->vec1_global, 0.0));
284: PetscCall(VecSetRandom(pcis->vec1_B, NULL));
285: PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_B, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
286: PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_B, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
287: PetscCall(VecDuplicate(pcis->vec1_global, &vec2_global));
288: PetscCall(VecCopy(pcis->vec1_global, vec2_global));
289: PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_global, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
290: PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_global, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
291: if (pcbddc->benign_n) {
292: IS is_dummy;
294: PetscCall(ISCreateStride(PETSC_COMM_SELF, pcbddc->benign_n, 0, 1, &is_dummy));
295: PetscCall(MatCreateSubMatrix(pcbddc->benign_B0, is_dummy, pcis->is_B_local, MAT_INITIAL_MATRIX, &B0_B));
296: PetscCall(ISDestroy(&is_dummy));
297: PetscCall(MatCreateVecs(B0_B, NULL, &B0_Bv));
298: PetscCall(VecDuplicate(B0_Bv, &B0_Bv2));
299: PetscCall(MatMult(B0_B, pcis->vec1_B, B0_Bv));
300: }
301: PetscCall(PCBDDCScalingExtension(pc, pcis->vec1_B, pcis->vec1_global));
302: if (pcbddc->benign_saddle_point) {
303: error = 0.;
304: PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_global, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
305: PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_global, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
306: if (pcbddc->benign_n) {
307: PetscCall(MatMult(B0_B, pcis->vec1_B, B0_Bv2));
308: PetscCall(VecAXPY(B0_Bv, -1.0, B0_Bv2));
309: PetscCall(VecNorm(B0_Bv, NORM_INFINITY, &error));
310: }
311: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &error, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)pc)));
312: PetscCall(PetscViewerASCIIPrintf(viewer, "Error benign extension %1.14e\n", (double)error));
313: }
314: PetscCall(VecAXPY(pcis->vec1_global, -1.0, vec2_global));
315: PetscCall(VecNorm(pcis->vec1_global, NORM_INFINITY, &error));
316: PetscCall(PetscViewerASCIIPrintf(viewer, "Error scaling extension %1.14e\n", (double)error));
317: PetscCall(VecDestroy(&vec2_global));
319: /* restriction -> from parallel to local */
320: PetscCall(VecSet(pcis->vec1_global, 0.0));
321: PetscCall(VecSetRandom(pcis->vec1_B, NULL));
322: PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_B, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
323: PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_B, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
324: PetscCall(PCBDDCScalingRestriction(pc, pcis->vec1_global, pcis->vec1_B));
325: PetscCall(VecScale(pcis->vec1_B, -1.0));
326: PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_B, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
327: PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_B, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
328: PetscCall(VecNorm(pcis->vec1_global, NORM_INFINITY, &error));
329: PetscCall(PetscViewerASCIIPrintf(viewer, "Error scaling restriction %1.14e\n", (double)error));
330: PetscCall(MatDestroy(&B0_B));
331: PetscCall(VecDestroy(&B0_Bv));
332: PetscCall(VecDestroy(&B0_Bv2));
333: }
334: PetscFunctionReturn(PETSC_SUCCESS);
335: }
337: PetscErrorCode PCBDDCScalingDestroy(PC pc)
338: {
339: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
341: PetscFunctionBegin;
342: if (pcbddc->deluxe_ctx) PetscCall(PCBDDCScalingDestroy_Deluxe(pc));
343: PetscCall(VecDestroy(&pcbddc->work_scaling));
344: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCScalingRestriction_C", NULL));
345: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCScalingExtension_C", NULL));
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }
349: static PetscErrorCode PCBDDCScalingCreate_Deluxe(PC pc)
350: {
351: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
352: PCBDDCDeluxeScaling deluxe_ctx;
354: PetscFunctionBegin;
355: PetscCall(PetscNew(&deluxe_ctx));
356: pcbddc->deluxe_ctx = deluxe_ctx;
357: PetscFunctionReturn(PETSC_SUCCESS);
358: }
360: static PetscErrorCode PCBDDCScalingDestroy_Deluxe(PC pc)
361: {
362: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
364: PetscFunctionBegin;
365: PetscCall(PCBDDCScalingReset_Deluxe_Solvers(pcbddc->deluxe_ctx));
366: PetscCall(PetscFree(pcbddc->deluxe_ctx));
367: PetscFunctionReturn(PETSC_SUCCESS);
368: }
370: static PetscErrorCode PCBDDCScalingReset_Deluxe_Solvers(PCBDDCDeluxeScaling deluxe_ctx)
371: {
372: PetscInt i;
374: PetscFunctionBegin;
375: PetscCall(PetscFree(deluxe_ctx->idx_simple_B));
376: deluxe_ctx->n_simple = 0;
377: for (i = 0; i < deluxe_ctx->seq_n; i++) {
378: PetscCall(VecScatterDestroy(&deluxe_ctx->seq_scctx[i]));
379: PetscCall(VecDestroy(&deluxe_ctx->seq_work1[i]));
380: PetscCall(VecDestroy(&deluxe_ctx->seq_work2[i]));
381: PetscCall(MatDestroy(&deluxe_ctx->seq_mat[i]));
382: PetscCall(MatDestroy(&deluxe_ctx->seq_mat_inv_sum[i]));
383: }
384: PetscCall(PetscFree5(deluxe_ctx->seq_scctx, deluxe_ctx->seq_work1, deluxe_ctx->seq_work2, deluxe_ctx->seq_mat, deluxe_ctx->seq_mat_inv_sum));
385: PetscCall(PetscFree(deluxe_ctx->workspace));
386: deluxe_ctx->seq_n = 0;
387: PetscFunctionReturn(PETSC_SUCCESS);
388: }
390: static PetscErrorCode PCBDDCScalingSetUp_Deluxe(PC pc)
391: {
392: PC_IS *pcis = (PC_IS *)pc->data;
393: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
394: PCBDDCDeluxeScaling deluxe_ctx = pcbddc->deluxe_ctx;
395: PCBDDCSubSchurs sub_schurs = pcbddc->sub_schurs;
397: PetscFunctionBegin;
398: /* reset data structures if the topology has changed */
399: if (pcbddc->recompute_topography) PetscCall(PCBDDCScalingReset_Deluxe_Solvers(deluxe_ctx));
401: /* Compute data structures to solve sequential problems */
402: PetscCall(PCBDDCScalingSetUp_Deluxe_Private(pc));
404: /* diagonal scaling on interface dofs not contained in cc */
405: if (sub_schurs->is_vertices || sub_schurs->is_dir) {
406: PetscInt n_com, n_dir;
407: n_com = 0;
408: if (sub_schurs->is_vertices) PetscCall(ISGetLocalSize(sub_schurs->is_vertices, &n_com));
409: n_dir = 0;
410: if (sub_schurs->is_dir) PetscCall(ISGetLocalSize(sub_schurs->is_dir, &n_dir));
411: if (!deluxe_ctx->n_simple) {
412: deluxe_ctx->n_simple = n_dir + n_com;
413: PetscCall(PetscMalloc1(deluxe_ctx->n_simple, &deluxe_ctx->idx_simple_B));
414: if (n_com) {
415: PetscInt nmap;
416: const PetscInt *idxs;
418: PetscCall(ISGetIndices(sub_schurs->is_vertices, &idxs));
419: PetscCall(ISGlobalToLocalMappingApply(pcis->BtoNmap, IS_GTOLM_DROP, n_com, idxs, &nmap, deluxe_ctx->idx_simple_B));
420: PetscCheck(nmap == n_com, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error when mapping simply scaled dofs (is_vertices)! %" PetscInt_FMT " != %" PetscInt_FMT, nmap, n_com);
421: PetscCall(ISRestoreIndices(sub_schurs->is_vertices, &idxs));
422: }
423: if (n_dir) {
424: PetscInt nmap;
425: const PetscInt *idxs;
427: PetscCall(ISGetIndices(sub_schurs->is_dir, &idxs));
428: PetscCall(ISGlobalToLocalMappingApply(pcis->BtoNmap, IS_GTOLM_DROP, n_dir, idxs, &nmap, deluxe_ctx->idx_simple_B + n_com));
429: PetscCheck(nmap == n_dir, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error when mapping simply scaled dofs (sub_schurs->is_dir)! %" PetscInt_FMT " != %" PetscInt_FMT, nmap, n_dir);
430: PetscCall(ISRestoreIndices(sub_schurs->is_dir, &idxs));
431: }
432: PetscCall(PetscSortInt(deluxe_ctx->n_simple, deluxe_ctx->idx_simple_B));
433: } else {
434: PetscCheck(deluxe_ctx->n_simple == n_dir + n_com, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of simply scaled dofs %" PetscInt_FMT " is different from the previous one computed %" PetscInt_FMT, n_dir + n_com, deluxe_ctx->n_simple);
435: }
436: } else {
437: deluxe_ctx->n_simple = 0;
438: deluxe_ctx->idx_simple_B = NULL;
439: }
440: PetscFunctionReturn(PETSC_SUCCESS);
441: }
443: static PetscErrorCode PCBDDCScalingSetUp_Deluxe_Private(PC pc)
444: {
445: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
446: PCBDDCDeluxeScaling deluxe_ctx = pcbddc->deluxe_ctx;
447: PCBDDCSubSchurs sub_schurs = pcbddc->sub_schurs;
448: PetscScalar *matdata, *matdata2;
449: PetscInt i, max_subset_size, cum, cum2;
450: const PetscInt *idxs;
451: PetscBool newsetup = PETSC_FALSE;
453: PetscFunctionBegin;
454: PetscCheck(sub_schurs, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Missing PCBDDCSubSchurs");
455: if (!sub_schurs->n_subs) PetscFunctionReturn(PETSC_SUCCESS);
457: /* Allocate arrays for subproblems */
458: if (!deluxe_ctx->seq_n) {
459: deluxe_ctx->seq_n = sub_schurs->n_subs;
460: PetscCall(PetscCalloc5(deluxe_ctx->seq_n, &deluxe_ctx->seq_scctx, deluxe_ctx->seq_n, &deluxe_ctx->seq_work1, deluxe_ctx->seq_n, &deluxe_ctx->seq_work2, deluxe_ctx->seq_n, &deluxe_ctx->seq_mat, deluxe_ctx->seq_n, &deluxe_ctx->seq_mat_inv_sum));
461: newsetup = PETSC_TRUE;
462: } else PetscCheck(deluxe_ctx->seq_n == sub_schurs->n_subs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of deluxe subproblems %" PetscInt_FMT " is different from the sub_schurs %" PetscInt_FMT, deluxe_ctx->seq_n, sub_schurs->n_subs);
464: /* the change of basis is just a reference to sub_schurs->change (if any) */
465: deluxe_ctx->change = sub_schurs->change;
466: deluxe_ctx->change_with_qr = sub_schurs->change_with_qr;
468: /* Create objects for deluxe */
469: max_subset_size = 0;
470: for (i = 0; i < sub_schurs->n_subs; i++) {
471: PetscInt subset_size;
472: PetscCall(ISGetLocalSize(sub_schurs->is_subs[i], &subset_size));
473: max_subset_size = PetscMax(subset_size, max_subset_size);
474: }
475: if (newsetup) PetscCall(PetscMalloc1(2 * max_subset_size, &deluxe_ctx->workspace));
476: cum = cum2 = 0;
477: PetscCall(ISGetIndices(sub_schurs->is_Ej_all, &idxs));
478: PetscCall(MatSeqAIJGetArray(sub_schurs->S_Ej_all, &matdata));
479: PetscCall(MatSeqAIJGetArray(sub_schurs->sum_S_Ej_all, &matdata2));
480: for (i = 0; i < deluxe_ctx->seq_n; i++) {
481: PetscInt subset_size;
483: PetscCall(ISGetLocalSize(sub_schurs->is_subs[i], &subset_size));
484: if (newsetup) {
485: IS sub;
486: /* work vectors */
487: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, subset_size, deluxe_ctx->workspace, &deluxe_ctx->seq_work1[i]));
488: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, subset_size, deluxe_ctx->workspace + subset_size, &deluxe_ctx->seq_work2[i]));
490: /* scatters */
491: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, subset_size, idxs + cum, PETSC_COPY_VALUES, &sub));
492: PetscCall(VecScatterCreate(pcbddc->work_scaling, sub, deluxe_ctx->seq_work1[i], NULL, &deluxe_ctx->seq_scctx[i]));
493: PetscCall(ISDestroy(&sub));
494: }
496: /* S_E_j */
497: PetscCall(MatDestroy(&deluxe_ctx->seq_mat[i]));
498: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, subset_size, subset_size, matdata + cum2, &deluxe_ctx->seq_mat[i]));
500: /* \sum_k S^k_E_j */
501: PetscCall(MatDestroy(&deluxe_ctx->seq_mat_inv_sum[i]));
502: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, subset_size, subset_size, matdata2 + cum2, &deluxe_ctx->seq_mat_inv_sum[i]));
503: PetscCall(MatSetOption(deluxe_ctx->seq_mat_inv_sum[i], MAT_SPD, sub_schurs->is_posdef));
504: PetscCall(MatSetOption(deluxe_ctx->seq_mat_inv_sum[i], MAT_HERMITIAN, sub_schurs->is_hermitian));
505: switch (sub_schurs->mat_factor_type) {
506: case MAT_FACTOR_CHOLESKY:
507: PetscCall(MatCholeskyFactor(deluxe_ctx->seq_mat_inv_sum[i], NULL, NULL));
508: break;
509: case MAT_FACTOR_LU:
510: PetscCall(MatLUFactor(deluxe_ctx->seq_mat_inv_sum[i], NULL, NULL, NULL));
511: break;
512: default:
513: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported factor type %s", MatFactorTypes[sub_schurs->mat_factor_type]);
514: }
515: if (pcbddc->deluxe_singlemat) {
516: Mat X, Y;
517: if (!sub_schurs->is_hermitian) {
518: PetscCall(MatTranspose(deluxe_ctx->seq_mat[i], MAT_INITIAL_MATRIX, &X));
519: } else {
520: PetscCall(PetscObjectReference((PetscObject)deluxe_ctx->seq_mat[i]));
521: X = deluxe_ctx->seq_mat[i];
522: }
523: PetscCall(MatDuplicate(X, MAT_DO_NOT_COPY_VALUES, &Y));
524: if (!sub_schurs->is_hermitian) {
525: PetscCall(PCBDDCMatTransposeMatSolve_SeqDense(deluxe_ctx->seq_mat_inv_sum[i], X, Y));
526: } else {
527: PetscCall(MatMatSolve(deluxe_ctx->seq_mat_inv_sum[i], X, Y));
528: }
529: PetscCall(MatDestroy(&deluxe_ctx->seq_mat_inv_sum[i]));
530: PetscCall(MatDestroy(&deluxe_ctx->seq_mat[i]));
531: PetscCall(MatDestroy(&X));
532: if (deluxe_ctx->change) {
533: Mat C, CY;
534: PetscCheck(deluxe_ctx->change_with_qr, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only QR based change of basis");
535: PetscCall(KSPGetOperators(deluxe_ctx->change[i], &C, NULL));
536: PetscCall(MatMatMult(C, Y, MAT_INITIAL_MATRIX, PETSC_CURRENT, &CY));
537: PetscCall(MatMatTransposeMult(CY, C, MAT_REUSE_MATRIX, PETSC_CURRENT, &Y));
538: PetscCall(MatDestroy(&CY));
539: PetscCall(MatProductClear(Y)); /* clear internal matproduct structure of Y since CY is destroyed */
540: }
541: PetscCall(MatTranspose(Y, MAT_INPLACE_MATRIX, &Y));
542: deluxe_ctx->seq_mat[i] = Y;
543: }
544: cum += subset_size;
545: cum2 += subset_size * subset_size;
546: }
547: PetscCall(ISRestoreIndices(sub_schurs->is_Ej_all, &idxs));
548: PetscCall(MatSeqAIJRestoreArray(sub_schurs->S_Ej_all, &matdata));
549: PetscCall(MatSeqAIJRestoreArray(sub_schurs->sum_S_Ej_all, &matdata2));
550: if (pcbddc->deluxe_singlemat) {
551: deluxe_ctx->change = NULL;
552: deluxe_ctx->change_with_qr = PETSC_FALSE;
553: }
555: if (deluxe_ctx->change && !deluxe_ctx->change_with_qr) {
556: for (i = 0; i < deluxe_ctx->seq_n; i++) {
557: if (newsetup) {
558: PC pc;
560: PetscCall(KSPGetPC(deluxe_ctx->change[i], &pc));
561: PetscCall(PCSetType(pc, PCLU));
562: PetscCall(KSPSetFromOptions(deluxe_ctx->change[i]));
563: }
564: PetscCall(KSPSetUp(deluxe_ctx->change[i]));
565: }
566: }
567: PetscFunctionReturn(PETSC_SUCCESS);
568: }