Actual source code: htoolm.c
1: #include <petsc/private/kspimpl.h>
2: #include <petsc/private/matimpl.h>
4: static PetscErrorCode MatFactorReset_NestHtool(Mat F)
5: {
6: PetscFunctionBegin;
7: PetscCall(PetscObjectCompose((PetscObject)F, "MatFactorKSP", NULL));
8: PetscFunctionReturn(PETSC_SUCCESS);
9: }
11: static PetscErrorCode MatFactorDestroy_NestHtool(Mat F)
12: {
13: PetscFunctionBegin;
14: PetscCall(MatFactorReset_NestHtool(F));
15: PetscCall(PetscObjectComposeFunction((PetscObject)F, "MatFactorGetSolverType_C", NULL));
16: PetscFunctionReturn(PETSC_SUCCESS);
17: }
19: static PetscErrorCode MatFactorGetSolverType_NestHtool(PETSC_UNUSED Mat F, MatSolverType *type)
20: {
21: PetscFunctionBegin;
22: *type = MATSOLVERHTOOL;
23: PetscFunctionReturn(PETSC_SUCCESS);
24: }
26: static PetscErrorCode MatFactorSymbolic_NestHtool_FieldSplit(Mat F, Mat A)
27: {
28: KSP ksp;
29: PC pc;
30: Mat **mats;
31: IS *rowis, *colis, *ises = NULL;
32: IS is_h, is_r = NULL;
33: PetscBool found = PETSC_FALSE, same;
34: PetscInt nr, nc, hidx = -1, nrest = 0;
36: PetscFunctionBegin;
37: PetscCall(MatNestGetSubMats(A, &nr, &nc, &mats));
38: PetscCheck(nr == nc, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only square MATNEST matrices are supported, got (%" PetscInt_FMT ",%" PetscInt_FMT ")", nr, nc);
39: PetscCall(PetscMalloc2(nr, &rowis, nc, &colis));
40: PetscCall(MatNestGetISs(A, rowis, colis));
41: for (PetscInt i = 0; i < nr; ++i) {
42: PetscCall(ISEqualUnsorted(rowis[i], colis[i], &same));
43: PetscCheck(same, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only congruent MATNEST row/column layouts are supported");
44: }
45: for (PetscInt i = 0; i < nr; ++i) {
46: PetscBool flg;
48: PetscCall(PetscObjectTypeCompare((PetscObject)mats[i][i], MATHTOOL, &flg));
49: if (flg) {
50: PetscCheck(!found, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Exactly one diagonal MATHTOOL block is required");
51: found = PETSC_TRUE;
52: hidx = i;
53: } else {
54: PetscCall(PetscObjectTypeCompareAny((PetscObject)mats[i][i], &flg, MATSEQDENSE, MATMPIDENSE, NULL));
55: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Diagonal block %" PetscInt_FMT " must be MATHTOOL or MATDENSE", i);
56: }
57: }
58: PetscCheck(found, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No diagonal MATHTOOL block found");
59: PetscCall(ISDuplicate(rowis[hidx], &is_h));
60: nrest = nr - 1;
61: if (nrest > 0) {
62: PetscCall(PetscMalloc1(nrest, &ises));
63: for (PetscInt i = 0, j = 0; i < nr; ++i) {
64: if (i == hidx) continue;
65: ises[j++] = rowis[i];
66: }
67: PetscCall(ISConcatenate(PetscObjectComm((PetscObject)A), nrest, ises, &is_r));
68: PetscCall(PetscFree(ises));
69: }
70: PetscCall(PetscFree2(rowis, colis));
71: PetscCall(PetscObjectCompose((PetscObject)F, "MatFactorKSP", NULL));
72: PetscCall(KSPCreate(PetscObjectComm((PetscObject)A), &ksp));
73: PetscCall(KSPSetType(ksp, KSPPREONLY));
74: PetscCall(KSPGetPC(ksp, &pc));
75: PetscCall(PCSetType(pc, PCFIELDSPLIT));
76: PetscCall(PCFieldSplitSetIS(pc, NULL, is_h));
77: PetscCall(ISDestroy(&is_h));
78: if (nrest > 0) {
79: PetscCall(PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR));
80: PetscCall(PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_FULL, NULL));
81: PetscCall(PCFieldSplitSetIS(pc, NULL, is_r));
82: PetscCall(ISDestroy(&is_r));
83: }
84: PetscCall(PetscObjectCompose((PetscObject)F, "MatFactorKSP", (PetscObject)ksp));
85: PetscCall(KSPDestroy(&ksp));
86: PetscFunctionReturn(PETSC_SUCCESS);
87: }
89: static PetscErrorCode MatFactorNumeric_NestHtool(Mat F, Mat A, MatFactorType factor_type)
90: {
91: KSP ksp;
92: KSP *subksp;
93: PC pc;
94: PetscInt nsplits;
96: PetscFunctionBegin;
97: PetscCall(PetscObjectQuery((PetscObject)F, "MatFactorKSP", (PetscObject *)&ksp));
98: PetscCheck(ksp, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Outer PCFIELDSPLIT has not been created yet, call symbolic factorization first");
99: PetscCall(KSPSetOperators(ksp, A, A));
100: PetscCall(KSPSetUp(ksp));
101: PetscCall(KSPGetPC(ksp, &pc));
102: PetscCall(PCFieldSplitGetSubKSP(pc, &nsplits, &subksp));
103: PetscCheck(nsplits == 1 || nsplits == 2, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Unexpected number of splits %" PetscInt_FMT " (!= 1 or 2)", nsplits);
104: for (PetscInt i = 0; i < nsplits; ++i) {
105: Mat subA;
106: PetscBool flg;
108: PetscCall(KSPSetType(subksp[i], KSPPREONLY));
109: PetscCall(KSPGetPC(subksp[i], &pc));
110: PetscCall(KSPGetOperators(subksp[i], &subA, NULL));
111: PetscCall(PetscObjectTypeCompare((PetscObject)subA, MATHTOOL, &flg));
112: if (flg) flg = (PetscBool)(subA->symmetric == PETSC_BOOL3_TRUE || subA->hermitian == PETSC_BOOL3_TRUE);
113: PetscCall(PCSetType(pc, factor_type == MAT_FACTOR_CHOLESKY || flg ? PCCHOLESKY : PCLU));
114: }
115: PetscCall(PetscFree(subksp));
116: PetscFunctionReturn(PETSC_SUCCESS);
117: }
119: static PetscErrorCode MatSolve_NestHtool(Mat F, Vec b, Vec x)
120: {
121: KSP ksp;
123: PetscFunctionBegin;
124: PetscCall(PetscObjectQuery((PetscObject)F, "MatFactorKSP", (PetscObject *)&ksp));
125: PetscCall(KSPSolve(ksp, b, x));
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: static PetscErrorCode MatSolveTranspose_NestHtool(Mat F, Vec b, Vec x)
130: {
131: KSP ksp;
133: PetscFunctionBegin;
134: PetscCall(PetscObjectQuery((PetscObject)F, "MatFactorKSP", (PetscObject *)&ksp));
135: PetscCall(KSPSolveTranspose(ksp, b, x));
136: PetscFunctionReturn(PETSC_SUCCESS);
137: }
139: static PetscErrorCode MatMatSolve_NestHtool(Mat F, Mat B, Mat X)
140: {
141: KSP ksp;
143: PetscFunctionBegin;
144: PetscCall(PetscObjectQuery((PetscObject)F, "MatFactorKSP", (PetscObject *)&ksp));
145: PetscCall(KSPMatSolve(ksp, B, X));
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: static PetscErrorCode MatMatSolveTranspose_NestHtool(Mat F, Mat B, Mat X)
150: {
151: KSP ksp;
153: PetscFunctionBegin;
154: PetscCall(PetscObjectQuery((PetscObject)F, "MatFactorKSP", (PetscObject *)&ksp));
155: PetscCall(KSPMatSolveTranspose(ksp, B, X));
156: PetscFunctionReturn(PETSC_SUCCESS);
157: }
159: static PetscErrorCode MatLUFactorNumeric_NestHtool(Mat F, Mat A, PETSC_UNUSED const MatFactorInfo *info)
160: {
161: PetscFunctionBegin;
162: PetscCall(MatFactorNumeric_NestHtool(F, A, MAT_FACTOR_LU));
163: F->ops->solve = MatSolve_NestHtool;
164: F->ops->matsolve = MatMatSolve_NestHtool;
165: F->ops->solvetranspose = MatSolveTranspose_NestHtool;
166: F->ops->matsolvetranspose = MatMatSolveTranspose_NestHtool;
167: PetscFunctionReturn(PETSC_SUCCESS);
168: }
170: static PetscErrorCode MatLUFactorSymbolic_NestHtool(Mat F, Mat A, PETSC_UNUSED IS row, PETSC_UNUSED IS col, PETSC_UNUSED const MatFactorInfo *info)
171: {
172: PetscFunctionBegin;
173: PetscCall(MatFactorSymbolic_NestHtool_FieldSplit(F, A));
174: F->ops->lufactornumeric = MatLUFactorNumeric_NestHtool;
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: static PetscErrorCode MatCholeskyFactorNumeric_NestHtool(Mat F, Mat A, PETSC_UNUSED const MatFactorInfo *info)
179: {
180: PetscFunctionBegin;
181: PetscCall(MatFactorNumeric_NestHtool(F, A, MAT_FACTOR_CHOLESKY));
182: F->ops->solve = MatSolve_NestHtool;
183: F->ops->matsolve = MatMatSolve_NestHtool;
184: F->ops->solvetranspose = MatSolveTranspose_NestHtool;
185: F->ops->matsolvetranspose = MatMatSolveTranspose_NestHtool;
186: PetscFunctionReturn(PETSC_SUCCESS);
187: }
189: static PetscErrorCode MatCholeskyFactorSymbolic_NestHtool(Mat F, Mat A, PETSC_UNUSED IS row, PETSC_UNUSED const MatFactorInfo *info)
190: {
191: PetscFunctionBegin;
192: PetscCall(MatFactorSymbolic_NestHtool_FieldSplit(F, A));
193: F->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_NestHtool;
194: PetscFunctionReturn(PETSC_SUCCESS);
195: }
197: static PetscErrorCode MatGetFactor_nest_htool(Mat A, MatFactorType ftype, Mat *F)
198: {
199: Mat B;
200: PetscMPIInt size;
202: PetscFunctionBegin;
203: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
204: PetscCheck(size == 1, PetscObjectComm((PetscObject)A), PETSC_ERR_WRONG_MPI_SIZE, "Unsupported parallel MatGetFactor()");
205: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
206: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
207: PetscCall(PetscStrallocpy(MATSOLVERHTOOL, &((PetscObject)B)->type_name));
208: PetscCall(MatSetUp(B));
210: B->ops->destroy = MatFactorDestroy_NestHtool;
211: B->ops->getinfo = MatGetInfo_External;
212: B->factortype = ftype;
213: B->trivialsymbolic = PETSC_FALSE;
214: B->preallocated = PETSC_TRUE;
215: B->assembled = PETSC_TRUE;
217: PetscCheck(ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_CHOLESKY, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only MAT_FACTOR_LU and MAT_FACTOR_CHOLESKY are supported");
218: if (ftype == MAT_FACTOR_LU) B->ops->lufactorsymbolic = MatLUFactorSymbolic_NestHtool;
219: else B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_NestHtool;
221: PetscCall(PetscFree(B->solvertype));
222: PetscCall(PetscStrallocpy(MATSOLVERHTOOL, &B->solvertype));
224: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_NestHtool));
225: *F = B;
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }
229: PETSC_INTERN PetscErrorCode MatSolverTypeRegister_NestHtool(void)
230: {
231: PetscFunctionBegin;
232: PetscCall(MatSolverTypeRegister(MATSOLVERHTOOL, MATNEST, MAT_FACTOR_LU, MatGetFactor_nest_htool));
233: PetscCall(MatSolverTypeRegister(MATSOLVERHTOOL, MATNEST, MAT_FACTOR_CHOLESKY, MatGetFactor_nest_htool));
234: PetscFunctionReturn(PETSC_SUCCESS);
235: }