Actual source code: factorschur.c
petsc-3.13.6 2020-09-29
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;
25: MatDestroy(&St);
26: return(0);
27: }
29: PETSC_INTERN PetscErrorCode MatFactorUpdateSchurStatus_Private(Mat F)
30: {
31: Mat S = F->schur;
35: switch(F->schur_status) {
36: case MAT_FACTOR_SCHUR_UNFACTORED:
37: case MAT_FACTOR_SCHUR_INVERTED:
38: if (S) {
39: S->ops->solve = NULL;
40: S->ops->matsolve = NULL;
41: S->ops->solvetranspose = NULL;
42: S->ops->matsolvetranspose = NULL;
43: S->ops->solveadd = NULL;
44: S->ops->solvetransposeadd = NULL;
45: S->factortype = MAT_FACTOR_NONE;
46: PetscFree(S->solvertype);
47: }
48: break;
49: case MAT_FACTOR_SCHUR_FACTORED:
50: break;
51: default:
52: SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_SUP,"Unhandled MatFactorSchurStatus %D",F->schur_status);
53: }
54: return(0);
55: }
57: /* Schur status updated in the interface */
58: PETSC_INTERN PetscErrorCode MatFactorFactorizeSchurComplement_Private(Mat F)
59: {
60: MatFactorInfo info;
64: PetscLogEventBegin(MAT_FactorFactS,F,0,0,0);
65: if (F->factortype == MAT_FACTOR_CHOLESKY) { /* LDL^t regarded as Cholesky */
66: MatCholeskyFactor(F->schur,NULL,&info);
67: } else {
68: MatLUFactor(F->schur,NULL,NULL,&info);
69: }
70: PetscLogEventEnd(MAT_FactorFactS,F,0,0,0);
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,isdensecuda;
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: PetscObjectTypeCompare((PetscObject)S,MATSEQDENSECUDA,&isdensecuda);
89: PetscLogEventBegin(MAT_FactorInvS,F,0,0,0);
90: if (isdense) {
91: MatSeqDenseInvertFactors_Private(S);
92: #if defined(PETSC_HAVE_CUDA)
93: } else if (isdensecuda) {
94: MatSeqDenseCUDAInvertFactors_Private(S);
95: #endif
96: } else SETERRQ1(PetscObjectComm((PetscObject)S),PETSC_ERR_SUP,"Not implemented for type %s",((PetscObject)S)->type_name);
97: PetscLogEventEnd(MAT_FactorInvS,F,0,0,0);
98: }
99: return(0);
100: }