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: }