Actual source code: bddcschurs.c
1: #include <petsc/private/pcbddcimpl.h>
2: #include <petsc/private/pcbddcprivateimpl.h>
3: #include <../src/mat/impls/dense/seq/dense.h>
4: #include <petscblaslapack.h>
6: static inline PetscErrorCode PCBDDCAdjGetNextLayer_Private(PetscInt *, PetscInt, PetscBT, PetscInt *, PetscInt *, PetscInt *);
7: static PetscErrorCode PCBDDCComputeExplicitSchur(Mat, PetscBool, MatReuse, Mat *);
8: static PetscErrorCode PCBDDCReuseSolvers_Interior(PC, Vec, Vec);
9: static PetscErrorCode PCBDDCReuseSolvers_Correction(PC, Vec, Vec);
11: /* if v2 is not present, correction is done in-place */
12: PetscErrorCode PCBDDCReuseSolversBenignAdapt(PCBDDCReuseSolvers ctx, Vec v, Vec v2, PetscBool sol, PetscBool full)
13: {
14: PetscScalar *array;
15: PetscScalar *array2;
17: PetscFunctionBegin;
18: if (!ctx->benign_n) PetscFunctionReturn(PETSC_SUCCESS);
19: if (sol && full) {
20: PetscInt n_I, size_schur;
22: /* get sizes */
23: PetscCall(MatGetSize(ctx->benign_csAIB, &size_schur, NULL));
24: PetscCall(VecGetSize(v, &n_I));
25: n_I = n_I - size_schur;
26: /* get schur sol from array */
27: PetscCall(VecGetArray(v, &array));
28: PetscCall(VecPlaceArray(ctx->benign_dummy_schur_vec, array + n_I));
29: PetscCall(VecRestoreArray(v, &array));
30: /* apply interior sol correction */
31: PetscCall(MatMultTranspose(ctx->benign_csAIB, ctx->benign_dummy_schur_vec, ctx->benign_corr_work));
32: PetscCall(VecResetArray(ctx->benign_dummy_schur_vec));
33: PetscCall(MatMultAdd(ctx->benign_AIIm1ones, ctx->benign_corr_work, v, v));
34: }
35: if (v2) {
36: PetscInt nl;
38: PetscCall(VecGetArrayRead(v, (const PetscScalar **)&array));
39: PetscCall(VecGetLocalSize(v2, &nl));
40: PetscCall(VecGetArray(v2, &array2));
41: PetscCall(PetscArraycpy(array2, array, nl));
42: } else {
43: PetscCall(VecGetArray(v, &array));
44: array2 = array;
45: }
46: if (!sol) { /* change rhs */
47: PetscInt n;
48: for (n = 0; n < ctx->benign_n; n++) {
49: PetscScalar sum = 0.;
50: const PetscInt *cols;
51: PetscInt nz, i;
53: PetscCall(ISGetLocalSize(ctx->benign_zerodiag_subs[n], &nz));
54: PetscCall(ISGetIndices(ctx->benign_zerodiag_subs[n], &cols));
55: for (i = 0; i < nz - 1; i++) sum += array[cols[i]];
56: #if PetscDefined(USE_COMPLEX)
57: sum = -(PetscRealPart(sum) / nz + PETSC_i * (PetscImaginaryPart(sum) / nz));
58: #else
59: sum = -sum / nz;
60: #endif
61: for (i = 0; i < nz - 1; i++) array2[cols[i]] += sum;
62: ctx->benign_save_vals[n] = array2[cols[nz - 1]];
63: array2[cols[nz - 1]] = sum;
64: PetscCall(ISRestoreIndices(ctx->benign_zerodiag_subs[n], &cols));
65: }
66: } else {
67: for (PetscInt n = 0; n < ctx->benign_n; n++) {
68: PetscScalar sum = 0.;
69: const PetscInt *cols;
70: PetscInt nz;
71: PetscCall(ISGetLocalSize(ctx->benign_zerodiag_subs[n], &nz));
72: PetscCall(ISGetIndices(ctx->benign_zerodiag_subs[n], &cols));
73: for (PetscInt i = 0; i < nz - 1; i++) sum += array[cols[i]];
74: #if PetscDefined(USE_COMPLEX)
75: sum = -(PetscRealPart(sum) / nz + PETSC_i * (PetscImaginaryPart(sum) / nz));
76: #else
77: sum = -sum / nz;
78: #endif
79: for (PetscInt i = 0; i < nz - 1; i++) array2[cols[i]] += sum;
80: array2[cols[nz - 1]] = ctx->benign_save_vals[n];
81: PetscCall(ISRestoreIndices(ctx->benign_zerodiag_subs[n], &cols));
82: }
83: }
84: if (v2) {
85: PetscCall(VecRestoreArrayRead(v, (const PetscScalar **)&array));
86: PetscCall(VecRestoreArray(v2, &array2));
87: } else {
88: PetscCall(VecRestoreArray(v, &array));
89: }
90: if (!sol && full) {
91: Vec usedv;
92: PetscInt n_I, size_schur;
94: /* get sizes */
95: PetscCall(MatGetSize(ctx->benign_csAIB, &size_schur, NULL));
96: PetscCall(VecGetSize(v, &n_I));
97: n_I = n_I - size_schur;
98: /* compute schur rhs correction */
99: if (v2) {
100: usedv = v2;
101: } else {
102: usedv = v;
103: }
104: /* apply schur rhs correction */
105: PetscCall(MatMultTranspose(ctx->benign_AIIm1ones, usedv, ctx->benign_corr_work));
106: PetscCall(VecGetArrayRead(usedv, (const PetscScalar **)&array));
107: PetscCall(VecPlaceArray(ctx->benign_dummy_schur_vec, array + n_I));
108: PetscCall(VecRestoreArrayRead(usedv, (const PetscScalar **)&array));
109: PetscCall(MatMultAdd(ctx->benign_csAIB, ctx->benign_corr_work, ctx->benign_dummy_schur_vec, ctx->benign_dummy_schur_vec));
110: PetscCall(VecResetArray(ctx->benign_dummy_schur_vec));
111: }
112: PetscFunctionReturn(PETSC_SUCCESS);
113: }
115: static PetscErrorCode PCBDDCReuseSolvers_Solve_Private(PC pc, Vec rhs, Vec sol, PetscBool transpose, PetscBool full)
116: {
117: PCBDDCReuseSolvers ctx;
118: PetscBool copy = PETSC_FALSE;
120: PetscFunctionBegin;
121: PetscCall(PCShellGetContext(pc, &ctx));
122: if (full) {
123: PetscCall(MatMumpsSetIcntl(ctx->F, 26, -1));
124: #if PetscDefined(HAVE_MKL_PARDISO)
125: PetscCall(MatMkl_PardisoSetCntl(ctx->F, 70, 0));
126: #endif
127: copy = ctx->has_vertices;
128: } else { /* interior solver */
129: PetscCall(MatMumpsSetIcntl(ctx->F, 26, 0));
130: #if PetscDefined(HAVE_MKL_PARDISO)
131: PetscCall(MatMkl_PardisoSetCntl(ctx->F, 70, 1));
132: #endif
133: copy = PETSC_TRUE;
134: }
135: /* copy rhs into factored matrix workspace */
136: if (copy) {
137: PetscInt n;
138: PetscScalar *array, *array_solver;
140: PetscCall(VecGetLocalSize(rhs, &n));
141: PetscCall(VecGetArrayRead(rhs, (const PetscScalar **)&array));
142: PetscCall(VecGetArray(ctx->rhs, &array_solver));
143: PetscCall(PetscArraycpy(array_solver, array, n));
144: PetscCall(VecRestoreArray(ctx->rhs, &array_solver));
145: PetscCall(VecRestoreArrayRead(rhs, (const PetscScalar **)&array));
147: PetscCall(PCBDDCReuseSolversBenignAdapt(ctx, ctx->rhs, NULL, PETSC_FALSE, full));
148: if (transpose) {
149: PetscCall(MatSolveTranspose(ctx->F, ctx->rhs, ctx->sol));
150: } else {
151: PetscCall(MatSolve(ctx->F, ctx->rhs, ctx->sol));
152: }
153: PetscCall(PCBDDCReuseSolversBenignAdapt(ctx, ctx->sol, NULL, PETSC_TRUE, full));
155: /* get back data to caller worskpace */
156: PetscCall(VecGetArrayRead(ctx->sol, (const PetscScalar **)&array_solver));
157: PetscCall(VecGetArray(sol, &array));
158: PetscCall(PetscArraycpy(array, array_solver, n));
159: PetscCall(VecRestoreArray(sol, &array));
160: PetscCall(VecRestoreArrayRead(ctx->sol, (const PetscScalar **)&array_solver));
161: } else {
162: if (ctx->benign_n) {
163: PetscCall(PCBDDCReuseSolversBenignAdapt(ctx, rhs, ctx->rhs, PETSC_FALSE, full));
164: if (transpose) {
165: PetscCall(MatSolveTranspose(ctx->F, ctx->rhs, sol));
166: } else {
167: PetscCall(MatSolve(ctx->F, ctx->rhs, sol));
168: }
169: PetscCall(PCBDDCReuseSolversBenignAdapt(ctx, sol, NULL, PETSC_TRUE, full));
170: } else {
171: if (transpose) {
172: PetscCall(MatSolveTranspose(ctx->F, rhs, sol));
173: } else {
174: PetscCall(MatSolve(ctx->F, rhs, sol));
175: }
176: }
177: }
178: /* restore defaults */
179: PetscCall(MatMumpsSetIcntl(ctx->F, 26, -1));
180: #if PetscDefined(HAVE_MKL_PARDISO)
181: PetscCall(MatMkl_PardisoSetCntl(ctx->F, 70, 0));
182: #endif
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }
186: static PetscErrorCode PCBDDCReuseSolvers_Correction(PC pc, Vec rhs, Vec sol)
187: {
188: PetscFunctionBegin;
189: PetscCall(PCBDDCReuseSolvers_Solve_Private(pc, rhs, sol, PETSC_FALSE, PETSC_TRUE));
190: PetscFunctionReturn(PETSC_SUCCESS);
191: }
193: static PetscErrorCode PCBDDCReuseSolvers_CorrectionTranspose(PC pc, Vec rhs, Vec sol)
194: {
195: PetscFunctionBegin;
196: PetscCall(PCBDDCReuseSolvers_Solve_Private(pc, rhs, sol, PETSC_TRUE, PETSC_TRUE));
197: PetscFunctionReturn(PETSC_SUCCESS);
198: }
200: static PetscErrorCode PCBDDCReuseSolvers_Interior(PC pc, Vec rhs, Vec sol)
201: {
202: PetscFunctionBegin;
203: PetscCall(PCBDDCReuseSolvers_Solve_Private(pc, rhs, sol, PETSC_FALSE, PETSC_FALSE));
204: PetscFunctionReturn(PETSC_SUCCESS);
205: }
207: static PetscErrorCode PCBDDCReuseSolvers_InteriorTranspose(PC pc, Vec rhs, Vec sol)
208: {
209: PetscFunctionBegin;
210: PetscCall(PCBDDCReuseSolvers_Solve_Private(pc, rhs, sol, PETSC_TRUE, PETSC_FALSE));
211: PetscFunctionReturn(PETSC_SUCCESS);
212: }
214: static PetscErrorCode PCBDDCReuseSolvers_View(PC pc, PetscViewer viewer)
215: {
216: PCBDDCReuseSolvers ctx;
217: PetscBool isascii;
219: PetscFunctionBegin;
220: PetscCall(PCShellGetContext(pc, &ctx));
221: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
222: if (isascii) PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
223: PetscCall(MatView(ctx->F, viewer));
224: if (isascii) PetscCall(PetscViewerPopFormat(viewer));
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }
228: static PetscErrorCode PCBDDCReuseSolversReset(PCBDDCReuseSolvers reuse)
229: {
230: PetscInt i;
232: PetscFunctionBegin;
233: PetscCall(MatDestroy(&reuse->F));
234: PetscCall(VecDestroy(&reuse->sol));
235: PetscCall(VecDestroy(&reuse->rhs));
236: PetscCall(PCDestroy(&reuse->interior_solver));
237: PetscCall(PCDestroy(&reuse->correction_solver));
238: PetscCall(ISDestroy(&reuse->is_R));
239: PetscCall(ISDestroy(&reuse->is_B));
240: PetscCall(VecScatterDestroy(&reuse->correction_scatter_B));
241: PetscCall(VecDestroy(&reuse->sol_B));
242: PetscCall(VecDestroy(&reuse->rhs_B));
243: for (i = 0; i < reuse->benign_n; i++) PetscCall(ISDestroy(&reuse->benign_zerodiag_subs[i]));
244: PetscCall(PetscFree(reuse->benign_zerodiag_subs));
245: PetscCall(PetscFree(reuse->benign_save_vals));
246: PetscCall(MatDestroy(&reuse->benign_csAIB));
247: PetscCall(MatDestroy(&reuse->benign_AIIm1ones));
248: PetscCall(VecDestroy(&reuse->benign_corr_work));
249: PetscCall(VecDestroy(&reuse->benign_dummy_schur_vec));
250: PetscFunctionReturn(PETSC_SUCCESS);
251: }
253: static PetscErrorCode PCBDDCReuseSolvers_Destroy(PC pc)
254: {
255: PCBDDCReuseSolvers ctx;
257: PetscFunctionBegin;
258: PetscCall(PCShellGetContext(pc, &ctx));
259: PetscCall(PCBDDCReuseSolversReset(ctx));
260: PetscCall(PetscFree(ctx));
261: PetscCall(PCShellSetContext(pc, NULL));
262: PetscFunctionReturn(PETSC_SUCCESS);
263: }
265: static PetscErrorCode PCBDDCComputeExplicitSchur(Mat M, PetscBool issym, MatReuse reuse, Mat *S)
266: {
267: Mat B, C, D, Bd, Cd, AinvBd;
268: KSP ksp;
269: PC pc;
270: PetscBool isLU, isILU, isCHOL, Bdense, Cdense;
271: PetscReal fill = 2.0;
272: PetscInt n_I;
273: PetscMPIInt size;
275: PetscFunctionBegin;
276: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)M), &size));
277: PetscCheck(size == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for parallel matrices");
278: if (reuse == MAT_REUSE_MATRIX) {
279: PetscBool Sdense;
281: PetscCall(PetscObjectTypeCompare((PetscObject)*S, MATSEQDENSE, &Sdense));
282: PetscCheck(Sdense, PetscObjectComm((PetscObject)M), PETSC_ERR_SUP, "S should dense");
283: }
284: PetscCall(MatSchurComplementGetSubMatrices(M, NULL, NULL, &B, &C, &D));
285: PetscCall(MatSchurComplementGetKSP(M, &ksp));
286: PetscCall(KSPGetPC(ksp, &pc));
287: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLU, &isLU));
288: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCILU, &isILU));
289: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCCHOLESKY, &isCHOL));
290: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATSEQDENSE, &Bdense));
291: PetscCall(PetscObjectTypeCompare((PetscObject)C, MATSEQDENSE, &Cdense));
292: PetscCall(MatGetSize(B, &n_I, NULL));
293: if (n_I) {
294: if (!Bdense) {
295: PetscCall(MatConvert(B, MATSEQDENSE, MAT_INITIAL_MATRIX, &Bd));
296: } else {
297: Bd = B;
298: }
300: if (isLU || isILU || isCHOL) {
301: Mat fact;
302: PetscCall(KSPSetUp(ksp));
303: PetscCall(PCFactorGetMatrix(pc, &fact));
304: PetscCall(MatDuplicate(Bd, MAT_DO_NOT_COPY_VALUES, &AinvBd));
305: PetscCall(MatMatSolve(fact, Bd, AinvBd));
306: } else {
307: PetscBool ex = PETSC_TRUE;
309: if (ex) {
310: Mat Ainvd;
312: PetscCall(PCComputeOperator(pc, MATDENSE, &Ainvd));
313: PetscCall(MatMatMult(Ainvd, Bd, MAT_INITIAL_MATRIX, fill, &AinvBd));
314: PetscCall(MatDestroy(&Ainvd));
315: } else {
316: Vec sol, rhs;
317: PetscScalar *arrayrhs, *arraysol;
318: PetscInt i, nrhs, n;
320: PetscCall(MatDuplicate(Bd, MAT_DO_NOT_COPY_VALUES, &AinvBd));
321: PetscCall(MatGetSize(Bd, &n, &nrhs));
322: PetscCall(MatDenseGetArray(Bd, &arrayrhs));
323: PetscCall(MatDenseGetArray(AinvBd, &arraysol));
324: PetscCall(KSPGetSolution(ksp, &sol));
325: PetscCall(KSPGetRhs(ksp, &rhs));
326: for (i = 0; i < nrhs; i++) {
327: PetscCall(VecPlaceArray(rhs, arrayrhs + i * n));
328: PetscCall(VecPlaceArray(sol, arraysol + i * n));
329: PetscCall(KSPSolve(ksp, rhs, sol));
330: PetscCall(VecResetArray(rhs));
331: PetscCall(VecResetArray(sol));
332: }
333: PetscCall(MatDenseRestoreArray(Bd, &arrayrhs));
334: PetscCall(MatDenseRestoreArray(AinvBd, &arrayrhs));
335: }
336: }
337: if (!Bdense & !issym) PetscCall(MatDestroy(&Bd));
339: if (!issym) {
340: if (!Cdense) {
341: PetscCall(MatConvert(C, MATSEQDENSE, MAT_INITIAL_MATRIX, &Cd));
342: } else {
343: Cd = C;
344: }
345: PetscCall(MatMatMult(Cd, AinvBd, reuse, fill, S));
346: if (!Cdense) PetscCall(MatDestroy(&Cd));
347: } else {
348: PetscCall(MatTransposeMatMult(Bd, AinvBd, reuse, fill, S));
349: if (!Bdense) PetscCall(MatDestroy(&Bd));
350: }
351: PetscCall(MatDestroy(&AinvBd));
352: }
354: if (D) {
355: Mat Dd;
356: PetscBool Ddense;
358: PetscCall(PetscObjectTypeCompare((PetscObject)D, MATSEQDENSE, &Ddense));
359: if (!Ddense) {
360: PetscCall(MatConvert(D, MATSEQDENSE, MAT_INITIAL_MATRIX, &Dd));
361: } else {
362: Dd = D;
363: }
364: if (n_I) PetscCall(MatAYPX(*S, -1.0, Dd, SAME_NONZERO_PATTERN));
365: else {
366: if (reuse == MAT_INITIAL_MATRIX) {
367: PetscCall(MatDuplicate(Dd, MAT_COPY_VALUES, S));
368: } else {
369: PetscCall(MatCopy(Dd, *S, SAME_NONZERO_PATTERN));
370: }
371: }
372: if (!Ddense) PetscCall(MatDestroy(&Dd));
373: } else {
374: PetscCall(MatScale(*S, -1.0));
375: }
376: PetscFunctionReturn(PETSC_SUCCESS);
377: }
379: PetscErrorCode PCBDDCSubSchursSetUp(PCBDDCSubSchurs sub_schurs, Mat Ain, Mat Sin, PetscBool exact_schur, PetscInt xadj[], PetscInt adjncy[], PetscInt nlayers, Vec scaling, PetscBool compute_Stilda, PetscBool reuse_solvers, PetscBool benign_trick, PetscInt benign_n, PetscInt benign_p0_lidx[], IS benign_zerodiag_subs[], Mat change, IS change_primal)
380: {
381: Mat F, A_II, A_IB, A_BI, A_BB, AE_II;
382: Mat S_all;
383: Vec gstash, lstash;
384: VecScatter sstash;
385: IS is_I, is_I_layer;
386: IS all_subsets, all_subsets_mult, all_subsets_n;
387: PetscScalar *stasharray, *Bwork;
388: PetscInt *all_local_idx_N, *all_local_subid_N = NULL;
389: PetscInt *auxnum1, *auxnum2;
390: PetscInt *local_subs = sub_schurs->graph->local_subs;
391: PetscInt i, subset_size, max_subset_size, n_local_subs = sub_schurs->graph->n_local_subs;
392: PetscInt n_B, extra, local_size, global_size;
393: PetscInt local_stash_size;
394: PetscBLASInt B_N, B_lwork, *pivots;
395: MPI_Comm comm_n;
396: PetscBool deluxe = PETSC_TRUE;
397: PetscBool use_potr = PETSC_FALSE, use_sytr = PETSC_FALSE;
398: PetscViewer matl_dbg_viewer = NULL;
399: PetscBool flg, multi_element = sub_schurs->graph->multi_element;
401: PetscFunctionBegin;
402: PetscCall(MatDestroy(&sub_schurs->A));
403: PetscCall(MatDestroy(&sub_schurs->S));
404: if (Ain) {
405: PetscCall(PetscObjectReference((PetscObject)Ain));
406: sub_schurs->A = Ain;
407: }
409: PetscCall(PetscObjectReference((PetscObject)Sin));
410: sub_schurs->S = Sin;
411: if (sub_schurs->schur_explicit) sub_schurs->schur_explicit = (PetscBool)(!!sub_schurs->A);
413: /* preliminary checks */
414: PetscCheck(sub_schurs->schur_explicit || !compute_Stilda, PetscObjectComm((PetscObject)sub_schurs->l2gmap), PETSC_ERR_SUP, "Adaptive selection of constraints requires MUMPS and/or MKL_PARDISO");
416: if (benign_trick) sub_schurs->is_posdef = PETSC_FALSE;
418: /* debug (MATLAB) */
419: if (sub_schurs->debug) {
420: PetscMPIInt size, rank;
421: PetscInt nr, *print_schurs_ranks, print_schurs = PETSC_FALSE;
423: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)sub_schurs->l2gmap), &size));
424: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)sub_schurs->l2gmap), &rank));
425: nr = size;
426: PetscCall(PetscMalloc1(nr, &print_schurs_ranks));
427: PetscOptionsBegin(PetscObjectComm((PetscObject)sub_schurs->l2gmap), sub_schurs->prefix, "BDDC sub_schurs options", "PC");
428: PetscCall(PetscOptionsIntArray("-sub_schurs_debug_ranks", "Ranks to debug (all if the option is not used)", NULL, print_schurs_ranks, &nr, &flg));
429: if (!flg) print_schurs = PETSC_TRUE;
430: else {
431: print_schurs = PETSC_FALSE;
432: for (i = 0; i < nr; i++)
433: if (print_schurs_ranks[i] == rank) {
434: print_schurs = PETSC_TRUE;
435: break;
436: }
437: }
438: PetscOptionsEnd();
439: PetscCall(PetscFree(print_schurs_ranks));
440: if (print_schurs) {
441: char filename[256];
443: PetscCall(PetscSNPrintf(filename, sizeof(filename), "sub_schurs_Schur_r%d.m", PetscGlobalRank));
444: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_SELF, filename, &matl_dbg_viewer));
445: PetscCall(PetscViewerPushFormat(matl_dbg_viewer, PETSC_VIEWER_ASCII_MATLAB));
446: }
447: }
449: /* DEBUG: turn on/off multi-element code path */
450: PetscCall(PetscOptionsGetBool(NULL, sub_schurs->prefix, "-sub_schurs_multielement_code", &multi_element, NULL));
451: if (n_local_subs == 0) multi_element = PETSC_FALSE;
453: /* restrict work on active processes */
454: if (sub_schurs->restrict_comm) {
455: PetscSubcomm subcomm;
456: PetscMPIInt color, rank;
458: color = 0;
459: if (!sub_schurs->n_subs) color = 1; /* this can happen if we are in a multilevel case or if the subdomain is disconnected */
460: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)sub_schurs->l2gmap), &rank));
461: PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)sub_schurs->l2gmap), &subcomm));
462: PetscCall(PetscSubcommSetNumber(subcomm, 2));
463: PetscCall(PetscSubcommSetTypeGeneral(subcomm, color, rank));
464: PetscCall(PetscCommDuplicate(PetscSubcommChild(subcomm), &comm_n, NULL));
465: PetscCall(PetscSubcommDestroy(&subcomm));
466: if (!sub_schurs->n_subs) {
467: PetscCall(PetscCommDestroy(&comm_n));
468: PetscFunctionReturn(PETSC_SUCCESS);
469: }
470: } else {
471: PetscCall(PetscCommDuplicate(PetscObjectComm((PetscObject)sub_schurs->l2gmap), &comm_n, NULL));
472: }
474: /* get Schur complement matrices */
475: if (!sub_schurs->schur_explicit) {
476: Mat tA_IB, tA_BI, tA_BB;
477: PetscBool isseqsbaij;
478: PetscCall(MatSchurComplementGetSubMatrices(sub_schurs->S, &A_II, NULL, &tA_IB, &tA_BI, &tA_BB));
479: PetscCall(PetscObjectTypeCompare((PetscObject)tA_BB, MATSEQSBAIJ, &isseqsbaij));
480: if (isseqsbaij) {
481: PetscCall(MatConvert(tA_BB, MATSEQAIJ, MAT_INITIAL_MATRIX, &A_BB));
482: PetscCall(MatConvert(tA_IB, MATSEQAIJ, MAT_INITIAL_MATRIX, &A_IB));
483: PetscCall(MatConvert(tA_BI, MATSEQAIJ, MAT_INITIAL_MATRIX, &A_BI));
484: } else {
485: PetscCall(PetscObjectReference((PetscObject)tA_BB));
486: A_BB = tA_BB;
487: PetscCall(PetscObjectReference((PetscObject)tA_IB));
488: A_IB = tA_IB;
489: PetscCall(PetscObjectReference((PetscObject)tA_BI));
490: A_BI = tA_BI;
491: }
492: } else {
493: A_II = NULL;
494: A_IB = NULL;
495: A_BI = NULL;
496: A_BB = NULL;
497: }
498: S_all = NULL;
500: /* determine interior problems */
501: PetscCall(ISGetLocalSize(sub_schurs->is_I, &i));
502: if (nlayers >= 0 && i) { /* Interior problems can be different from the original one */
503: PetscBT touched;
504: const PetscInt *idx_B;
505: PetscInt n_I, n_B, n_local_dofs, n_prev_added, j, layer, *local_numbering;
507: PetscCheck(xadj, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot request layering without adjacency");
508: /* get sizes */
509: PetscCall(ISGetLocalSize(sub_schurs->is_I, &n_I));
510: PetscCall(ISGetLocalSize(sub_schurs->is_B, &n_B));
512: PetscCall(PetscMalloc1(n_I + n_B, &local_numbering));
513: PetscCall(PetscBTCreate(n_I + n_B, &touched));
514: PetscCall(PetscBTMemzero(n_I + n_B, touched));
516: /* all boundary dofs must be skipped when adding layers */
517: PetscCall(ISGetIndices(sub_schurs->is_B, &idx_B));
518: for (j = 0; j < n_B; j++) PetscCall(PetscBTSet(touched, idx_B[j]));
519: PetscCall(PetscArraycpy(local_numbering, idx_B, n_B));
520: PetscCall(ISRestoreIndices(sub_schurs->is_B, &idx_B));
522: /* add prescribed number of layers of dofs */
523: n_local_dofs = n_B;
524: n_prev_added = n_B;
525: for (layer = 0; layer < nlayers; layer++) {
526: PetscInt n_added = 0;
527: if (n_local_dofs == n_I + n_B) break;
528: PetscCheck(n_local_dofs <= n_I + n_B, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error querying layer %" PetscInt_FMT ". Out of bound access (%" PetscInt_FMT " > %" PetscInt_FMT ")", layer, n_local_dofs, n_I + n_B);
529: PetscCall(PCBDDCAdjGetNextLayer_Private(local_numbering + n_local_dofs, n_prev_added, touched, xadj, adjncy, &n_added));
530: n_prev_added = n_added;
531: n_local_dofs += n_added;
532: if (!n_added) break;
533: }
534: PetscCall(PetscBTDestroy(&touched));
536: /* IS for I layer dofs in original numbering */
537: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)sub_schurs->is_I), n_local_dofs - n_B, local_numbering + n_B, PETSC_COPY_VALUES, &is_I_layer));
538: PetscCall(PetscFree(local_numbering));
539: PetscCall(ISSort(is_I_layer));
540: /* IS for I layer dofs in I numbering */
541: if (!sub_schurs->schur_explicit) {
542: ISLocalToGlobalMapping ItoNmap;
543: PetscCall(ISLocalToGlobalMappingCreateIS(sub_schurs->is_I, &ItoNmap));
544: PetscCall(ISGlobalToLocalMappingApplyIS(ItoNmap, IS_GTOLM_DROP, is_I_layer, &is_I));
545: PetscCall(ISLocalToGlobalMappingDestroy(&ItoNmap));
547: /* II block */
548: PetscCall(MatCreateSubMatrix(A_II, is_I, is_I, MAT_INITIAL_MATRIX, &AE_II));
549: }
550: } else {
551: PetscInt n_I;
553: /* IS for I dofs in original numbering */
554: PetscCall(PetscObjectReference((PetscObject)sub_schurs->is_I));
555: is_I_layer = sub_schurs->is_I;
557: /* IS for I dofs in I numbering (strided 1) */
558: if (!sub_schurs->schur_explicit) {
559: PetscCall(ISGetSize(sub_schurs->is_I, &n_I));
560: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)sub_schurs->is_I), n_I, 0, 1, &is_I));
562: /* II block is the same */
563: PetscCall(PetscObjectReference((PetscObject)A_II));
564: AE_II = A_II;
565: }
566: }
568: /* Get info on subset sizes and sum of all subsets sizes */
569: max_subset_size = 0;
570: local_size = 0;
571: for (i = 0; i < sub_schurs->n_subs; i++) {
572: PetscCall(ISGetLocalSize(sub_schurs->is_subs[i], &subset_size));
573: max_subset_size = PetscMax(subset_size, max_subset_size);
574: local_size += subset_size;
575: }
577: /* Work arrays for local indices */
578: extra = 0;
579: PetscCall(ISGetLocalSize(sub_schurs->is_B, &n_B));
580: if (sub_schurs->schur_explicit && is_I_layer) PetscCall(ISGetLocalSize(is_I_layer, &extra));
581: PetscCall(PetscMalloc1(n_B + extra, &all_local_idx_N));
582: if (multi_element) PetscCall(PetscMalloc1(n_B + extra, &all_local_subid_N));
583: if (extra) {
584: const PetscInt *idxs;
585: PetscCall(ISGetIndices(is_I_layer, &idxs));
586: PetscCall(PetscArraycpy(all_local_idx_N, idxs, extra));
587: if (multi_element)
588: for (PetscInt j = 0; j < extra; j++) all_local_subid_N[j] = local_subs[idxs[j]];
589: PetscCall(ISRestoreIndices(is_I_layer, &idxs));
590: }
591: PetscCall(PetscMalloc1(sub_schurs->n_subs, &auxnum1));
592: PetscCall(PetscMalloc1(sub_schurs->n_subs, &auxnum2));
594: /* Get local indices in local numbering */
595: local_size = 0;
596: local_stash_size = 0;
597: for (i = 0; i < sub_schurs->n_subs; i++) {
598: const PetscInt *idxs;
600: PetscCall(ISGetLocalSize(sub_schurs->is_subs[i], &subset_size));
601: PetscCall(ISGetIndices(sub_schurs->is_subs[i], &idxs));
602: /* start (smallest in global ordering) and multiplicity */
603: auxnum1[i] = idxs[0];
604: auxnum2[i] = subset_size * subset_size;
605: /* subset indices in local numbering */
606: PetscCall(PetscArraycpy(all_local_idx_N + local_size + extra, idxs, subset_size));
607: if (multi_element)
608: for (PetscInt j = 0; j < subset_size; j++) all_local_subid_N[j + local_size + extra] = local_subs[idxs[j]];
609: PetscCall(ISRestoreIndices(sub_schurs->is_subs[i], &idxs));
610: local_size += subset_size;
611: local_stash_size += subset_size * subset_size;
612: }
614: /* allocate extra workspace needed only for GETRI or SYTRF when inverting the blocks or the entire Schur complement */
615: use_potr = use_sytr = PETSC_FALSE;
616: if (benign_trick || (sub_schurs->is_hermitian && sub_schurs->is_posdef)) {
617: use_potr = PETSC_TRUE;
618: } else if (sub_schurs->is_symmetric) {
619: use_sytr = PETSC_TRUE;
620: }
621: if (local_size && !use_potr && compute_Stilda) {
622: PetscScalar lwork, dummyscalar = 0.;
623: PetscBLASInt dummyint = 0;
625: B_lwork = -1;
626: PetscCall(PetscBLASIntCast(local_size, &B_N));
627: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
628: if (use_sytr) {
629: PetscCallLAPACKInfo("LAPACKsytrf", LAPACKsytrf_("L", &B_N, &dummyscalar, &B_N, &dummyint, &lwork, &B_lwork, &info));
630: } else {
631: PetscCallLAPACKInfo("LAPACKgetri", LAPACKgetri_(&B_N, &dummyscalar, &B_N, &dummyint, &lwork, &B_lwork, &info));
632: }
633: PetscCall(PetscFPTrapPop());
634: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(lwork), &B_lwork));
635: PetscCall(PetscMalloc2(B_lwork, &Bwork, B_N, &pivots));
636: } else {
637: Bwork = NULL;
638: pivots = NULL;
639: }
641: /* prepare data for summing up properly schurs on subsets */
642: PetscCall(ISCreateGeneral(comm_n, sub_schurs->n_subs, auxnum1, PETSC_OWN_POINTER, &all_subsets_n));
643: PetscCall(ISLocalToGlobalMappingApplyIS(sub_schurs->l2gmap, all_subsets_n, &all_subsets));
644: PetscCall(ISDestroy(&all_subsets_n));
645: PetscCall(ISCreateGeneral(comm_n, sub_schurs->n_subs, auxnum2, PETSC_OWN_POINTER, &all_subsets_mult));
646: PetscCall(ISRenumber(all_subsets, all_subsets_mult, &global_size, &all_subsets_n));
647: PetscCall(ISDestroy(&all_subsets));
648: PetscCall(ISDestroy(&all_subsets_mult));
649: PetscCall(ISGetLocalSize(all_subsets_n, &i));
650: PetscCheck(i == local_stash_size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid size of new subset! %" PetscInt_FMT " != %" PetscInt_FMT, i, local_stash_size);
651: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, local_stash_size, NULL, &lstash));
652: PetscCall(VecCreateMPI(comm_n, PETSC_DECIDE, global_size, &gstash));
653: PetscCall(VecScatterCreate(lstash, NULL, gstash, all_subsets_n, &sstash));
654: PetscCall(ISDestroy(&all_subsets_n));
656: /* subset indices in local boundary numbering */
657: if (!sub_schurs->is_Ej_all) {
658: PetscInt *all_local_idx_B;
660: PetscCall(PetscMalloc1(local_size, &all_local_idx_B));
661: PetscCall(ISGlobalToLocalMappingApply(sub_schurs->BtoNmap, IS_GTOLM_DROP, local_size, all_local_idx_N + extra, &subset_size, all_local_idx_B));
662: PetscCheck(subset_size == local_size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in sub_schurs serial (BtoNmap)! %" PetscInt_FMT " != %" PetscInt_FMT, subset_size, local_size);
663: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, local_size, all_local_idx_B, PETSC_OWN_POINTER, &sub_schurs->is_Ej_all));
664: }
666: if (change) {
667: ISLocalToGlobalMapping BtoS;
668: IS change_primal_B;
669: IS change_primal_all;
671: PetscCheck(!sub_schurs->change_primal_sub, PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
672: PetscCheck(!sub_schurs->change, PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
673: PetscCall(PetscMalloc1(sub_schurs->n_subs, &sub_schurs->change_primal_sub));
674: for (i = 0; i < sub_schurs->n_subs; i++) {
675: ISLocalToGlobalMapping NtoS;
676: PetscCall(ISLocalToGlobalMappingCreateIS(sub_schurs->is_subs[i], &NtoS));
677: PetscCall(ISGlobalToLocalMappingApplyIS(NtoS, IS_GTOLM_DROP, change_primal, &sub_schurs->change_primal_sub[i]));
678: PetscCall(ISLocalToGlobalMappingDestroy(&NtoS));
679: }
680: PetscCall(ISGlobalToLocalMappingApplyIS(sub_schurs->BtoNmap, IS_GTOLM_DROP, change_primal, &change_primal_B));
681: PetscCall(ISLocalToGlobalMappingCreateIS(sub_schurs->is_Ej_all, &BtoS));
682: PetscCall(ISGlobalToLocalMappingApplyIS(BtoS, IS_GTOLM_DROP, change_primal_B, &change_primal_all));
683: PetscCall(ISLocalToGlobalMappingDestroy(&BtoS));
684: PetscCall(ISDestroy(&change_primal_B));
685: PetscCall(PetscMalloc1(sub_schurs->n_subs, &sub_schurs->change));
686: for (i = 0; i < sub_schurs->n_subs; i++) {
687: Mat change_sub;
689: PetscCall(ISGetLocalSize(sub_schurs->is_subs[i], &subset_size));
690: PetscCall(KSPCreate(PETSC_COMM_SELF, &sub_schurs->change[i]));
691: PetscCall(KSPSetNestLevel(sub_schurs->change[i], 1)); /* do not seem to have direct access to a PC from which to get the level of nests */
692: PetscCall(KSPSetType(sub_schurs->change[i], KSPPREONLY));
693: if (!sub_schurs->change_with_qr) {
694: PetscCall(MatCreateSubMatrix(change, sub_schurs->is_subs[i], sub_schurs->is_subs[i], MAT_INITIAL_MATRIX, &change_sub));
695: } else {
696: Mat change_subt;
697: PetscCall(MatCreateSubMatrix(change, sub_schurs->is_subs[i], sub_schurs->is_subs[i], MAT_INITIAL_MATRIX, &change_subt));
698: PetscCall(MatConvert(change_subt, MATSEQDENSE, MAT_INITIAL_MATRIX, &change_sub));
699: PetscCall(MatDestroy(&change_subt));
700: }
701: PetscCall(KSPSetOperators(sub_schurs->change[i], change_sub, change_sub));
702: PetscCall(MatDestroy(&change_sub));
703: PetscCall(KSPSetOptionsPrefix(sub_schurs->change[i], sub_schurs->prefix));
704: PetscCall(KSPAppendOptionsPrefix(sub_schurs->change[i], "sub_schurs_change_"));
705: }
706: PetscCall(ISDestroy(&change_primal_all));
707: }
709: /* Local matrix of all local Schur on subsets (transposed) */
710: if (!sub_schurs->S_Ej_all) {
711: Mat T;
712: PetscScalar *v;
713: PetscInt *ii, *jj;
714: PetscInt cum, i, j, k;
716: /* MatSeqAIJSetPreallocation + MatSetValues is slow for these kind of matrices (may have large blocks)
717: Allocate properly a representative matrix and duplicate */
718: PetscCall(PetscMalloc3(local_size + 1, &ii, local_stash_size, &jj, local_stash_size, &v));
719: ii[0] = 0;
720: cum = 0;
721: for (i = 0; i < sub_schurs->n_subs; i++) {
722: PetscCall(ISGetLocalSize(sub_schurs->is_subs[i], &subset_size));
723: for (j = 0; j < subset_size; j++) {
724: const PetscInt row = cum + j;
725: PetscInt col = cum;
727: ii[row + 1] = ii[row] + subset_size;
728: for (k = ii[row]; k < ii[row + 1]; k++) {
729: jj[k] = col;
730: col++;
731: }
732: }
733: cum += subset_size;
734: }
735: PetscCall(MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, local_size, local_size, ii, jj, v, &T));
736: PetscCall(MatDuplicate(T, MAT_DO_NOT_COPY_VALUES, &sub_schurs->S_Ej_all));
737: PetscCall(MatDestroy(&T));
738: PetscCall(PetscFree3(ii, jj, v));
739: }
740: /* matrices for deluxe scaling and adaptive selection */
741: if (compute_Stilda) {
742: if (!sub_schurs->sum_S_Ej_tilda_all) PetscCall(MatDuplicate(sub_schurs->S_Ej_all, MAT_DO_NOT_COPY_VALUES, &sub_schurs->sum_S_Ej_tilda_all));
743: if (!sub_schurs->sum_S_Ej_inv_all && deluxe) PetscCall(MatDuplicate(sub_schurs->S_Ej_all, MAT_DO_NOT_COPY_VALUES, &sub_schurs->sum_S_Ej_inv_all));
744: }
746: /* Compute Schur complements explicitly */
747: F = NULL;
748: if (!sub_schurs->schur_explicit) {
749: /* this code branch is used when MatFactor with Schur complement support is not present or when explicitly requested;
750: it is not efficient, unless the economic version of the scaling is used */
751: Mat S_Ej_expl;
752: PetscScalar *work;
753: PetscInt j, *dummy_idx;
754: PetscBool Sdense;
756: PetscCall(PetscMalloc2(max_subset_size, &dummy_idx, max_subset_size * max_subset_size, &work));
757: local_size = 0;
758: for (i = 0; i < sub_schurs->n_subs; i++) {
759: IS is_subset_B;
760: Mat AE_EE, AE_IE, AE_EI, S_Ej;
762: /* subsets in original and boundary numbering */
763: PetscCall(ISGlobalToLocalMappingApplyIS(sub_schurs->BtoNmap, IS_GTOLM_DROP, sub_schurs->is_subs[i], &is_subset_B));
764: /* EE block */
765: PetscCall(MatCreateSubMatrix(A_BB, is_subset_B, is_subset_B, MAT_INITIAL_MATRIX, &AE_EE));
766: /* IE block */
767: PetscCall(MatCreateSubMatrix(A_IB, is_I, is_subset_B, MAT_INITIAL_MATRIX, &AE_IE));
768: /* EI block */
769: if (sub_schurs->is_symmetric) {
770: PetscCall(MatCreateTranspose(AE_IE, &AE_EI));
771: } else if (sub_schurs->is_hermitian) {
772: PetscCall(MatCreateHermitianTranspose(AE_IE, &AE_EI));
773: } else {
774: PetscCall(MatCreateSubMatrix(A_BI, is_subset_B, is_I, MAT_INITIAL_MATRIX, &AE_EI));
775: }
776: PetscCall(ISDestroy(&is_subset_B));
777: PetscCall(MatCreateSchurComplement(AE_II, AE_II, AE_IE, AE_EI, AE_EE, &S_Ej));
778: PetscCall(MatDestroy(&AE_EE));
779: PetscCall(MatDestroy(&AE_IE));
780: PetscCall(MatDestroy(&AE_EI));
781: if (AE_II == A_II) { /* we can reuse the same ksp */
782: KSP ksp;
783: PetscCall(MatSchurComplementGetKSP(sub_schurs->S, &ksp));
784: PetscCall(MatSchurComplementSetKSP(S_Ej, ksp));
785: } else { /* build new ksp object which inherits ksp and pc types from the original one */
786: KSP origksp, schurksp;
787: PC origpc, schurpc;
788: KSPType ksp_type;
789: PetscInt n_internal;
790: PetscBool ispcnone;
792: PetscCall(MatSchurComplementGetKSP(sub_schurs->S, &origksp));
793: PetscCall(MatSchurComplementGetKSP(S_Ej, &schurksp));
794: PetscCall(KSPGetType(origksp, &ksp_type));
795: PetscCall(KSPSetType(schurksp, ksp_type));
796: PetscCall(KSPGetPC(schurksp, &schurpc));
797: PetscCall(KSPGetPC(origksp, &origpc));
798: PetscCall(PetscObjectTypeCompare((PetscObject)origpc, PCNONE, &ispcnone));
799: if (!ispcnone) {
800: PCType pc_type;
801: PetscCall(PCGetType(origpc, &pc_type));
802: PetscCall(PCSetType(schurpc, pc_type));
803: } else {
804: PetscCall(PCSetType(schurpc, PCLU));
805: }
806: PetscCall(ISGetSize(is_I, &n_internal));
807: if (!n_internal) { /* UMFPACK gives error with 0 sized problems */
808: MatSolverType solver = NULL;
809: PetscCall(PCFactorGetMatSolverType(origpc, &solver));
810: if (solver) PetscCall(PCFactorSetMatSolverType(schurpc, solver));
811: }
812: PetscCall(KSPSetUp(schurksp));
813: }
814: PetscCall(ISGetLocalSize(sub_schurs->is_subs[i], &subset_size));
815: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, subset_size, subset_size, work, &S_Ej_expl));
816: PetscCall(PCBDDCComputeExplicitSchur(S_Ej, sub_schurs->is_symmetric, MAT_REUSE_MATRIX, &S_Ej_expl));
817: PetscCall(PetscObjectTypeCompare((PetscObject)S_Ej_expl, MATSEQDENSE, &Sdense));
818: PetscCheck(Sdense, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not yet implemented for sparse matrices");
819: for (j = 0; j < subset_size; j++) dummy_idx[j] = local_size + j;
820: PetscCall(MatSetValues(sub_schurs->S_Ej_all, subset_size, dummy_idx, subset_size, dummy_idx, work, INSERT_VALUES));
821: PetscCall(MatDestroy(&S_Ej));
822: PetscCall(MatDestroy(&S_Ej_expl));
823: local_size += subset_size;
824: }
825: PetscCall(PetscFree2(dummy_idx, work));
826: /* free */
827: PetscCall(ISDestroy(&is_I));
828: PetscCall(MatDestroy(&AE_II));
829: PetscCall(PetscFree(all_local_idx_N));
830: } else {
831: Mat A, cs_AIB_mat = NULL, benign_AIIm1_ones_mat = NULL;
832: Mat *gdswA;
833: Vec Dall = NULL;
834: IS is_A_all, *is_p_r = NULL, is_schur;
835: MatType Stype;
836: PetscScalar *work, *S_data, *schur_factor, infty = PETSC_MAX_REAL;
837: PetscScalar *SEj_arr = NULL, *SEjinv_arr = NULL;
838: const PetscScalar *rS_data;
839: PetscInt n, n_I, size_schur, size_active_schur, cum, cum2;
840: PetscBool economic, solver_S, S_lower_triangular = PETSC_FALSE;
841: PetscBool schur_has_vertices, factor_workaround;
842: PetscBool use_cholesky;
843: #if PetscDefined(HAVE_VIENNACL) || PetscDefined(HAVE_CUDA)
844: PetscBool oldpin;
845: #endif
846: /* multi-element */
847: IS *is_sub_all = NULL, *is_sub_schur_all = NULL, *is_sub_schur = NULL;
849: /* get sizes */
850: n_I = 0;
851: if (is_I_layer) PetscCall(ISGetLocalSize(is_I_layer, &n_I));
852: economic = PETSC_FALSE;
853: PetscCall(ISGetLocalSize(sub_schurs->is_I, &cum));
854: if (cum != n_I) economic = PETSC_TRUE;
855: PetscCall(MatGetLocalSize(sub_schurs->A, &n, NULL));
856: size_active_schur = local_size;
858: /* import scaling vector (wrong formulation if we have 3D edges) */
859: if (scaling && compute_Stilda) {
860: const PetscScalar *array;
861: PetscScalar *array2;
862: const PetscInt *idxs;
864: PetscCall(ISGetIndices(sub_schurs->is_Ej_all, &idxs));
865: PetscCall(VecCreateSeq(PETSC_COMM_SELF, size_active_schur, &Dall));
866: PetscCall(VecGetArrayRead(scaling, &array));
867: PetscCall(VecGetArray(Dall, &array2));
868: for (PetscInt i = 0; i < size_active_schur; i++) array2[i] = array[idxs[i]];
869: PetscCall(VecRestoreArray(Dall, &array2));
870: PetscCall(VecRestoreArrayRead(scaling, &array));
871: PetscCall(ISRestoreIndices(sub_schurs->is_Ej_all, &idxs));
872: deluxe = PETSC_FALSE;
873: }
875: /* size active schurs does not count any dirichlet or vertex dof on the interface */
876: factor_workaround = PETSC_FALSE;
877: schur_has_vertices = PETSC_FALSE;
878: cum = n_I + size_active_schur;
879: if (sub_schurs->is_dir) {
880: const PetscInt *idxs;
881: PetscInt n_dir;
883: PetscCall(ISGetLocalSize(sub_schurs->is_dir, &n_dir));
884: PetscCall(ISGetIndices(sub_schurs->is_dir, &idxs));
885: PetscCall(PetscArraycpy(all_local_idx_N + cum, idxs, n_dir));
886: if (multi_element)
887: for (PetscInt j = 0; j < n_dir; j++) all_local_subid_N[j + cum] = local_subs[idxs[j]];
888: PetscCall(ISRestoreIndices(sub_schurs->is_dir, &idxs));
889: cum += n_dir;
890: if (!sub_schurs->gdsw) factor_workaround = PETSC_TRUE;
891: }
892: /* include the primal vertices in the Schur complement */
893: if (exact_schur && sub_schurs->is_vertices && (compute_Stilda || benign_n)) {
894: PetscInt n_v;
896: PetscCall(ISGetLocalSize(sub_schurs->is_vertices, &n_v));
897: if (n_v) {
898: const PetscInt *idxs;
900: PetscCall(ISGetIndices(sub_schurs->is_vertices, &idxs));
901: PetscCall(PetscArraycpy(all_local_idx_N + cum, idxs, n_v));
902: if (multi_element)
903: for (PetscInt j = 0; j < n_v; j++) all_local_subid_N[j + cum] = local_subs[idxs[j]];
904: PetscCall(ISRestoreIndices(sub_schurs->is_vertices, &idxs));
905: cum += n_v;
906: if (!sub_schurs->gdsw) factor_workaround = PETSC_TRUE;
907: schur_has_vertices = PETSC_TRUE;
908: }
909: }
910: size_schur = cum - n_I;
911: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, cum, all_local_idx_N, PETSC_OWN_POINTER, &is_A_all));
912: #if PetscDefined(HAVE_VIENNACL) || PetscDefined(HAVE_CUDA)
913: oldpin = sub_schurs->A->boundtocpu;
914: PetscCall(MatBindToCPU(sub_schurs->A, PETSC_TRUE));
915: #endif
916: if (cum == n) {
917: PetscCall(ISSetPermutation(is_A_all));
918: PetscCall(MatPermute(sub_schurs->A, is_A_all, is_A_all, &A));
919: } else {
920: PetscCall(MatCreateSubMatrix(sub_schurs->A, is_A_all, is_A_all, MAT_INITIAL_MATRIX, &A));
921: }
922: #if PetscDefined(HAVE_VIENNACL) || PetscDefined(HAVE_CUDA)
923: PetscCall(MatBindToCPU(sub_schurs->A, oldpin));
924: #endif
925: PetscCall(MatSetOptionsPrefixFactor(A, sub_schurs->prefix));
926: PetscCall(MatAppendOptionsPrefixFactor(A, "sub_schurs_"));
927: /* subsets ordered last */
928: PetscCall(ISCreateStride(PETSC_COMM_SELF, size_schur, n_I, 1, &is_schur));
930: if (multi_element) {
931: PetscInt *idx_sub;
933: PetscCall(PetscMalloc3(n_local_subs, &is_sub_all, n_local_subs, &is_sub_schur_all, n_local_subs, &is_sub_schur));
934: PetscCall(PetscMalloc1(n + size_schur, &idx_sub));
935: for (PetscInt sub = 0; sub < n_local_subs; sub++) {
936: PetscInt size_sub = 0, size_schur_sub = 0, size_I_sub;
938: for (PetscInt j = 0; j < n_I; j++)
939: if (all_local_subid_N[j] == sub) idx_sub[size_sub++] = j;
940: size_I_sub = size_sub;
941: for (PetscInt j = n_I; j < n_I + size_schur; j++)
942: if (all_local_subid_N[j] == sub) {
943: idx_sub[size_sub++] = j;
944: idx_sub[n + size_schur_sub++] = j - n_I;
945: }
947: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size_sub, idx_sub, PETSC_COPY_VALUES, &is_sub_all[sub]));
948: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size_schur_sub, idx_sub + n, PETSC_COPY_VALUES, &is_sub_schur[sub]));
949: PetscCall(ISCreateStride(PETSC_COMM_SELF, size_schur_sub, size_I_sub, 1, &is_sub_schur_all[sub]));
950: }
951: PetscCall(PetscFree(idx_sub));
952: }
954: /* if we actually change the basis for the pressures, LDL^T factors will use a lot of memory
955: this is a workaround */
956: if (benign_n) {
957: Vec v, benign_AIIm1_ones;
958: ISLocalToGlobalMapping N_to_reor;
959: IS is_p0, is_p0_p;
960: PetscScalar *cs_AIB, *AIIm1_data;
961: PetscInt sizeA;
963: PetscCall(ISLocalToGlobalMappingCreateIS(is_A_all, &N_to_reor));
964: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, benign_n, benign_p0_lidx, PETSC_COPY_VALUES, &is_p0));
965: PetscCall(ISGlobalToLocalMappingApplyIS(N_to_reor, IS_GTOLM_DROP, is_p0, &is_p0_p));
966: PetscCall(ISDestroy(&is_p0));
967: PetscCall(MatCreateVecs(A, &v, &benign_AIIm1_ones));
968: PetscCall(VecGetSize(v, &sizeA));
969: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, sizeA, benign_n, NULL, &benign_AIIm1_ones_mat));
970: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, size_schur, benign_n, NULL, &cs_AIB_mat));
971: PetscCall(MatDenseGetArray(cs_AIB_mat, &cs_AIB));
972: PetscCall(MatDenseGetArray(benign_AIIm1_ones_mat, &AIIm1_data));
973: PetscCall(PetscMalloc1(benign_n, &is_p_r));
974: /* compute colsum of A_IB restricted to pressures */
975: for (PetscInt i = 0; i < benign_n; i++) {
976: const PetscScalar *array;
977: const PetscInt *idxs;
978: PetscInt nz;
980: PetscCall(ISGlobalToLocalMappingApplyIS(N_to_reor, IS_GTOLM_DROP, benign_zerodiag_subs[i], &is_p_r[i]));
981: PetscCall(ISGetLocalSize(is_p_r[i], &nz));
982: PetscCall(ISGetIndices(is_p_r[i], &idxs));
983: for (PetscInt j = 0; j < nz; j++) AIIm1_data[idxs[j] + sizeA * i] = 1.;
984: PetscCall(ISRestoreIndices(is_p_r[i], &idxs));
985: PetscCall(VecPlaceArray(benign_AIIm1_ones, AIIm1_data + sizeA * i));
986: PetscCall(MatMult(A, benign_AIIm1_ones, v));
987: PetscCall(VecResetArray(benign_AIIm1_ones));
988: PetscCall(VecGetArrayRead(v, &array));
989: for (PetscInt j = 0; j < size_schur; j++) {
990: #if PetscDefined(USE_COMPLEX)
991: cs_AIB[i * size_schur + j] = (PetscRealPart(array[j + n_I]) / nz + PETSC_i * (PetscImaginaryPart(array[j + n_I]) / nz));
992: #else
993: cs_AIB[i * size_schur + j] = array[j + n_I] / nz;
994: #endif
995: }
996: PetscCall(VecRestoreArrayRead(v, &array));
997: }
998: PetscCall(MatDenseRestoreArray(cs_AIB_mat, &cs_AIB));
999: PetscCall(MatDenseRestoreArray(benign_AIIm1_ones_mat, &AIIm1_data));
1000: PetscCall(VecDestroy(&v));
1001: PetscCall(VecDestroy(&benign_AIIm1_ones));
1002: PetscCall(MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, PETSC_FALSE));
1003: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE));
1004: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
1005: PetscCall(MatZeroRowsColumnsIS(A, is_p0_p, 1.0, NULL, NULL));
1006: PetscCall(ISDestroy(&is_p0_p));
1007: PetscCall(ISLocalToGlobalMappingDestroy(&N_to_reor));
1008: }
1009: PetscCall(MatSetOption(A, MAT_SYMMETRIC, sub_schurs->is_symmetric));
1010: PetscCall(MatSetOption(A, MAT_HERMITIAN, sub_schurs->is_hermitian));
1011: PetscCall(MatSetOption(A, MAT_SPD, sub_schurs->is_posdef));
1013: /* for complexes, symmetric and hermitian at the same time implies null imaginary part */
1014: use_cholesky = (PetscBool)((use_potr || use_sytr) && sub_schurs->is_hermitian && sub_schurs->is_symmetric);
1016: /* when using the benign subspace trick, the local Schur complements are SPD */
1017: /* MKL_PARDISO does not handle well the computation of a Schur complement from a symmetric indefinite factorization
1018: Use LU and adapt pivoting perturbation (still, solution is not as accurate as with using MUMPS) */
1019: if (benign_trick) {
1020: sub_schurs->is_posdef = PETSC_TRUE;
1021: PetscCall(PetscStrcmp(sub_schurs->mat_solver_type, MATSOLVERMKL_PARDISO, &flg));
1022: if (flg) use_cholesky = PETSC_FALSE;
1023: }
1024: if (sub_schurs->mat_factor_type == MAT_FACTOR_NONE) sub_schurs->mat_factor_type = use_cholesky ? MAT_FACTOR_CHOLESKY : MAT_FACTOR_LU;
1026: if (n_I && !multi_element) {
1027: char stype[64];
1028: PetscBool gpu = PETSC_FALSE;
1030: PetscCall(MatGetFactor(A, sub_schurs->mat_solver_type, sub_schurs->mat_factor_type, &F));
1031: PetscCheck(F, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MatGetFactor not supported by matrix instance of type %s. Rerun with \"-info :mat | grep MatGetFactor_\" for additional information", ((PetscObject)A)->type_name);
1032: PetscCall(MatSetErrorIfFailure(A, PETSC_TRUE));
1033: #if PetscDefined(HAVE_MKL_PARDISO)
1034: if (benign_trick) PetscCall(MatMkl_PardisoSetCntl(F, 10, 10));
1035: #endif
1036: PetscCall(MatFactorSetSchurIS(F, is_schur));
1038: /* factorization step */
1039: switch (sub_schurs->mat_factor_type) {
1040: case MAT_FACTOR_CHOLESKY:
1041: PetscCall(MatCholeskyFactorSymbolic(F, A, NULL, NULL));
1042: /* be sure that icntl 19 is not set by command line */
1043: PetscCall(MatMumpsSetIcntl(F, 19, 2));
1044: PetscCall(MatCholeskyFactorNumeric(F, A, NULL));
1045: S_lower_triangular = PETSC_TRUE;
1046: break;
1047: case MAT_FACTOR_LU:
1048: PetscCall(MatLUFactorSymbolic(F, A, NULL, NULL, NULL));
1049: /* be sure that icntl 19 is not set by command line */
1050: PetscCall(MatMumpsSetIcntl(F, 19, 3));
1051: PetscCall(MatLUFactorNumeric(F, A, NULL));
1052: break;
1053: default:
1054: SETERRQ(PetscObjectComm((PetscObject)F), PETSC_ERR_SUP, "Unsupported factor type %s", MatFactorTypes[sub_schurs->mat_factor_type]);
1055: }
1056: PetscCall(MatViewFromOptions(F, (PetscObject)A, "-mat_factor_view"));
1058: if (matl_dbg_viewer) {
1059: Mat S;
1060: IS is;
1062: PetscCall(PetscObjectSetName((PetscObject)A, "A"));
1063: PetscCall(MatView(A, matl_dbg_viewer));
1064: PetscCall(MatFactorCreateSchurComplement(F, &S, NULL));
1065: PetscCall(PetscObjectSetName((PetscObject)S, "S"));
1066: PetscCall(MatView(S, matl_dbg_viewer));
1067: PetscCall(MatDestroy(&S));
1068: PetscCall(ISCreateStride(PETSC_COMM_SELF, n_I, 0, 1, &is));
1069: PetscCall(PetscObjectSetName((PetscObject)is, "I"));
1070: PetscCall(ISView(is, matl_dbg_viewer));
1071: PetscCall(ISDestroy(&is));
1072: PetscCall(ISCreateStride(PETSC_COMM_SELF, size_schur, n_I, 1, &is));
1073: PetscCall(PetscObjectSetName((PetscObject)is, "B"));
1074: PetscCall(ISView(is, matl_dbg_viewer));
1075: PetscCall(ISDestroy(&is));
1076: PetscCall(PetscObjectSetName((PetscObject)is_A_all, "IA"));
1077: PetscCall(ISView(is_A_all, matl_dbg_viewer));
1078: for (i = 0, cum = 0; i < sub_schurs->n_subs; i++) {
1079: IS is;
1080: char name[16];
1082: PetscCall(PetscSNPrintf(name, sizeof(name), "IE%" PetscInt_FMT, i));
1083: PetscCall(ISGetLocalSize(sub_schurs->is_subs[i], &subset_size));
1084: PetscCall(ISCreateStride(PETSC_COMM_SELF, subset_size, cum, 1, &is));
1085: PetscCall(PetscObjectSetName((PetscObject)is, name));
1086: PetscCall(ISView(is, matl_dbg_viewer));
1087: PetscCall(ISDestroy(&is));
1088: if (sub_schurs->change) {
1089: Mat T;
1091: PetscCall(PetscSNPrintf(name, sizeof(name), "TE%" PetscInt_FMT, i));
1092: PetscCall(KSPGetOperators(sub_schurs->change[i], &T, NULL));
1093: PetscCall(PetscObjectSetName((PetscObject)T, name));
1094: PetscCall(MatView(T, matl_dbg_viewer));
1095: PetscCall(PetscSNPrintf(name, sizeof(name), "ITE%" PetscInt_FMT, i));
1096: PetscCall(PetscObjectSetName((PetscObject)sub_schurs->change_primal_sub[i], name));
1097: PetscCall(ISView(sub_schurs->change_primal_sub[i], matl_dbg_viewer));
1098: }
1099: cum += subset_size;
1100: }
1101: PetscCall(PetscViewerFlush(matl_dbg_viewer));
1102: }
1104: /* get explicit Schur Complement computed during numeric factorization */
1105: PetscCall(MatFactorGetSchurComplement(F, &S_all, NULL));
1106: PetscCall(PetscStrncpy(stype, MATSEQDENSE, sizeof(stype)));
1107: #if PetscDefined(HAVE_CUDA)
1108: PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &gpu, MATSEQAIJVIENNACL, MATSEQAIJCUSPARSE, ""));
1109: #endif
1110: if (gpu) PetscCall(PetscStrncpy(stype, MATSEQDENSECUDA, sizeof(stype)));
1111: PetscCall(PetscOptionsGetString(NULL, sub_schurs->prefix, "-sub_schurs_schur_mat_type", stype, sizeof(stype), NULL));
1112: PetscCall(MatConvert(S_all, stype, MAT_INPLACE_MATRIX, &S_all));
1113: PetscCall(MatSetOption(S_all, MAT_SPD, sub_schurs->is_posdef));
1114: PetscCall(MatSetOption(S_all, MAT_HERMITIAN, sub_schurs->is_hermitian));
1115: PetscCall(MatGetType(S_all, &Stype));
1117: /* we can reuse the solvers if we are not using the economic version */
1118: reuse_solvers = (PetscBool)(reuse_solvers && !economic && !sub_schurs->graph->multi_element);
1119: if (!sub_schurs->gdsw) {
1120: factor_workaround = (PetscBool)(reuse_solvers && factor_workaround);
1121: if (!sub_schurs->is_posdef && factor_workaround && compute_Stilda && size_active_schur) reuse_solvers = factor_workaround = PETSC_FALSE;
1122: }
1123: solver_S = PETSC_TRUE;
1125: /* update the Schur complement with the change of basis on the pressures */
1126: if (benign_n) {
1127: const PetscScalar *cs_AIB;
1128: PetscScalar *S_data, *AIIm1_data;
1129: Mat S2 = NULL, S3 = NULL; /* dbg */
1130: PetscScalar *S2_data, *S3_data; /* dbg */
1131: Vec v, benign_AIIm1_ones;
1132: PetscInt sizeA;
1134: PetscCall(MatDenseGetArray(S_all, &S_data));
1135: PetscCall(MatCreateVecs(A, &v, &benign_AIIm1_ones));
1136: PetscCall(VecGetSize(v, &sizeA));
1137: PetscCall(MatMumpsSetIcntl(F, 26, 0));
1138: #if PetscDefined(HAVE_MKL_PARDISO)
1139: PetscCall(MatMkl_PardisoSetCntl(F, 70, 1));
1140: #endif
1141: PetscCall(MatDenseGetArrayRead(cs_AIB_mat, &cs_AIB));
1142: PetscCall(MatDenseGetArray(benign_AIIm1_ones_mat, &AIIm1_data));
1143: if (matl_dbg_viewer) {
1144: PetscCall(MatDuplicate(S_all, MAT_DO_NOT_COPY_VALUES, &S2));
1145: PetscCall(MatDuplicate(S_all, MAT_DO_NOT_COPY_VALUES, &S3));
1146: PetscCall(MatDenseGetArray(S2, &S2_data));
1147: PetscCall(MatDenseGetArray(S3, &S3_data));
1148: }
1149: for (i = 0; i < benign_n; i++) {
1150: PetscScalar *array, sum = 0., one = 1., *sums;
1151: const PetscInt *idxs;
1152: PetscInt k, j, nz;
1153: PetscBLASInt B_k, B_n;
1155: PetscCall(PetscCalloc1(benign_n, &sums));
1156: PetscCall(VecPlaceArray(benign_AIIm1_ones, AIIm1_data + sizeA * i));
1157: PetscCall(VecCopy(benign_AIIm1_ones, v));
1158: PetscCall(MatSolve(F, v, benign_AIIm1_ones));
1159: PetscCall(MatMult(A, benign_AIIm1_ones, v));
1160: PetscCall(VecResetArray(benign_AIIm1_ones));
1161: /* p0 dofs (eliminated) are excluded from the sums */
1162: for (k = 0; k < benign_n; k++) {
1163: PetscCall(ISGetLocalSize(is_p_r[k], &nz));
1164: PetscCall(ISGetIndices(is_p_r[k], &idxs));
1165: for (j = 0; j < nz - 1; j++) sums[k] -= AIIm1_data[idxs[j] + sizeA * i];
1166: PetscCall(ISRestoreIndices(is_p_r[k], &idxs));
1167: }
1168: PetscCall(VecGetArrayRead(v, (const PetscScalar **)&array));
1169: if (matl_dbg_viewer) {
1170: Vec vv;
1171: char name[16];
1173: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size_schur, array + n_I, &vv));
1174: PetscCall(PetscSNPrintf(name, sizeof(name), "Pvs%" PetscInt_FMT, i));
1175: PetscCall(PetscObjectSetName((PetscObject)vv, name));
1176: PetscCall(VecView(vv, matl_dbg_viewer));
1177: }
1178: /* perform sparse rank updates on symmetric Schur (TODO: move outside of the loop?) */
1179: /* cs_AIB already scaled by 1./nz */
1180: B_k = 1;
1181: PetscCall(PetscBLASIntCast(size_schur, &B_n));
1182: for (k = 0; k < benign_n; k++) {
1183: sum = sums[k];
1185: if (PetscAbsScalar(sum) == 0.0) continue;
1186: if (k == i) {
1187: if (B_n) PetscCallBLAS("BLASsyrk", BLASsyrk_("L", "N", &B_n, &B_k, &sum, cs_AIB + i * size_schur, &B_n, &one, S_data, &B_n));
1188: if (matl_dbg_viewer && B_n) PetscCallBLAS("BLASsyrk", BLASsyrk_("L", "N", &B_n, &B_k, &sum, cs_AIB + i * size_schur, &B_n, &one, S3_data, &B_n));
1189: } else { /* XXX Is it correct to use symmetric rank-2 update with half of the sum? */
1190: sum /= 2.0;
1191: if (B_n) PetscCallBLAS("BLASsyr2k", BLASsyr2k_("L", "N", &B_n, &B_k, &sum, cs_AIB + k * size_schur, &B_n, cs_AIB + i * size_schur, &B_n, &one, S_data, &B_n));
1192: if (matl_dbg_viewer && B_n) PetscCallBLAS("BLASsyr2k", BLASsyr2k_("L", "N", &B_n, &B_k, &sum, cs_AIB + k * size_schur, &B_n, cs_AIB + i * size_schur, &B_n, &one, S3_data, &B_n));
1193: }
1194: }
1195: sum = 1.;
1196: if (B_n) PetscCallBLAS("BLASsyr2k", BLASsyr2k_("L", "N", &B_n, &B_k, &sum, array + n_I, &B_n, cs_AIB + i * size_schur, &B_n, &one, S_data, &B_n));
1197: if (matl_dbg_viewer && B_n) PetscCallBLAS("BLASsyr2k", BLASsyr2k_("L", "N", &B_n, &B_k, &sum, array + n_I, &B_n, cs_AIB + i * size_schur, &B_n, &one, S2_data, &B_n));
1198: PetscCall(VecRestoreArrayRead(v, (const PetscScalar **)&array));
1199: /* set p0 entry of AIIm1_ones to zero */
1200: PetscCall(ISGetLocalSize(is_p_r[i], &nz));
1201: PetscCall(ISGetIndices(is_p_r[i], &idxs));
1202: for (j = 0; j < benign_n; j++) AIIm1_data[idxs[nz - 1] + sizeA * j] = 0.;
1203: PetscCall(ISRestoreIndices(is_p_r[i], &idxs));
1204: PetscCall(PetscFree(sums));
1205: }
1206: PetscCall(VecDestroy(&benign_AIIm1_ones));
1207: if (matl_dbg_viewer) {
1208: PetscCall(MatDenseRestoreArray(S2, &S2_data));
1209: PetscCall(MatDenseRestoreArray(S3, &S3_data));
1210: }
1211: if (!S_lower_triangular) { /* I need to expand the upper triangular data (column-oriented) */
1212: for (PetscInt k = 0; k < size_schur; k++) {
1213: for (PetscInt j = k; j < size_schur; j++) S_data[j * size_schur + k] = PetscConj(S_data[k * size_schur + j]);
1214: }
1215: }
1217: /* restore defaults */
1218: PetscCall(MatMumpsSetIcntl(F, 26, -1));
1219: #if PetscDefined(HAVE_MKL_PARDISO)
1220: PetscCall(MatMkl_PardisoSetCntl(F, 70, 0));
1221: #endif
1222: PetscCall(MatDenseRestoreArrayRead(cs_AIB_mat, &cs_AIB));
1223: PetscCall(MatDenseRestoreArray(benign_AIIm1_ones_mat, &AIIm1_data));
1224: PetscCall(VecDestroy(&v));
1225: PetscCall(MatDenseRestoreArray(S_all, &S_data));
1226: if (matl_dbg_viewer) {
1227: Mat S;
1229: PetscCall(MatFactorRestoreSchurComplement(F, &S_all, MAT_FACTOR_SCHUR_UNFACTORED));
1230: PetscCall(MatFactorCreateSchurComplement(F, &S, NULL));
1231: PetscCall(PetscObjectSetName((PetscObject)S, "Sb"));
1232: PetscCall(MatView(S, matl_dbg_viewer));
1233: PetscCall(MatDestroy(&S));
1234: PetscCall(PetscObjectSetName((PetscObject)S2, "S2P"));
1235: PetscCall(MatView(S2, matl_dbg_viewer));
1236: PetscCall(PetscObjectSetName((PetscObject)S3, "S3P"));
1237: PetscCall(MatView(S3, matl_dbg_viewer));
1238: PetscCall(PetscObjectSetName((PetscObject)cs_AIB_mat, "cs"));
1239: PetscCall(MatView(cs_AIB_mat, matl_dbg_viewer));
1240: PetscCall(MatFactorGetSchurComplement(F, &S_all, NULL));
1241: }
1242: PetscCall(MatDestroy(&S2));
1243: PetscCall(MatDestroy(&S3));
1244: }
1245: if (!reuse_solvers) {
1246: for (i = 0; i < benign_n; i++) PetscCall(ISDestroy(&is_p_r[i]));
1247: PetscCall(PetscFree(is_p_r));
1248: PetscCall(MatDestroy(&cs_AIB_mat));
1249: PetscCall(MatDestroy(&benign_AIIm1_ones_mat));
1250: }
1251: } else if (multi_element) { /* MUMPS does not support sparse Schur complements. Loop over local subs */
1252: PetscInt *nnz;
1253: const PetscInt *idxs;
1254: PetscInt size_schur_sub;
1256: PetscCall(PetscCalloc1(size_schur, &nnz));
1257: for (PetscInt sub = 0; sub < n_local_subs; sub++) {
1258: PetscCall(ISGetLocalSize(is_sub_schur[sub], &size_schur_sub));
1259: PetscCall(ISGetIndices(is_sub_schur[sub], &idxs));
1260: for (PetscInt j = 0; j < size_schur_sub; j++) nnz[idxs[j]] = size_schur_sub;
1261: PetscCall(ISRestoreIndices(is_sub_schur[sub], &idxs));
1262: }
1263: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, size_schur, size_schur, 0, nnz, &S_all));
1264: PetscCall(MatSetOption(S_all, MAT_ROW_ORIENTED, sub_schurs->is_hermitian));
1265: PetscCall(PetscFree(nnz));
1267: for (PetscInt sub = 0; sub < n_local_subs; sub++) {
1268: Mat Asub, Ssub;
1269: const PetscScalar *vals;
1270: PetscInt size_all_sub;
1272: F = NULL;
1273: PetscCall(ISGetLocalSize(is_sub_schur[sub], &size_schur_sub));
1274: PetscCall(ISGetLocalSize(is_sub_all[sub], &size_all_sub));
1275: PetscCall(MatCreateSubMatrix(A, is_sub_all[sub], is_sub_all[sub], MAT_INITIAL_MATRIX, &Asub));
1276: if (size_schur_sub == size_all_sub) {
1277: /* we can't use MatFactor when size_schur == size_of_the_problem */
1278: PetscCall(MatConvert(Asub, MATDENSE, MAT_INITIAL_MATRIX, &Ssub));
1279: } else {
1280: PetscCall(MatGetFactor(Asub, sub_schurs->mat_solver_type, sub_schurs->mat_factor_type, &F));
1281: PetscCheck(F, PetscObjectComm((PetscObject)Asub), PETSC_ERR_SUP, "MatGetFactor not supported by matrix instance of type %s. Rerun with \"-info :mat | grep MatGetFactor_\" for additional information", ((PetscObject)Asub)->type_name);
1282: PetscCall(MatSetErrorIfFailure(Asub, PETSC_TRUE));
1283: #if PetscDefined(HAVE_MKL_PARDISO)
1284: if (benign_trick) PetscCall(MatMkl_PardisoSetCntl(F, 10, 10));
1285: #endif
1286: /* subsets ordered last */
1287: PetscCall(MatFactorSetSchurIS(F, is_sub_schur_all[sub]));
1289: /* factorization step */
1290: switch (sub_schurs->mat_factor_type) {
1291: case MAT_FACTOR_CHOLESKY:
1292: PetscCall(MatCholeskyFactorSymbolic(F, Asub, NULL, NULL));
1293: /* be sure that icntl 19 is not set by command line */
1294: PetscCall(MatMumpsSetIcntl(F, 19, 2));
1295: PetscCall(MatCholeskyFactorNumeric(F, Asub, NULL));
1296: S_lower_triangular = PETSC_TRUE;
1297: break;
1298: case MAT_FACTOR_LU:
1299: PetscCall(MatLUFactorSymbolic(F, Asub, NULL, NULL, NULL));
1300: /* be sure that icntl 19 is not set by command line */
1301: PetscCall(MatMumpsSetIcntl(F, 19, 3));
1302: PetscCall(MatLUFactorNumeric(F, Asub, NULL));
1303: break;
1304: default:
1305: SETERRQ(PetscObjectComm((PetscObject)F), PETSC_ERR_SUP, "Unsupported factor type %s", MatFactorTypes[sub_schurs->mat_factor_type]);
1306: }
1307: PetscCall(MatFactorCreateSchurComplement(F, &Ssub, NULL));
1308: }
1309: PetscCall(MatDestroy(&Asub));
1310: PetscCall(MatDenseGetArrayRead(Ssub, &vals));
1311: PetscCall(ISGetIndices(is_sub_schur[sub], &idxs));
1312: PetscCall(MatSetValues(S_all, size_schur_sub, idxs, size_schur_sub, idxs, vals, INSERT_VALUES));
1313: PetscCall(ISRestoreIndices(is_sub_schur[sub], &idxs));
1314: PetscCall(MatDenseRestoreArrayRead(Ssub, &vals));
1315: PetscCall(MatDestroy(&Ssub));
1316: PetscCall(MatDestroy(&F));
1317: }
1318: PetscCall(MatAssemblyBegin(S_all, MAT_FINAL_ASSEMBLY));
1319: PetscCall(MatAssemblyEnd(S_all, MAT_FINAL_ASSEMBLY));
1320: PetscCall(MatSetOption(S_all, MAT_SPD, sub_schurs->is_posdef));
1321: PetscCall(MatSetOption(S_all, MAT_HERMITIAN, sub_schurs->is_hermitian));
1322: Stype = MATDENSE;
1323: reuse_solvers = PETSC_FALSE;
1324: factor_workaround = PETSC_FALSE;
1325: solver_S = PETSC_FALSE;
1326: } else { /* we can't use MatFactor when size_schur == size_of_the_problem */
1327: PetscCall(MatConvert(A, MATSEQDENSE, MAT_INITIAL_MATRIX, &S_all));
1328: PetscCall(MatGetType(S_all, &Stype));
1329: reuse_solvers = PETSC_FALSE; /* TODO: why we can't reuse the solvers here? */
1330: factor_workaround = PETSC_FALSE;
1331: solver_S = PETSC_FALSE;
1332: }
1334: if (reuse_solvers) {
1335: Mat A_II, pA_II, Afake;
1336: Vec vec1_B;
1337: PCBDDCReuseSolvers msolv_ctx;
1338: PetscInt n_R;
1340: if (sub_schurs->reuse_solver) {
1341: PetscCall(PCBDDCReuseSolversReset(sub_schurs->reuse_solver));
1342: } else {
1343: PetscCall(PetscNew(&sub_schurs->reuse_solver));
1344: }
1345: msolv_ctx = sub_schurs->reuse_solver;
1346: PetscCall(MatSchurComplementGetSubMatrices(sub_schurs->S, &A_II, &pA_II, NULL, NULL, NULL));
1347: PetscCall(PetscObjectReference((PetscObject)F));
1348: msolv_ctx->F = F;
1349: PetscCall(MatCreateVecs(F, &msolv_ctx->sol, NULL));
1350: /* currently PETSc has no support for MatSolve(F,x,x), so cheat and let rhs and sol share the same memory */
1351: {
1352: PetscScalar *array;
1353: PetscInt n;
1355: PetscCall(VecGetLocalSize(msolv_ctx->sol, &n));
1356: PetscCall(VecGetArray(msolv_ctx->sol, &array));
1357: PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)msolv_ctx->sol), 1, n, array, &msolv_ctx->rhs));
1358: PetscCall(VecRestoreArray(msolv_ctx->sol, &array));
1359: }
1360: msolv_ctx->has_vertices = schur_has_vertices;
1362: /* interior solver */
1363: PetscCall(PCCreate(PetscObjectComm((PetscObject)A_II), &msolv_ctx->interior_solver));
1364: PetscCall(PCSetOperators(msolv_ctx->interior_solver, A_II, pA_II));
1365: PetscCall(PCSetType(msolv_ctx->interior_solver, PCSHELL));
1366: PetscCall(PCShellSetName(msolv_ctx->interior_solver, "Interior solver (w/o Schur factorization)"));
1367: PetscCall(PCShellSetContext(msolv_ctx->interior_solver, msolv_ctx));
1368: PetscCall(PCShellSetView(msolv_ctx->interior_solver, PCBDDCReuseSolvers_View));
1369: PetscCall(PCShellSetApply(msolv_ctx->interior_solver, PCBDDCReuseSolvers_Interior));
1370: PetscCall(PCShellSetApplyTranspose(msolv_ctx->interior_solver, PCBDDCReuseSolvers_InteriorTranspose));
1371: if (sub_schurs->gdsw) PetscCall(PCShellSetDestroy(msolv_ctx->interior_solver, PCBDDCReuseSolvers_Destroy));
1373: /* correction solver */
1374: if (!sub_schurs->gdsw) {
1375: PetscCall(PCCreate(PetscObjectComm((PetscObject)A_II), &msolv_ctx->correction_solver));
1376: PetscCall(PCSetType(msolv_ctx->correction_solver, PCSHELL));
1377: PetscCall(PCShellSetName(msolv_ctx->correction_solver, "Correction solver (with Schur factorization)"));
1378: PetscCall(PCShellSetContext(msolv_ctx->correction_solver, msolv_ctx));
1379: PetscCall(PCShellSetView(msolv_ctx->interior_solver, PCBDDCReuseSolvers_View));
1380: PetscCall(PCShellSetApply(msolv_ctx->correction_solver, PCBDDCReuseSolvers_Correction));
1381: PetscCall(PCShellSetApplyTranspose(msolv_ctx->correction_solver, PCBDDCReuseSolvers_CorrectionTranspose));
1383: /* scatter and vecs for Schur complement solver */
1384: PetscCall(MatCreateVecs(S_all, &msolv_ctx->sol_B, &msolv_ctx->rhs_B));
1385: PetscCall(MatCreateVecs(sub_schurs->S, &vec1_B, NULL));
1386: if (!schur_has_vertices) {
1387: PetscCall(ISGlobalToLocalMappingApplyIS(sub_schurs->BtoNmap, IS_GTOLM_DROP, is_A_all, &msolv_ctx->is_B));
1388: PetscCall(VecScatterCreate(vec1_B, msolv_ctx->is_B, msolv_ctx->sol_B, NULL, &msolv_ctx->correction_scatter_B));
1389: PetscCall(PetscObjectReference((PetscObject)is_A_all));
1390: msolv_ctx->is_R = is_A_all;
1391: } else {
1392: IS is_B_all;
1393: const PetscInt *idxs;
1394: PetscInt dual, n_v, n;
1396: PetscCall(ISGetLocalSize(sub_schurs->is_vertices, &n_v));
1397: dual = size_schur - n_v;
1398: PetscCall(ISGetLocalSize(is_A_all, &n));
1399: PetscCall(ISGetIndices(is_A_all, &idxs));
1400: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is_A_all), dual, idxs + n_I, PETSC_COPY_VALUES, &is_B_all));
1401: PetscCall(ISGlobalToLocalMappingApplyIS(sub_schurs->BtoNmap, IS_GTOLM_DROP, is_B_all, &msolv_ctx->is_B));
1402: PetscCall(ISDestroy(&is_B_all));
1403: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)is_A_all), dual, 0, 1, &is_B_all));
1404: PetscCall(VecScatterCreate(vec1_B, msolv_ctx->is_B, msolv_ctx->sol_B, is_B_all, &msolv_ctx->correction_scatter_B));
1405: PetscCall(ISDestroy(&is_B_all));
1406: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is_A_all), n - n_v, idxs, PETSC_COPY_VALUES, &msolv_ctx->is_R));
1407: PetscCall(ISRestoreIndices(is_A_all, &idxs));
1408: }
1409: PetscCall(ISGetLocalSize(msolv_ctx->is_R, &n_R));
1410: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, n_R, n_R, 0, NULL, &Afake));
1411: PetscCall(MatAssemblyBegin(Afake, MAT_FINAL_ASSEMBLY));
1412: PetscCall(MatAssemblyEnd(Afake, MAT_FINAL_ASSEMBLY));
1413: PetscCall(PCSetOperators(msolv_ctx->correction_solver, Afake, Afake));
1414: PetscCall(MatDestroy(&Afake));
1415: PetscCall(VecDestroy(&vec1_B));
1416: }
1417: /* communicate benign info to solver context */
1418: if (benign_n) {
1419: PetscScalar *array;
1421: msolv_ctx->benign_n = benign_n;
1422: msolv_ctx->benign_zerodiag_subs = is_p_r;
1423: PetscCall(PetscMalloc1(benign_n, &msolv_ctx->benign_save_vals));
1424: msolv_ctx->benign_csAIB = cs_AIB_mat;
1425: PetscCall(MatCreateVecs(cs_AIB_mat, &msolv_ctx->benign_corr_work, NULL));
1426: PetscCall(VecGetArray(msolv_ctx->benign_corr_work, &array));
1427: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size_schur, array, &msolv_ctx->benign_dummy_schur_vec));
1428: PetscCall(VecRestoreArray(msolv_ctx->benign_corr_work, &array));
1429: msolv_ctx->benign_AIIm1ones = benign_AIIm1_ones_mat;
1430: }
1431: } else {
1432: if (sub_schurs->reuse_solver) PetscCall(PCBDDCReuseSolversReset(sub_schurs->reuse_solver));
1433: PetscCall(PetscFree(sub_schurs->reuse_solver));
1434: }
1435: PetscCall(MatDestroy(&A));
1436: PetscCall(ISDestroy(&is_A_all));
1438: /* Work arrays */
1439: PetscCall(PetscMalloc1(max_subset_size * max_subset_size, &work));
1441: /* S_Ej_all */
1442: PetscInt *idx_work = NULL;
1443: cum = cum2 = 0;
1444: if (!multi_element) PetscCall(MatDenseGetArrayRead(S_all, &rS_data));
1445: else PetscCall(PetscMalloc1(max_subset_size, &idx_work));
1446: PetscCall(MatSeqAIJGetArray(sub_schurs->S_Ej_all, &SEj_arr));
1447: if (sub_schurs->sum_S_Ej_inv_all) PetscCall(MatSeqAIJGetArray(sub_schurs->sum_S_Ej_inv_all, &SEjinv_arr));
1448: if (sub_schurs->gdsw) PetscCall(MatCreateSubMatrices(sub_schurs->A, sub_schurs->n_subs, sub_schurs->is_subs, sub_schurs->is_subs, MAT_INITIAL_MATRIX, &gdswA));
1449: for (i = 0; i < sub_schurs->n_subs; i++) {
1450: /* get S_E (or K^i_EE for GDSW) */
1451: PetscCall(ISGetLocalSize(sub_schurs->is_subs[i], &subset_size));
1452: if (sub_schurs->gdsw) {
1453: Mat T;
1455: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, subset_size, subset_size, work, &T));
1456: PetscCall(MatConvert(gdswA[i], MATDENSE, MAT_REUSE_MATRIX, &T));
1457: PetscCall(MatDestroy(&T));
1458: } else {
1459: if (multi_element) { /* transpose copy to workspace */
1460: // XXX CSR directly?
1461: for (PetscInt j = 0; j < subset_size; j++) idx_work[j] = cum + j;
1462: PetscCall(MatGetValues(S_all, subset_size, idx_work, subset_size, idx_work, work));
1463: if (!sub_schurs->is_hermitian) {
1464: for (PetscInt k = 0; k < subset_size; k++) {
1465: for (PetscInt j = k; j < subset_size; j++) {
1466: PetscScalar t = work[k * subset_size + j];
1467: work[k * subset_size + j] = work[j * subset_size + k];
1468: work[j * subset_size + k] = t;
1469: }
1470: }
1471: }
1472: } else if (S_lower_triangular) { /* I need to expand the upper triangular data (column-oriented) */
1473: for (PetscInt k = 0; k < subset_size; k++) {
1474: for (PetscInt j = k; j < subset_size; j++) {
1475: work[k * subset_size + j] = rS_data[cum2 + k * size_schur + j];
1476: work[j * subset_size + k] = PetscConj(rS_data[cum2 + k * size_schur + j]);
1477: }
1478: }
1479: } else { /* just copy to workspace */
1480: for (PetscInt k = 0; k < subset_size; k++) {
1481: for (PetscInt j = 0; j < subset_size; j++) work[k * subset_size + j] = rS_data[cum2 + k * size_schur + j];
1482: }
1483: }
1484: }
1485: /* insert S_E values */
1486: if (sub_schurs->change) {
1487: Mat change_sub, SEj, T;
1489: /* change basis */
1490: PetscCall(KSPGetOperators(sub_schurs->change[i], &change_sub, NULL));
1491: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, subset_size, subset_size, work, &SEj));
1492: if (!sub_schurs->change_with_qr) { /* currently there's no support for PtAP with P SeqAIJ */
1493: Mat T2;
1494: PetscCall(MatTransposeMatMult(change_sub, SEj, MAT_INITIAL_MATRIX, 1.0, &T2));
1495: PetscCall(MatMatMult(T2, change_sub, MAT_INITIAL_MATRIX, 1.0, &T));
1496: PetscCall(MatConvert(T, MATSEQDENSE, MAT_INPLACE_MATRIX, &T));
1497: PetscCall(MatDestroy(&T2));
1498: } else {
1499: PetscCall(MatPtAP(SEj, change_sub, MAT_INITIAL_MATRIX, 1.0, &T));
1500: }
1501: PetscCall(MatCopy(T, SEj, SAME_NONZERO_PATTERN));
1502: PetscCall(MatDestroy(&T));
1503: PetscCall(MatZeroRowsColumnsIS(SEj, sub_schurs->change_primal_sub[i], 1.0, NULL, NULL));
1504: PetscCall(MatDestroy(&SEj));
1505: }
1506: PetscCall(PetscArraycpy(SEj_arr, work, subset_size * subset_size));
1507: if (compute_Stilda) {
1508: if (deluxe) { /* if adaptivity is requested, invert S_E blocks */
1509: Mat M;
1510: const PetscScalar *vals;
1511: PetscBool isdense, isdensecuda;
1513: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, subset_size, subset_size, work, &M));
1514: PetscCall(MatSetOption(M, MAT_SPD, sub_schurs->is_posdef));
1515: PetscCall(MatSetOption(M, MAT_HERMITIAN, sub_schurs->is_hermitian));
1516: if (!PetscBTLookup(sub_schurs->is_edge, i)) PetscCall(MatSetType(M, Stype));
1517: PetscCall(PetscObjectTypeCompare((PetscObject)M, MATSEQDENSE, &isdense));
1518: PetscCall(PetscObjectTypeCompare((PetscObject)M, MATSEQDENSECUDA, &isdensecuda));
1519: switch (sub_schurs->mat_factor_type) {
1520: case MAT_FACTOR_CHOLESKY:
1521: PetscCall(MatCholeskyFactor(M, NULL, NULL));
1522: break;
1523: case MAT_FACTOR_LU:
1524: PetscCall(MatLUFactor(M, NULL, NULL, NULL));
1525: break;
1526: default:
1527: SETERRQ(PetscObjectComm((PetscObject)F), PETSC_ERR_SUP, "Unsupported factor type %s", MatFactorTypes[sub_schurs->mat_factor_type]);
1528: }
1529: if (isdense) {
1530: PetscCall(MatSeqDenseInvertFactors_Private(M));
1531: #if PetscDefined(HAVE_CUDA)
1532: } else if (isdensecuda) {
1533: PetscCall(MatSeqDenseCUDAInvertFactors_Internal(M));
1534: #endif
1535: } else SETERRQ(PetscObjectComm((PetscObject)M), PETSC_ERR_SUP, "Not implemented for type %s", Stype);
1536: PetscCall(MatDenseGetArrayRead(M, &vals));
1537: PetscCall(PetscArraycpy(SEjinv_arr, vals, subset_size * subset_size));
1538: PetscCall(MatDenseRestoreArrayRead(M, &vals));
1539: PetscCall(MatDestroy(&M));
1540: } else if (scaling) { /* not using deluxe */
1541: Mat SEj;
1542: Vec D;
1543: PetscScalar *array;
1545: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, subset_size, subset_size, work, &SEj));
1546: PetscCall(VecGetArray(Dall, &array));
1547: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, subset_size, array + cum, &D));
1548: PetscCall(VecRestoreArray(Dall, &array));
1549: PetscCall(VecShift(D, -1.));
1550: PetscCall(MatDiagonalScale(SEj, D, D));
1551: PetscCall(MatDestroy(&SEj));
1552: PetscCall(VecDestroy(&D));
1553: PetscCall(PetscArraycpy(SEj_arr, work, subset_size * subset_size));
1554: }
1555: }
1556: cum += subset_size;
1557: cum2 += subset_size * (size_schur + 1);
1558: SEj_arr += subset_size * subset_size;
1559: if (SEjinv_arr) SEjinv_arr += subset_size * subset_size;
1560: }
1561: if (sub_schurs->gdsw) PetscCall(MatDestroySubMatrices(sub_schurs->n_subs, &gdswA));
1562: if (!multi_element) PetscCall(MatDenseRestoreArrayRead(S_all, &rS_data));
1563: PetscCall(MatSeqAIJRestoreArray(sub_schurs->S_Ej_all, &SEj_arr));
1564: if (sub_schurs->sum_S_Ej_inv_all) PetscCall(MatSeqAIJRestoreArray(sub_schurs->sum_S_Ej_inv_all, &SEjinv_arr));
1565: if (solver_S) PetscCall(MatFactorRestoreSchurComplement(F, &S_all, MAT_FACTOR_SCHUR_UNFACTORED));
1567: /* may prevent from unneeded copies, since MUMPS or MKL_Pardiso always use CPU memory
1568: however, preliminary tests indicate using GPUs is still faster in the solve phase */
1569: #if PetscDefined(HAVE_VIENNACL) || PetscDefined(HAVE_CUDA)
1570: if (reuse_solvers) {
1571: Mat St;
1572: MatFactorSchurStatus st;
1574: flg = PETSC_FALSE;
1575: PetscCall(PetscOptionsGetBool(NULL, sub_schurs->prefix, "-sub_schurs_schur_pin_to_cpu", &flg, NULL));
1576: PetscCall(MatFactorGetSchurComplement(F, &St, &st));
1577: PetscCall(MatBindToCPU(St, flg));
1578: PetscCall(MatFactorRestoreSchurComplement(F, &St, st));
1579: }
1580: #endif
1582: schur_factor = NULL;
1583: if (compute_Stilda && size_active_schur) {
1584: if (sub_schurs->n_subs == 1 && size_schur == size_active_schur && deluxe) { /* we already computed the inverse */
1585: PetscCall(MatSeqAIJGetArrayWrite(sub_schurs->sum_S_Ej_tilda_all, &SEjinv_arr));
1586: PetscCall(PetscArraycpy(SEjinv_arr, work, size_schur * size_schur));
1587: PetscCall(MatSeqAIJRestoreArrayWrite(sub_schurs->sum_S_Ej_tilda_all, &SEjinv_arr));
1588: } else {
1589: Mat S_all_inv = NULL;
1591: if (solver_S && !sub_schurs->gdsw) {
1592: /* for adaptive selection we need S^-1; for solver reusage we need S_\Delta\Delta^-1.
1593: The latter is not the principal subminor for S^-1. However, the factors can be reused since S_\Delta\Delta is the leading principal submatrix of S */
1594: if (factor_workaround) { /* invert without calling MatFactorInvertSchurComplement, since we are hacking */
1595: PetscScalar *data;
1596: PetscInt nd = 0;
1598: PetscCheck(use_potr, PETSC_COMM_SELF, PETSC_ERR_SUP, "Factor update not yet implemented for non SPD matrices");
1599: PetscCall(MatFactorGetSchurComplement(F, &S_all_inv, NULL));
1600: PetscCall(MatDenseGetArray(S_all_inv, &data));
1601: if (sub_schurs->is_dir) { /* dirichlet dofs could have different scalings */
1602: PetscCall(ISGetLocalSize(sub_schurs->is_dir, &nd));
1603: }
1605: /* factor and invert activedofs and vertices (dirichlet dofs does not contribute) */
1606: if (schur_has_vertices) {
1607: Mat M;
1608: PetscScalar *tdata;
1609: PetscInt nv = 0, news;
1611: PetscCall(ISGetLocalSize(sub_schurs->is_vertices, &nv));
1612: news = size_active_schur + nv;
1613: PetscCall(PetscCalloc1(news * news, &tdata));
1614: for (i = 0; i < size_active_schur; i++) {
1615: PetscCall(PetscArraycpy(tdata + i * (news + 1), data + i * (size_schur + 1), size_active_schur - i));
1616: PetscCall(PetscArraycpy(tdata + i * (news + 1) + size_active_schur - i, data + i * size_schur + size_active_schur + nd, nv));
1617: }
1618: for (i = 0; i < nv; i++) {
1619: PetscInt k = i + size_active_schur;
1620: PetscCall(PetscArraycpy(tdata + k * (news + 1), data + (k + nd) * (size_schur + 1), nv - i));
1621: }
1623: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, news, news, tdata, &M));
1624: PetscCall(MatSetOption(M, MAT_SPD, PETSC_TRUE));
1625: PetscCall(MatCholeskyFactor(M, NULL, NULL));
1626: /* save the factors */
1627: cum = 0;
1628: PetscCall(PetscMalloc1((size_active_schur * (size_active_schur + 1)) / 2 + nd, &schur_factor));
1629: for (i = 0; i < size_active_schur; i++) {
1630: PetscCall(PetscArraycpy(schur_factor + cum, tdata + i * (news + 1), size_active_schur - i));
1631: cum += size_active_schur - i;
1632: }
1633: for (i = 0; i < nd; i++) schur_factor[cum + i] = PetscSqrtReal(PetscRealPart(data[(i + size_active_schur) * (size_schur + 1)]));
1634: S_all_inv->ops->solve = M->ops->solve;
1635: S_all_inv->ops->matsolve = M->ops->matsolve;
1636: S_all_inv->ops->solvetranspose = M->ops->solvetranspose;
1637: S_all_inv->ops->matsolvetranspose = M->ops->matsolvetranspose;
1638: S_all_inv->factortype = MAT_FACTOR_CHOLESKY;
1639: PetscCall(MatSeqDenseInvertFactors_Private(M));
1640: /* move back just the active dofs to the Schur complement */
1641: for (i = 0; i < size_active_schur; i++) PetscCall(PetscArraycpy(data + i * size_schur, tdata + i * news, size_active_schur));
1642: PetscCall(PetscFree(tdata));
1643: PetscCall(MatDestroy(&M));
1644: } else { /* we can factorize and invert just the activedofs */
1645: Mat M;
1646: PetscScalar *aux;
1648: PetscCall(PetscMalloc1(nd, &aux));
1649: for (i = 0; i < nd; i++) aux[i] = 1.0 / data[(i + size_active_schur) * (size_schur + 1)];
1650: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, size_active_schur, size_active_schur, data, &M));
1651: PetscCall(MatDenseSetLDA(M, size_schur));
1652: PetscCall(MatSetOption(M, MAT_SPD, PETSC_TRUE));
1653: PetscCall(MatCholeskyFactor(M, NULL, NULL));
1654: PetscCall(MatSeqDenseInvertFactors_Private(M));
1655: PetscCall(MatDestroy(&M));
1656: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, size_schur, nd, data + size_active_schur * size_schur, &M));
1657: PetscCall(MatZeroEntries(M));
1658: PetscCall(MatDestroy(&M));
1659: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, nd, size_schur, data + size_active_schur, &M));
1660: PetscCall(MatDenseSetLDA(M, size_schur));
1661: PetscCall(MatZeroEntries(M));
1662: PetscCall(MatDestroy(&M));
1663: for (i = 0; i < nd; i++) data[(i + size_active_schur) * (size_schur + 1)] = aux[i];
1664: PetscCall(PetscFree(aux));
1665: }
1666: PetscCall(MatDenseRestoreArray(S_all_inv, &data));
1667: } else { /* use MatFactor calls to invert S */
1668: PetscCall(MatFactorInvertSchurComplement(F));
1669: PetscCall(MatFactorGetSchurComplement(F, &S_all_inv, NULL));
1670: }
1671: } else if (!sub_schurs->gdsw) { /* we need to invert explicitly since we are not using MatFactor for S */
1672: if (multi_element) {
1673: PetscCall(MatDuplicate(S_all, MAT_DO_NOT_COPY_VALUES, &S_all_inv));
1674: PetscCall(MatSetOption(S_all_inv, MAT_ROW_ORIENTED, sub_schurs->is_hermitian));
1675: for (PetscInt sub = 0; sub < n_local_subs; sub++) {
1676: const PetscScalar *vals;
1677: const PetscInt *idxs;
1678: PetscInt size_schur_sub;
1679: Mat M;
1681: PetscCall(MatCreateSubMatrix(S_all, is_sub_schur[sub], is_sub_schur[sub], MAT_INITIAL_MATRIX, &M));
1682: PetscCall(MatConvert(M, MATDENSE, MAT_INPLACE_MATRIX, &M));
1683: PetscCall(MatSetOption(M, MAT_SPD, sub_schurs->is_posdef));
1684: PetscCall(MatSetOption(M, MAT_HERMITIAN, sub_schurs->is_hermitian));
1685: switch (sub_schurs->mat_factor_type) {
1686: case MAT_FACTOR_CHOLESKY:
1687: PetscCall(MatCholeskyFactor(M, NULL, NULL));
1688: break;
1689: case MAT_FACTOR_LU:
1690: PetscCall(MatLUFactor(M, NULL, NULL, NULL));
1691: break;
1692: default:
1693: SETERRQ(PetscObjectComm((PetscObject)F), PETSC_ERR_SUP, "Unsupported factor type %s", MatFactorTypes[sub_schurs->mat_factor_type]);
1694: }
1695: PetscCall(MatSeqDenseInvertFactors_Private(M));
1696: PetscCall(MatDenseGetArrayRead(M, &vals));
1697: PetscCall(ISGetLocalSize(is_sub_schur[sub], &size_schur_sub));
1698: PetscCall(ISGetIndices(is_sub_schur[sub], &idxs));
1699: PetscCall(MatSetValues(S_all_inv, size_schur_sub, idxs, size_schur_sub, idxs, vals, INSERT_VALUES));
1700: PetscCall(ISRestoreIndices(is_sub_schur[sub], &idxs));
1701: PetscCall(MatDenseRestoreArrayRead(M, &vals));
1702: PetscCall(MatDestroy(&M));
1703: }
1704: PetscCall(MatAssemblyBegin(S_all_inv, MAT_FINAL_ASSEMBLY));
1705: PetscCall(MatAssemblyEnd(S_all_inv, MAT_FINAL_ASSEMBLY));
1706: } else {
1707: PetscCall(PetscObjectReference((PetscObject)S_all));
1708: S_all_inv = S_all;
1709: PetscCall(MatDenseGetArray(S_all_inv, &S_data));
1710: PetscCall(PetscBLASIntCast(size_schur, &B_N));
1711: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1712: if (use_potr) {
1713: PetscCallLAPACKInfo("LAPACKpotrf", LAPACKpotrf_("L", &B_N, S_data, &B_N, &info));
1714: PetscCallLAPACKInfo("LAPACKpotri", LAPACKpotri_("L", &B_N, S_data, &B_N, &info));
1715: } else if (use_sytr) {
1716: PetscCallLAPACKInfo("LAPACKsytrf", LAPACKsytrf_("L", &B_N, S_data, &B_N, pivots, Bwork, &B_lwork, &info));
1717: PetscCallLAPACKInfo("LAPACKsytri", LAPACKsytri_("L", &B_N, S_data, &B_N, pivots, Bwork, &info));
1718: } else {
1719: PetscCallLAPACKInfo("LAPACKgetrf", LAPACKgetrf_(&B_N, &B_N, S_data, &B_N, pivots, &info));
1720: PetscCallLAPACKInfo("LAPACKgetri", LAPACKgetri_(&B_N, S_data, &B_N, pivots, Bwork, &B_lwork, &info));
1721: }
1722: PetscCall(PetscLogFlops(1.0 * size_schur * size_schur * size_schur));
1723: PetscCall(PetscFPTrapPop());
1724: PetscCall(MatDenseRestoreArray(S_all_inv, &S_data));
1725: }
1726: } else if (sub_schurs->gdsw) {
1727: Mat tS, tX, SEj, S_II, S_IE, S_EE;
1728: KSP pS_II;
1729: PC pS_II_pc;
1730: IS EE, II;
1731: PetscInt nS;
1733: PetscCall(MatFactorCreateSchurComplement(F, &tS, NULL));
1734: PetscCall(MatGetSize(tS, &nS, NULL));
1735: PetscCall(MatSeqAIJGetArray(sub_schurs->sum_S_Ej_tilda_all, &SEjinv_arr));
1736: for (i = 0, cum = 0; i < sub_schurs->n_subs; i++) { /* naive implementation */
1737: PetscCall(ISGetLocalSize(sub_schurs->is_subs[i], &subset_size));
1738: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, subset_size, subset_size, SEjinv_arr, &SEj));
1740: PetscCall(ISCreateStride(PETSC_COMM_SELF, subset_size, cum, 1, &EE));
1741: PetscCall(ISComplement(EE, 0, nS, &II));
1742: PetscCall(MatCreateSubMatrix(tS, II, II, MAT_INITIAL_MATRIX, &S_II));
1743: PetscCall(MatCreateSubMatrix(tS, II, EE, MAT_INITIAL_MATRIX, &S_IE));
1744: PetscCall(MatCreateSubMatrix(tS, EE, EE, MAT_INITIAL_MATRIX, &S_EE));
1745: PetscCall(ISDestroy(&II));
1746: PetscCall(ISDestroy(&EE));
1748: PetscCall(KSPCreate(PETSC_COMM_SELF, &pS_II));
1749: PetscCall(KSPSetNestLevel(pS_II, 1)); /* do not have direct access to a PC to provide the level of nesting of the KSP */
1750: PetscCall(KSPSetType(pS_II, KSPPREONLY));
1751: PetscCall(KSPGetPC(pS_II, &pS_II_pc));
1752: PetscCall(PCSetType(pS_II_pc, PCSVD));
1753: PetscCall(KSPSetOptionsPrefix(pS_II, sub_schurs->prefix));
1754: PetscCall(KSPAppendOptionsPrefix(pS_II, "pseudo_"));
1755: PetscCall(KSPSetOperators(pS_II, S_II, S_II));
1756: PetscCall(MatDestroy(&S_II));
1757: PetscCall(KSPSetFromOptions(pS_II));
1758: PetscCall(KSPSetUp(pS_II));
1759: PetscCall(MatDuplicate(S_IE, MAT_DO_NOT_COPY_VALUES, &tX));
1760: PetscCall(KSPMatSolve(pS_II, S_IE, tX));
1761: PetscCall(KSPDestroy(&pS_II));
1763: PetscCall(MatTransposeMatMult(S_IE, tX, MAT_REUSE_MATRIX, PETSC_DETERMINE, &SEj));
1764: PetscCall(MatDestroy(&S_IE));
1765: PetscCall(MatDestroy(&tX));
1766: PetscCall(MatAYPX(SEj, -1, S_EE, SAME_NONZERO_PATTERN));
1767: PetscCall(MatDestroy(&S_EE));
1769: PetscCall(MatDestroy(&SEj));
1770: cum += subset_size;
1771: SEjinv_arr += subset_size * subset_size;
1772: }
1773: PetscCall(MatDestroy(&tS));
1774: PetscCall(MatSeqAIJRestoreArray(sub_schurs->sum_S_Ej_tilda_all, &SEjinv_arr));
1775: }
1776: /* S_Ej_tilda_all */
1777: cum = cum2 = 0;
1778: rS_data = NULL;
1779: if (S_all_inv && !multi_element) PetscCall(MatDenseGetArrayRead(S_all_inv, &rS_data));
1780: PetscCall(MatSeqAIJGetArrayWrite(sub_schurs->sum_S_Ej_tilda_all, &SEjinv_arr));
1781: for (i = 0; i < sub_schurs->n_subs; i++) {
1782: PetscInt j;
1784: PetscCall(ISGetLocalSize(sub_schurs->is_subs[i], &subset_size));
1785: /* get (St^-1)_E */
1786: /* Unless we are changing the variables, I don't need to expand to upper triangular since St^-1
1787: will be properly accessed later during adaptive selection */
1788: if (multi_element) { /* transpose copy to workspace */
1789: // XXX CSR directly?
1790: for (PetscInt j = 0; j < subset_size; j++) idx_work[j] = cum + j;
1791: PetscCall(MatGetValues(S_all_inv, subset_size, idx_work, subset_size, idx_work, work));
1792: if (!sub_schurs->is_hermitian) {
1793: for (PetscInt k = 0; k < subset_size; k++) {
1794: for (PetscInt j = k; j < subset_size; j++) {
1795: PetscScalar t = work[k * subset_size + j];
1796: work[k * subset_size + j] = work[j * subset_size + k];
1797: work[j * subset_size + k] = t;
1798: }
1799: }
1800: }
1801: } else if (rS_data) {
1802: if (S_lower_triangular) {
1803: if (sub_schurs->change) {
1804: for (PetscInt k = 0; k < subset_size; k++) {
1805: for (j = k; j < subset_size; j++) {
1806: work[k * subset_size + j] = rS_data[cum2 + k * size_schur + j];
1807: work[j * subset_size + k] = work[k * subset_size + j];
1808: }
1809: }
1810: } else {
1811: for (PetscInt k = 0; k < subset_size; k++) {
1812: for (j = k; j < subset_size; j++) work[k * subset_size + j] = rS_data[cum2 + k * size_schur + j];
1813: }
1814: }
1815: } else {
1816: for (PetscInt k = 0; k < subset_size; k++) {
1817: for (j = 0; j < subset_size; j++) work[k * subset_size + j] = rS_data[cum2 + k * size_schur + j];
1818: }
1819: }
1820: }
1821: if (sub_schurs->change) {
1822: Mat change_sub, SEj, T;
1823: PetscScalar val = sub_schurs->gdsw ? PETSC_SMALL : 1. / PETSC_SMALL;
1825: /* change basis */
1826: PetscCall(KSPGetOperators(sub_schurs->change[i], &change_sub, NULL));
1827: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, subset_size, subset_size, (rS_data || multi_element) ? work : SEjinv_arr, &SEj));
1828: if (!sub_schurs->change_with_qr) { /* currently there's no support for PtAP with P SeqAIJ */
1829: Mat T2;
1830: PetscCall(MatTransposeMatMult(change_sub, SEj, MAT_INITIAL_MATRIX, 1.0, &T2));
1831: PetscCall(MatMatMult(T2, change_sub, MAT_INITIAL_MATRIX, 1.0, &T));
1832: PetscCall(MatDestroy(&T2));
1833: PetscCall(MatConvert(T, MATSEQDENSE, MAT_INPLACE_MATRIX, &T));
1834: } else {
1835: PetscCall(MatPtAP(SEj, change_sub, MAT_INITIAL_MATRIX, 1.0, &T));
1836: }
1837: PetscCall(MatCopy(T, SEj, SAME_NONZERO_PATTERN));
1838: PetscCall(MatDestroy(&T));
1839: PetscCall(MatZeroRowsColumnsIS(SEj, sub_schurs->change_primal_sub[i], val, NULL, NULL));
1840: PetscCall(MatDestroy(&SEj));
1841: }
1842: if (rS_data || multi_element) PetscCall(PetscArraycpy(SEjinv_arr, work, subset_size * subset_size));
1843: cum += subset_size;
1844: cum2 += subset_size * (size_schur + 1);
1845: SEjinv_arr += subset_size * subset_size;
1846: }
1847: PetscCall(MatSeqAIJRestoreArrayWrite(sub_schurs->sum_S_Ej_tilda_all, &SEjinv_arr));
1848: if (S_all_inv) {
1849: if (!multi_element) PetscCall(MatDenseRestoreArrayRead(S_all_inv, &rS_data));
1850: if (solver_S) {
1851: if (schur_has_vertices) {
1852: PetscCall(MatFactorRestoreSchurComplement(F, &S_all_inv, MAT_FACTOR_SCHUR_FACTORED));
1853: } else {
1854: PetscCall(MatFactorRestoreSchurComplement(F, &S_all_inv, MAT_FACTOR_SCHUR_INVERTED));
1855: }
1856: }
1857: }
1858: PetscCall(MatDestroy(&S_all_inv));
1859: }
1861: /* move back factors if needed */
1862: if (schur_has_vertices && factor_workaround && !sub_schurs->gdsw) {
1863: Mat S_tmp;
1864: PetscInt nd = 0;
1865: PetscScalar *data;
1867: PetscCheck(use_potr, PETSC_COMM_SELF, PETSC_ERR_SUP, "Factor update not yet implemented for non SPD matrices");
1868: PetscCheck(solver_S, PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1869: PetscCall(MatFactorGetSchurComplement(F, &S_tmp, NULL));
1870: PetscCall(MatDenseGetArray(S_tmp, &data));
1871: PetscCall(PetscArrayzero(data, size_schur * size_schur));
1873: if (S_lower_triangular) {
1874: cum = 0;
1875: for (i = 0; i < size_active_schur; i++) {
1876: PetscCall(PetscArraycpy(data + i * (size_schur + 1), schur_factor + cum, size_active_schur - i));
1877: cum += size_active_schur - i;
1878: }
1879: } else {
1880: PetscCall(PetscArraycpy(data, schur_factor, size_schur * size_schur));
1881: }
1882: if (sub_schurs->is_dir) {
1883: PetscCall(ISGetLocalSize(sub_schurs->is_dir, &nd));
1884: for (i = 0; i < nd; i++) data[(i + size_active_schur) * (size_schur + 1)] = schur_factor[cum + i];
1885: }
1886: /* workaround: since I cannot modify the matrices used inside the solvers for the forward and backward substitutions,
1887: set the diagonal entry of the Schur factor to a very large value */
1888: for (i = size_active_schur + nd; i < size_schur; i++) data[i * (size_schur + 1)] = infty;
1889: PetscCall(MatDenseRestoreArray(S_tmp, &data));
1890: PetscCall(MatFactorRestoreSchurComplement(F, &S_tmp, MAT_FACTOR_SCHUR_FACTORED));
1891: }
1892: } else if (factor_workaround && !sub_schurs->gdsw) { /* we need to eliminate any unneeded coupling */
1893: PetscScalar *data;
1894: PetscInt nd = 0;
1896: if (sub_schurs->is_dir) { /* dirichlet dofs could have different scalings */
1897: PetscCall(ISGetLocalSize(sub_schurs->is_dir, &nd));
1898: }
1899: PetscCall(MatFactorGetSchurComplement(F, &S_all, NULL));
1900: PetscCall(MatDenseGetArray(S_all, &data));
1901: for (i = 0; i < size_active_schur; i++) PetscCall(PetscArrayzero(data + i * size_schur + size_active_schur, size_schur - size_active_schur));
1902: for (i = size_active_schur + nd; i < size_schur; i++) {
1903: PetscCall(PetscArrayzero(data + i * size_schur + size_active_schur, size_schur - size_active_schur));
1904: data[i * (size_schur + 1)] = infty;
1905: }
1906: PetscCall(MatDenseRestoreArray(S_all, &data));
1907: PetscCall(MatFactorRestoreSchurComplement(F, &S_all, MAT_FACTOR_SCHUR_UNFACTORED));
1908: }
1909: PetscCall(PetscFree(idx_work));
1910: PetscCall(PetscFree(work));
1911: PetscCall(PetscFree(schur_factor));
1912: PetscCall(VecDestroy(&Dall));
1913: PetscCall(ISDestroy(&is_schur));
1914: if (multi_element) {
1915: for (PetscInt sub = 0; sub < n_local_subs; sub++) {
1916: PetscCall(ISDestroy(&is_sub_all[sub]));
1917: PetscCall(ISDestroy(&is_sub_schur_all[sub]));
1918: PetscCall(ISDestroy(&is_sub_schur[sub]));
1919: }
1920: PetscCall(PetscFree3(is_sub_all, is_sub_schur_all, is_sub_schur));
1921: }
1922: }
1923: PetscCall(ISDestroy(&is_I_layer));
1924: PetscCall(MatDestroy(&S_all));
1925: PetscCall(MatDestroy(&A_BB));
1926: PetscCall(MatDestroy(&A_IB));
1927: PetscCall(MatDestroy(&A_BI));
1928: PetscCall(MatDestroy(&F));
1930: PetscCall(MatAssemblyBegin(sub_schurs->S_Ej_all, MAT_FINAL_ASSEMBLY));
1931: PetscCall(MatAssemblyEnd(sub_schurs->S_Ej_all, MAT_FINAL_ASSEMBLY));
1932: if (compute_Stilda) {
1933: PetscCall(MatAssemblyBegin(sub_schurs->sum_S_Ej_tilda_all, MAT_FINAL_ASSEMBLY));
1934: PetscCall(MatAssemblyEnd(sub_schurs->sum_S_Ej_tilda_all, MAT_FINAL_ASSEMBLY));
1935: if (deluxe) {
1936: PetscCall(MatAssemblyBegin(sub_schurs->sum_S_Ej_inv_all, MAT_FINAL_ASSEMBLY));
1937: PetscCall(MatAssemblyEnd(sub_schurs->sum_S_Ej_inv_all, MAT_FINAL_ASSEMBLY));
1938: }
1939: }
1941: /* Get local part of (\sum_j S_Ej) */
1942: if (!sub_schurs->sum_S_Ej_all) PetscCall(MatDuplicate(sub_schurs->S_Ej_all, MAT_DO_NOT_COPY_VALUES, &sub_schurs->sum_S_Ej_all));
1943: PetscCall(VecSet(gstash, 0.0));
1944: PetscCall(MatSeqAIJGetArray(sub_schurs->S_Ej_all, &stasharray));
1945: PetscCall(VecPlaceArray(lstash, stasharray));
1946: PetscCall(VecScatterBegin(sstash, lstash, gstash, ADD_VALUES, SCATTER_FORWARD));
1947: PetscCall(VecScatterEnd(sstash, lstash, gstash, ADD_VALUES, SCATTER_FORWARD));
1948: PetscCall(MatSeqAIJRestoreArray(sub_schurs->S_Ej_all, &stasharray));
1949: PetscCall(VecResetArray(lstash));
1950: PetscCall(MatSeqAIJGetArray(sub_schurs->sum_S_Ej_all, &stasharray));
1951: PetscCall(VecPlaceArray(lstash, stasharray));
1952: PetscCall(VecScatterBegin(sstash, gstash, lstash, INSERT_VALUES, SCATTER_REVERSE));
1953: PetscCall(VecScatterEnd(sstash, gstash, lstash, INSERT_VALUES, SCATTER_REVERSE));
1954: PetscCall(MatSeqAIJRestoreArray(sub_schurs->sum_S_Ej_all, &stasharray));
1955: PetscCall(VecResetArray(lstash));
1957: /* Get local part of (\sum_j S^-1_Ej) (\sum_j St^-1_Ej) */
1958: if (compute_Stilda) {
1959: PetscCall(VecSet(gstash, 0.0));
1960: PetscCall(MatSeqAIJGetArray(sub_schurs->sum_S_Ej_tilda_all, &stasharray));
1961: PetscCall(VecPlaceArray(lstash, stasharray));
1962: PetscCall(VecScatterBegin(sstash, lstash, gstash, ADD_VALUES, SCATTER_FORWARD));
1963: PetscCall(VecScatterEnd(sstash, lstash, gstash, ADD_VALUES, SCATTER_FORWARD));
1964: PetscCall(VecScatterBegin(sstash, gstash, lstash, INSERT_VALUES, SCATTER_REVERSE));
1965: PetscCall(VecScatterEnd(sstash, gstash, lstash, INSERT_VALUES, SCATTER_REVERSE));
1966: PetscCall(MatSeqAIJRestoreArray(sub_schurs->sum_S_Ej_tilda_all, &stasharray));
1967: PetscCall(VecResetArray(lstash));
1968: if (deluxe) {
1969: PetscCall(VecSet(gstash, 0.0));
1970: PetscCall(MatSeqAIJGetArray(sub_schurs->sum_S_Ej_inv_all, &stasharray));
1971: PetscCall(VecPlaceArray(lstash, stasharray));
1972: PetscCall(VecScatterBegin(sstash, lstash, gstash, ADD_VALUES, SCATTER_FORWARD));
1973: PetscCall(VecScatterEnd(sstash, lstash, gstash, ADD_VALUES, SCATTER_FORWARD));
1974: PetscCall(VecScatterBegin(sstash, gstash, lstash, INSERT_VALUES, SCATTER_REVERSE));
1975: PetscCall(VecScatterEnd(sstash, gstash, lstash, INSERT_VALUES, SCATTER_REVERSE));
1976: PetscCall(MatSeqAIJRestoreArray(sub_schurs->sum_S_Ej_inv_all, &stasharray));
1977: PetscCall(VecResetArray(lstash));
1978: } else if (!sub_schurs->gdsw) {
1979: PetscScalar *array;
1980: PetscInt cum;
1982: PetscCall(MatSeqAIJGetArray(sub_schurs->sum_S_Ej_tilda_all, &array));
1983: cum = 0;
1984: for (i = 0; i < sub_schurs->n_subs; i++) {
1985: PetscCall(ISGetLocalSize(sub_schurs->is_subs[i], &subset_size));
1986: PetscCall(PetscBLASIntCast(subset_size, &B_N));
1987: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1988: if (use_potr) {
1989: PetscCallLAPACKInfo("LAPACKpotrf", LAPACKpotrf_("L", &B_N, array + cum, &B_N, &info));
1990: PetscCallLAPACKInfo("LAPACKpotri", LAPACKpotri_("L", &B_N, array + cum, &B_N, &info));
1991: } else if (use_sytr) {
1992: PetscCallLAPACKInfo("LAPACKsytrf", LAPACKsytrf_("L", &B_N, array + cum, &B_N, pivots, Bwork, &B_lwork, &info));
1993: PetscCallLAPACKInfo("LAPACKsytri", LAPACKsytri_("L", &B_N, array + cum, &B_N, pivots, Bwork, &info));
1994: } else {
1995: PetscCallLAPACKInfo("LAPACKgetrf", LAPACKgetrf_(&B_N, &B_N, array + cum, &B_N, pivots, &info));
1996: PetscCallLAPACKInfo("LAPACKgetri", LAPACKgetri_(&B_N, array + cum, &B_N, pivots, Bwork, &B_lwork, &info));
1997: }
1998: PetscCall(PetscLogFlops(1.0 * subset_size * subset_size * subset_size));
1999: PetscCall(PetscFPTrapPop());
2000: cum += subset_size * subset_size;
2001: }
2002: PetscCall(MatSeqAIJRestoreArray(sub_schurs->sum_S_Ej_tilda_all, &array));
2003: PetscCall(PetscObjectReference((PetscObject)sub_schurs->sum_S_Ej_all));
2004: PetscCall(MatDestroy(&sub_schurs->sum_S_Ej_inv_all));
2005: sub_schurs->sum_S_Ej_inv_all = sub_schurs->sum_S_Ej_all;
2006: }
2007: }
2008: PetscCall(VecDestroy(&lstash));
2009: PetscCall(VecDestroy(&gstash));
2010: PetscCall(VecScatterDestroy(&sstash));
2012: if (matl_dbg_viewer) {
2013: if (sub_schurs->S_Ej_all) {
2014: PetscCall(PetscObjectSetName((PetscObject)sub_schurs->S_Ej_all, "SE"));
2015: PetscCall(MatView(sub_schurs->S_Ej_all, matl_dbg_viewer));
2016: }
2017: if (sub_schurs->sum_S_Ej_all) {
2018: PetscCall(PetscObjectSetName((PetscObject)sub_schurs->sum_S_Ej_all, "SSE"));
2019: PetscCall(MatView(sub_schurs->sum_S_Ej_all, matl_dbg_viewer));
2020: }
2021: if (sub_schurs->sum_S_Ej_inv_all) {
2022: PetscCall(PetscObjectSetName((PetscObject)sub_schurs->sum_S_Ej_inv_all, "SSEm"));
2023: PetscCall(MatView(sub_schurs->sum_S_Ej_inv_all, matl_dbg_viewer));
2024: }
2025: if (sub_schurs->sum_S_Ej_tilda_all) {
2026: PetscCall(PetscObjectSetName((PetscObject)sub_schurs->sum_S_Ej_tilda_all, "SSEt"));
2027: PetscCall(MatView(sub_schurs->sum_S_Ej_tilda_all, matl_dbg_viewer));
2028: }
2029: }
2031: /* when not explicit, we need to set the factor type */
2032: if (sub_schurs->mat_factor_type == MAT_FACTOR_NONE) sub_schurs->mat_factor_type = sub_schurs->is_hermitian ? MAT_FACTOR_CHOLESKY : MAT_FACTOR_LU;
2034: /* free workspace */
2035: if (matl_dbg_viewer) PetscCall(PetscViewerFlush(matl_dbg_viewer));
2036: if (sub_schurs->debug) PetscCallMPI(MPI_Barrier(comm_n));
2037: PetscCall(PetscViewerDestroy(&matl_dbg_viewer));
2038: PetscCall(PetscFree2(Bwork, pivots));
2039: PetscCall(PetscCommDestroy(&comm_n));
2040: PetscCall(PetscFree(all_local_subid_N));
2041: PetscFunctionReturn(PETSC_SUCCESS);
2042: }
2044: PetscErrorCode PCBDDCSubSchursInit(PCBDDCSubSchurs sub_schurs, const char *prefix, IS is_I, IS is_B, PCBDDCGraph graph, ISLocalToGlobalMapping BtoNmap, PetscBool copycc, PetscBool gdsw)
2045: {
2046: IS *faces, *edges, *all_cc, vertices;
2047: PetscInt s, i, n_faces, n_edges, n_all_cc;
2048: PetscBool is_sorted, ispardiso, ismumps;
2050: PetscFunctionBegin;
2051: PetscCall(ISSorted(is_I, &is_sorted));
2052: PetscCheck(is_sorted, PetscObjectComm((PetscObject)is_I), PETSC_ERR_PLIB, "IS for I dofs should be shorted");
2053: PetscCall(ISSorted(is_B, &is_sorted));
2054: PetscCheck(is_sorted, PetscObjectComm((PetscObject)is_B), PETSC_ERR_PLIB, "IS for B dofs should be shorted");
2056: /* reset any previous data */
2057: PetscCall(PCBDDCSubSchursReset(sub_schurs));
2059: sub_schurs->gdsw = gdsw;
2060: sub_schurs->graph = graph;
2062: /* get index sets for faces and edges (already sorted by global ordering) */
2063: PetscCall(PCBDDCGraphGetCandidatesIS(graph, &n_faces, &faces, &n_edges, &edges, &vertices));
2064: n_all_cc = n_faces + n_edges;
2065: PetscCall(PetscBTCreate(n_all_cc, &sub_schurs->is_edge));
2066: PetscCall(PetscMalloc1(n_all_cc, &all_cc));
2067: n_all_cc = 0;
2068: for (i = 0; i < n_faces; i++) {
2069: PetscCall(ISGetSize(faces[i], &s));
2070: if (!s) continue;
2071: if (copycc) PetscCall(ISDuplicate(faces[i], &all_cc[n_all_cc]));
2072: else {
2073: PetscCall(PetscObjectReference((PetscObject)faces[i]));
2074: all_cc[n_all_cc] = faces[i];
2075: }
2076: n_all_cc++;
2077: }
2078: for (i = 0; i < n_edges; i++) {
2079: PetscCall(ISGetSize(edges[i], &s));
2080: if (!s) continue;
2081: if (copycc) PetscCall(ISDuplicate(edges[i], &all_cc[n_all_cc]));
2082: else {
2083: PetscCall(PetscObjectReference((PetscObject)edges[i]));
2084: all_cc[n_all_cc] = edges[i];
2085: }
2086: PetscCall(PetscBTSet(sub_schurs->is_edge, n_all_cc));
2087: n_all_cc++;
2088: }
2089: PetscCall(PetscObjectReference((PetscObject)vertices));
2090: sub_schurs->is_vertices = vertices;
2091: PetscCall(PCBDDCGraphRestoreCandidatesIS(graph, &n_faces, &faces, &n_edges, &edges, &vertices));
2092: sub_schurs->is_dir = NULL;
2093: PetscCall(PCBDDCGraphGetDirichletDofsB(graph, &sub_schurs->is_dir));
2095: /* Determine if MatFactor can be used */
2096: PetscCall(PetscStrallocpy(prefix, &sub_schurs->prefix));
2097: #if PetscDefined(HAVE_MUMPS)
2098: PetscCall(PetscStrncpy(sub_schurs->mat_solver_type, MATSOLVERMUMPS, sizeof(sub_schurs->mat_solver_type)));
2099: #elif PetscDefined(HAVE_MKL_PARDISO)
2100: PetscCall(PetscStrncpy(sub_schurs->mat_solver_type, MATSOLVERMKL_PARDISO, sizeof(sub_schurs->mat_solver_type)));
2101: #else
2102: PetscCall(PetscStrncpy(sub_schurs->mat_solver_type, MATSOLVERPETSC, sizeof(sub_schurs->mat_solver_type)));
2103: #endif
2104: sub_schurs->mat_factor_type = MAT_FACTOR_NONE;
2105: sub_schurs->is_hermitian = PetscDefined(USE_COMPLEX) ? PETSC_FALSE : PETSC_TRUE; /* Hermitian Cholesky is not supported by PETSc and external packages */
2106: sub_schurs->is_posdef = PETSC_TRUE;
2107: sub_schurs->is_symmetric = PETSC_TRUE;
2108: sub_schurs->debug = PETSC_FALSE;
2109: sub_schurs->restrict_comm = PETSC_FALSE;
2110: PetscOptionsBegin(PetscObjectComm((PetscObject)graph->l2gmap), sub_schurs->prefix, "BDDC sub_schurs options", "PC");
2111: PetscCall(PetscOptionsString("-sub_schurs_mat_solver_type", "Specific direct solver to use", NULL, sub_schurs->mat_solver_type, sub_schurs->mat_solver_type, sizeof(sub_schurs->mat_solver_type), NULL));
2112: PetscCall(PetscOptionsEnum("-sub_schurs_mat_factor_type", "Factor type to use. Use MAT_FACTOR_NONE for automatic selection", NULL, MatFactorTypes, (PetscEnum)sub_schurs->mat_factor_type, (PetscEnum *)&sub_schurs->mat_factor_type, NULL));
2113: PetscCall(PetscOptionsBool("-sub_schurs_symmetric", "Symmetric problem", NULL, sub_schurs->is_symmetric, &sub_schurs->is_symmetric, NULL));
2114: PetscCall(PetscOptionsBool("-sub_schurs_hermitian", "Hermitian problem", NULL, sub_schurs->is_hermitian, &sub_schurs->is_hermitian, NULL));
2115: PetscCall(PetscOptionsBool("-sub_schurs_posdef", "Positive definite problem", NULL, sub_schurs->is_posdef, &sub_schurs->is_posdef, NULL));
2116: PetscCall(PetscOptionsBool("-sub_schurs_restrictcomm", "Restrict communicator on active processes only", NULL, sub_schurs->restrict_comm, &sub_schurs->restrict_comm, NULL));
2117: PetscCall(PetscOptionsBool("-sub_schurs_debug", "Debug output", NULL, sub_schurs->debug, &sub_schurs->debug, NULL));
2118: PetscOptionsEnd();
2119: PetscCall(PetscStrcmp(sub_schurs->mat_solver_type, MATSOLVERMUMPS, &ismumps));
2120: PetscCall(PetscStrcmp(sub_schurs->mat_solver_type, MATSOLVERMKL_PARDISO, &ispardiso));
2121: sub_schurs->schur_explicit = (PetscBool)(ispardiso || ismumps);
2123: /* for reals, symmetric and Hermitian are synonyms */
2124: #if !PetscDefined(USE_COMPLEX)
2125: sub_schurs->is_symmetric = (PetscBool)(sub_schurs->is_symmetric && sub_schurs->is_hermitian);
2126: sub_schurs->is_hermitian = sub_schurs->is_symmetric;
2127: #endif
2129: PetscCall(PetscObjectReference((PetscObject)is_I));
2130: sub_schurs->is_I = is_I;
2131: PetscCall(PetscObjectReference((PetscObject)is_B));
2132: sub_schurs->is_B = is_B;
2133: PetscCall(PetscObjectReference((PetscObject)graph->l2gmap));
2134: sub_schurs->l2gmap = graph->l2gmap;
2135: PetscCall(PetscObjectReference((PetscObject)BtoNmap));
2136: sub_schurs->BtoNmap = BtoNmap;
2137: sub_schurs->n_subs = n_all_cc;
2138: sub_schurs->is_subs = all_cc;
2139: sub_schurs->S_Ej_all = NULL;
2140: sub_schurs->sum_S_Ej_all = NULL;
2141: sub_schurs->sum_S_Ej_inv_all = NULL;
2142: sub_schurs->sum_S_Ej_tilda_all = NULL;
2143: sub_schurs->is_Ej_all = NULL;
2144: PetscFunctionReturn(PETSC_SUCCESS);
2145: }
2147: PetscErrorCode PCBDDCSubSchursCreate(PCBDDCSubSchurs *sub_schurs)
2148: {
2149: PCBDDCSubSchurs schurs_ctx;
2151: PetscFunctionBegin;
2152: PetscCall(PetscNew(&schurs_ctx));
2153: schurs_ctx->n_subs = 0;
2154: *sub_schurs = schurs_ctx;
2155: PetscFunctionReturn(PETSC_SUCCESS);
2156: }
2158: PetscErrorCode PCBDDCSubSchursReset(PCBDDCSubSchurs sub_schurs)
2159: {
2160: PetscFunctionBegin;
2161: if (!sub_schurs) PetscFunctionReturn(PETSC_SUCCESS);
2162: sub_schurs->graph = NULL;
2163: PetscCall(PetscFree(sub_schurs->prefix));
2164: PetscCall(MatDestroy(&sub_schurs->A));
2165: PetscCall(MatDestroy(&sub_schurs->S));
2166: PetscCall(ISDestroy(&sub_schurs->is_I));
2167: PetscCall(ISDestroy(&sub_schurs->is_B));
2168: PetscCall(ISLocalToGlobalMappingDestroy(&sub_schurs->l2gmap));
2169: PetscCall(ISLocalToGlobalMappingDestroy(&sub_schurs->BtoNmap));
2170: PetscCall(MatDestroy(&sub_schurs->S_Ej_all));
2171: PetscCall(MatDestroy(&sub_schurs->sum_S_Ej_all));
2172: PetscCall(MatDestroy(&sub_schurs->sum_S_Ej_inv_all));
2173: PetscCall(MatDestroy(&sub_schurs->sum_S_Ej_tilda_all));
2174: PetscCall(ISDestroy(&sub_schurs->is_Ej_all));
2175: PetscCall(ISDestroy(&sub_schurs->is_vertices));
2176: PetscCall(ISDestroy(&sub_schurs->is_dir));
2177: PetscCall(PetscBTDestroy(&sub_schurs->is_edge));
2178: for (PetscInt i = 0; i < sub_schurs->n_subs; i++) PetscCall(ISDestroy(&sub_schurs->is_subs[i]));
2179: if (sub_schurs->n_subs) PetscCall(PetscFree(sub_schurs->is_subs));
2180: if (sub_schurs->reuse_solver) PetscCall(PCBDDCReuseSolversReset(sub_schurs->reuse_solver));
2181: PetscCall(PetscFree(sub_schurs->reuse_solver));
2182: if (sub_schurs->change) {
2183: for (PetscInt i = 0; i < sub_schurs->n_subs; i++) {
2184: PetscCall(KSPDestroy(&sub_schurs->change[i]));
2185: PetscCall(ISDestroy(&sub_schurs->change_primal_sub[i]));
2186: }
2187: }
2188: PetscCall(PetscFree(sub_schurs->change));
2189: PetscCall(PetscFree(sub_schurs->change_primal_sub));
2190: sub_schurs->n_subs = 0;
2191: PetscFunctionReturn(PETSC_SUCCESS);
2192: }
2194: PetscErrorCode PCBDDCSubSchursDestroy(PCBDDCSubSchurs *sub_schurs)
2195: {
2196: PetscFunctionBegin;
2197: PetscCall(PCBDDCSubSchursReset(*sub_schurs));
2198: PetscCall(PetscFree(*sub_schurs));
2199: PetscFunctionReturn(PETSC_SUCCESS);
2200: }
2202: static inline PetscErrorCode PCBDDCAdjGetNextLayer_Private(PetscInt *queue_tip, PetscInt n_prev, PetscBT touched, PetscInt *xadj, PetscInt *adjncy, PetscInt *n_added)
2203: {
2204: PetscInt i, j, n;
2206: PetscFunctionBegin;
2207: n = 0;
2208: for (i = -n_prev; i < 0; i++) {
2209: PetscInt start_dof = queue_tip[i];
2210: for (j = xadj[start_dof]; j < xadj[start_dof + 1]; j++) {
2211: PetscInt dof = adjncy[j];
2212: if (!PetscBTLookup(touched, dof)) {
2213: PetscCall(PetscBTSet(touched, dof));
2214: queue_tip[n] = dof;
2215: n++;
2216: }
2217: }
2218: }
2219: *n_added = n;
2220: PetscFunctionReturn(PETSC_SUCCESS);
2221: }