Actual source code: factorschur.c
petsc-3.10.5 2019-03-28
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;
8: PetscErrorCode ierr;
11: MatSetUnfactored(S);
12: MatGetFactor(S,S->solvertype ? S->solvertype : MATSOLVERPETSC,F->factortype,&St);
13: if (St->factortype == MAT_FACTOR_CHOLESKY) { /* LDL^t regarded as Cholesky */
14: MatCholeskyFactorSymbolic(St,S,NULL,&info);
15: } else {
16: MatLUFactorSymbolic(St,S,NULL,NULL,&info);
17: }
18: S->ops->solve = St->ops->solve;
19: S->ops->matsolve = St->ops->matsolve;
20: S->ops->solvetranspose = St->ops->solvetranspose;
21: S->ops->matsolvetranspose = St->ops->matsolvetranspose;
22: S->ops->solveadd = St->ops->solveadd;
23: S->ops->solvetransposeadd = St->ops->solvetransposeadd;
24: S->factortype = St->factortype;
26: MatDestroy(&St);
27: return(0);
28: }
30: PETSC_INTERN PetscErrorCode MatFactorUpdateSchurStatus_Private(Mat F)
31: {
32: Mat S = F->schur;
36: switch(F->schur_status) {
37: case MAT_FACTOR_SCHUR_UNFACTORED:
38: case MAT_FACTOR_SCHUR_INVERTED:
39: if (S) {
40: S->ops->solve = NULL;
41: S->ops->matsolve = NULL;
42: S->ops->solvetranspose = NULL;
43: S->ops->matsolvetranspose = NULL;
44: S->ops->solveadd = NULL;
45: S->ops->solvetransposeadd = NULL;
46: S->factortype = MAT_FACTOR_NONE;
47: PetscFree(S->solvertype);
48: }
49: break;
50: case MAT_FACTOR_SCHUR_FACTORED:
51: MatFactorSetUpInPlaceSchur_Private(F);
52: break;
53: default:
54: SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_SUP,"Unhandled MatFactorSchurStatus %D",F->schur_status);
55: }
56: return(0);
57: }
59: /* Schur status updated in the interface */
60: PETSC_INTERN PetscErrorCode MatFactorFactorizeSchurComplement_Private(Mat F)
61: {
62: MatFactorInfo info;
66: if (F->factortype == MAT_FACTOR_CHOLESKY) { /* LDL^t regarded as Cholesky */
67: MatCholeskyFactor(F->schur,NULL,&info);
68: } else {
69: MatLUFactor(F->schur,NULL,NULL,&info);
70: }
71: return(0);
72: }
74: /* Schur status updated in the interface */
75: PETSC_INTERN PetscErrorCode MatFactorInvertSchurComplement_Private(Mat F)
76: {
77: Mat S = F->schur;
80: if (S) {
81: PetscMPIInt size;
82: PetscBool isdense;
85: MPI_Comm_size(PetscObjectComm((PetscObject)S),&size);
86: if (size > 1) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_SUP,"Not yet implemented");
87: PetscObjectTypeCompare((PetscObject)S,MATSEQDENSE,&isdense);
88: if (!isdense) SETERRQ1(PetscObjectComm((PetscObject)S),PETSC_ERR_SUP,"Not implemented for type %s",((PetscObject)S)->type_name);
89: MatSeqDenseInvertFactors_Private(S);
90: }
91: return(0);
92: }