Actual source code: factorschur.c
1: #include <petsc/private/matimpl.h>
2: #include <../src/mat/impls/dense/seq/dense.h>
4: PETSC_INTERN PetscErrorCode MatFactorSetUpInPlaceSchur_Private(Mat F)
5: {
6: Mat St, S = F->schur;
7: MatFactorInfo info;
9: MatSetUnfactored(S);
10: MatGetFactor(S, S->solvertype ? S->solvertype : MATSOLVERPETSC, F->factortype, &St);
11: if (St->factortype == MAT_FACTOR_CHOLESKY) { /* LDL^t regarded as Cholesky */
12: MatCholeskyFactorSymbolic(St, S, NULL, &info);
13: } else {
14: MatLUFactorSymbolic(St, S, NULL, NULL, &info);
15: }
16: S->ops->solve = St->ops->solve;
17: S->ops->matsolve = St->ops->matsolve;
18: S->ops->solvetranspose = St->ops->solvetranspose;
19: S->ops->matsolvetranspose = St->ops->matsolvetranspose;
20: S->ops->solveadd = St->ops->solveadd;
21: S->ops->solvetransposeadd = St->ops->solvetransposeadd;
23: MatDestroy(&St);
24: return 0;
25: }
27: PETSC_INTERN PetscErrorCode MatFactorUpdateSchurStatus_Private(Mat F)
28: {
29: Mat S = F->schur;
31: switch (F->schur_status) {
32: case MAT_FACTOR_SCHUR_UNFACTORED:
33: case MAT_FACTOR_SCHUR_INVERTED:
34: if (S) {
35: S->ops->solve = NULL;
36: S->ops->matsolve = NULL;
37: S->ops->solvetranspose = NULL;
38: S->ops->matsolvetranspose = NULL;
39: S->ops->solveadd = NULL;
40: S->ops->solvetransposeadd = NULL;
41: S->factortype = MAT_FACTOR_NONE;
42: PetscFree(S->solvertype);
43: }
44: break;
45: case MAT_FACTOR_SCHUR_FACTORED:
46: break;
47: default:
48: SETERRQ(PetscObjectComm((PetscObject)F), PETSC_ERR_SUP, "Unhandled MatFactorSchurStatus %d", F->schur_status);
49: }
50: return 0;
51: }
53: /* Schur status updated in the interface */
54: PETSC_INTERN PetscErrorCode MatFactorFactorizeSchurComplement_Private(Mat F)
55: {
56: MatFactorInfo info;
58: PetscLogEventBegin(MAT_FactorFactS, F, 0, 0, 0);
59: if (F->factortype == MAT_FACTOR_CHOLESKY) { /* LDL^t regarded as Cholesky */
60: MatCholeskyFactor(F->schur, NULL, &info);
61: } else {
62: MatLUFactor(F->schur, NULL, NULL, &info);
63: }
64: PetscLogEventEnd(MAT_FactorFactS, F, 0, 0, 0);
65: return 0;
66: }
68: /* Schur status updated in the interface */
69: PETSC_INTERN PetscErrorCode MatFactorInvertSchurComplement_Private(Mat F)
70: {
71: Mat S = F->schur;
73: if (S) {
74: PetscMPIInt size;
75: PetscBool isdense, isdensecuda;
77: MPI_Comm_size(PetscObjectComm((PetscObject)S), &size);
79: PetscObjectTypeCompare((PetscObject)S, MATSEQDENSE, &isdense);
80: PetscObjectTypeCompare((PetscObject)S, MATSEQDENSECUDA, &isdensecuda);
81: PetscLogEventBegin(MAT_FactorInvS, F, 0, 0, 0);
82: if (isdense) {
83: MatSeqDenseInvertFactors_Private(S);
84: #if defined(PETSC_HAVE_CUDA)
85: } else if (isdensecuda) {
86: MatSeqDenseCUDAInvertFactors_Private(S);
87: #endif
88: } else SETERRQ(PetscObjectComm((PetscObject)S), PETSC_ERR_SUP, "Not implemented for type %s", ((PetscObject)S)->type_name);
89: PetscLogEventEnd(MAT_FactorInvS, F, 0, 0, 0);
90: }
91: return 0;
92: }