/*
   Defines a  (S)SOR  preconditioner for any Mat implementation
*/
#include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/

typedef struct {
  PetscInt   its;  /* inner iterations, number of sweeps */
  PetscInt   lits; /* local inner iterations, number of sweeps applied by the local matrix mat->A */
  MatSORType sym;  /* forward, reverse, symmetric etc. */
  PetscReal  omega;
  PetscReal  fshift;
} PC_SOR;

static PetscErrorCode PCDestroy_SOR(PC pc)
{
  PetscFunctionBegin;
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORSetSymmetric_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORSetOmega_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORSetIterations_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORGetSymmetric_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORGetOmega_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORGetIterations_C", NULL));
  PetscCall(PetscFree(pc->data));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCApply_SOR(PC pc, Vec x, Vec y)
{
  PC_SOR  *jac  = (PC_SOR *)pc->data;
  PetscInt flag = jac->sym | SOR_ZERO_INITIAL_GUESS;

  PetscFunctionBegin;
  PetscCall(MatSOR(pc->pmat, x, jac->omega, (MatSORType)flag, jac->fshift, jac->its, jac->lits, y));
  PetscCall(MatFactorGetError(pc->pmat, (MatFactorError *)&pc->failedreason));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCApplyTranspose_SOR(PC pc, Vec x, Vec y)
{
  PC_SOR   *jac  = (PC_SOR *)pc->data;
  PetscInt  flag = jac->sym | SOR_ZERO_INITIAL_GUESS;
  PetscBool set, sym;

  PetscFunctionBegin;
  PetscCall(MatIsSymmetricKnown(pc->pmat, &set, &sym));
  PetscCheck(set && sym && (jac->sym == SOR_SYMMETRIC_SWEEP || jac->sym == SOR_LOCAL_SYMMETRIC_SWEEP), PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Can only apply transpose of SOR if matrix is symmetric and sweep is symmetric");
  PetscCall(MatSOR(pc->pmat, x, jac->omega, (MatSORType)flag, jac->fshift, jac->its, jac->lits, y));
  PetscCall(MatFactorGetError(pc->pmat, (MatFactorError *)&pc->failedreason));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCApplyRichardson_SOR(PC pc, Vec b, Vec y, Vec w, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt its, PetscBool guesszero, PetscInt *outits, PCRichardsonConvergedReason *reason)
{
  PC_SOR    *jac   = (PC_SOR *)pc->data;
  MatSORType stype = jac->sym;

  PetscFunctionBegin;
  if (guesszero) stype = (MatSORType)(stype | SOR_ZERO_INITIAL_GUESS);
  PetscCall(MatSOR(pc->pmat, b, jac->omega, stype, jac->fshift, its * jac->its, jac->lits, y));
  PetscCall(MatFactorGetError(pc->pmat, (MatFactorError *)&pc->failedreason));
  *outits = its;
  *reason = PCRICHARDSON_CONVERGED_ITS;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCSetFromOptions_SOR(PC pc, PetscOptionItems PetscOptionsObject)
{
  PC_SOR   *jac = (PC_SOR *)pc->data;
  PetscBool flg;
  PetscReal omega;

  PetscFunctionBegin;
  PetscOptionsHeadBegin(PetscOptionsObject, "(S)SOR options");
  PetscCall(PetscOptionsReal("-pc_sor_omega", "relaxation factor (0 < omega < 2)", "PCSORSetOmega", jac->omega, &omega, &flg));
  if (flg) PetscCall(PCSORSetOmega(pc, omega));
  PetscCall(PetscOptionsReal("-pc_sor_diagonal_shift", "Add to the diagonal entries", "", jac->fshift, &jac->fshift, NULL));
  PetscCall(PetscOptionsInt("-pc_sor_its", "number of inner SOR iterations", "PCSORSetIterations", jac->its, &jac->its, NULL));
  PetscCall(PetscOptionsInt("-pc_sor_lits", "number of local inner SOR iterations", "PCSORSetIterations", jac->lits, &jac->lits, NULL));
  PetscCall(PetscOptionsBoolGroupBegin("-pc_sor_symmetric", "SSOR, not SOR", "PCSORSetSymmetric", &flg));
  if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_SYMMETRIC_SWEEP));
  PetscCall(PetscOptionsBoolGroup("-pc_sor_backward", "use backward sweep instead of forward", "PCSORSetSymmetric", &flg));
  if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_BACKWARD_SWEEP));
  PetscCall(PetscOptionsBoolGroup("-pc_sor_forward", "use forward sweep", "PCSORSetSymmetric", &flg));
  if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_FORWARD_SWEEP));
  PetscCall(PetscOptionsBoolGroup("-pc_sor_local_symmetric", "use SSOR separately on each processor", "PCSORSetSymmetric", &flg));
  if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_LOCAL_SYMMETRIC_SWEEP));
  PetscCall(PetscOptionsBoolGroup("-pc_sor_local_backward", "use backward sweep locally", "PCSORSetSymmetric", &flg));
  if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_LOCAL_BACKWARD_SWEEP));
  PetscCall(PetscOptionsBoolGroupEnd("-pc_sor_local_forward", "use forward sweep locally", "PCSORSetSymmetric", &flg));
  if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_LOCAL_FORWARD_SWEEP));
  PetscOptionsHeadEnd();
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCView_SOR(PC pc, PetscViewer viewer)
{
  PC_SOR     *jac = (PC_SOR *)pc->data;
  MatSORType  sym = jac->sym;
  const char *sortype;
  PetscBool   isascii;

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
  if (isascii) {
    if (sym & SOR_ZERO_INITIAL_GUESS) PetscCall(PetscViewerASCIIPrintf(viewer, "  zero initial guess\n"));
    if (sym == SOR_APPLY_UPPER) sortype = "apply_upper";
    else if (sym == SOR_APPLY_LOWER) sortype = "apply_lower";
    else if (sym & SOR_EISENSTAT) sortype = "Eisenstat";
    else if ((sym & SOR_SYMMETRIC_SWEEP) == SOR_SYMMETRIC_SWEEP) sortype = "symmetric";
    else if (sym & SOR_BACKWARD_SWEEP) sortype = "backward";
    else if (sym & SOR_FORWARD_SWEEP) sortype = "forward";
    else if ((sym & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) sortype = "local_symmetric";
    else if (sym & SOR_LOCAL_FORWARD_SWEEP) sortype = "local_forward";
    else if (sym & SOR_LOCAL_BACKWARD_SWEEP) sortype = "local_backward";
    else sortype = "unknown";
    PetscCall(PetscViewerASCIIPrintf(viewer, "  type = %s, iterations = %" PetscInt_FMT ", local iterations = %" PetscInt_FMT ", omega = %g\n", sortype, jac->its, jac->lits, (double)jac->omega));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCSORSetSymmetric_SOR(PC pc, MatSORType flag)
{
  PC_SOR *jac = (PC_SOR *)pc->data;

  PetscFunctionBegin;
  jac->sym = flag;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCSORSetOmega_SOR(PC pc, PetscReal omega)
{
  PC_SOR *jac = (PC_SOR *)pc->data;

  PetscFunctionBegin;
  PetscCheck(omega > 0.0 && omega < 2.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Relaxation out of range");
  jac->omega = omega;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCSORSetIterations_SOR(PC pc, PetscInt its, PetscInt lits)
{
  PC_SOR *jac = (PC_SOR *)pc->data;

  PetscFunctionBegin;
  jac->its  = its;
  jac->lits = lits;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCSORGetSymmetric_SOR(PC pc, MatSORType *flag)
{
  PC_SOR *jac = (PC_SOR *)pc->data;

  PetscFunctionBegin;
  *flag = jac->sym;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCSORGetOmega_SOR(PC pc, PetscReal *omega)
{
  PC_SOR *jac = (PC_SOR *)pc->data;

  PetscFunctionBegin;
  *omega = jac->omega;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCSORGetIterations_SOR(PC pc, PetscInt *its, PetscInt *lits)
{
  PC_SOR *jac = (PC_SOR *)pc->data;

  PetscFunctionBegin;
  if (its) *its = jac->its;
  if (lits) *lits = jac->lits;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCSORGetSymmetric - Gets the form the SOR preconditioner is using;   backward, or forward relaxation.  The local variants perform SOR on
  each processor.  By default forward relaxation is used.

  Logically Collective

  Input Parameter:
. pc - the preconditioner context

  Output Parameter:
. flag - one of the following
.vb
    SOR_FORWARD_SWEEP
    SOR_BACKWARD_SWEEP
    SOR_SYMMETRIC_SWEEP
    SOR_LOCAL_FORWARD_SWEEP
    SOR_LOCAL_BACKWARD_SWEEP
    SOR_LOCAL_SYMMETRIC_SWEEP
.ve

  Options Database Keys:
+ -pc_sor_symmetric       - Activates symmetric version
. -pc_sor_backward        - Activates backward version
. -pc_sor_local_forward   - Activates local forward version
. -pc_sor_local_symmetric - Activates local symmetric version
- -pc_sor_local_backward  - Activates local backward version

  Note:
  To use the Eisenstat trick with SSOR, employ the `PCEISENSTAT` preconditioner,
  which can be chosen with the option
.  -pc_type eisenstat - Activates Eisenstat trick

  Level: intermediate

.seealso: [](ch_ksp), `PCSOR`, `PCEisenstatSetOmega()`, `PCSORSetIterations()`, `PCSORSetOmega()`, `PCSORSetSymmetric()`
@*/
PetscErrorCode PCSORGetSymmetric(PC pc, MatSORType *flag)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscUseMethod(pc, "PCSORGetSymmetric_C", (PC, MatSORType *), (pc, flag));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCSORGetOmega - Gets the SOR relaxation coefficient, omega
  (where omega = 1.0 by default).

  Logically Collective

  Input Parameter:
. pc - the preconditioner context

  Output Parameter:
. omega - relaxation coefficient (0 < omega < 2).

  Options Database Key:
. -pc_sor_omega <omega> - Sets omega

  Level: intermediate

.seealso: [](ch_ksp), `PCSOR`, `PCSORSetSymmetric()`, `PCSORSetIterations()`, `PCEisenstatSetOmega()`, `PCSORSetOmega()`
@*/
PetscErrorCode PCSORGetOmega(PC pc, PetscReal *omega)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscUseMethod(pc, "PCSORGetOmega_C", (PC, PetscReal *), (pc, omega));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCSORGetIterations - Gets the number of inner iterations to
  be used by the SOR preconditioner. The default is 1.

  Logically Collective

  Input Parameter:
. pc - the preconditioner context

  Output Parameters:
+ lits - number of local iterations, smoothings over just variables on processor
- its  - number of parallel iterations to use; each parallel iteration has lits local iterations

  Options Database Keys:
+ -pc_sor_its <its>   - Sets number of iterations
- -pc_sor_lits <lits> - Sets number of local iterations

  Level: intermediate

  Note:
  When run on one processor the number of smoothings is lits*its

.seealso: [](ch_ksp), `PCSOR`, `PCSORSetOmega()`, `PCSORSetSymmetric()`, `PCSORSetIterations()`
@*/
PetscErrorCode PCSORGetIterations(PC pc, PetscInt *its, PetscInt *lits)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscUseMethod(pc, "PCSORGetIterations_C", (PC, PetscInt *, PetscInt *), (pc, its, lits));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCSORSetSymmetric - Sets the SOR preconditioner to use symmetric (SSOR),
  backward, or forward relaxation.  The local variants perform SOR on
  each processor.  By default forward relaxation is used.

  Logically Collective

  Input Parameters:
+ pc   - the preconditioner context
- flag - one of the following
.vb
    SOR_FORWARD_SWEEP
    SOR_BACKWARD_SWEEP
    SOR_SYMMETRIC_SWEEP
    SOR_LOCAL_FORWARD_SWEEP
    SOR_LOCAL_BACKWARD_SWEEP
    SOR_LOCAL_SYMMETRIC_SWEEP
.ve

  Options Database Keys:
+ -pc_sor_symmetric       - Activates symmetric version
. -pc_sor_backward        - Activates backward version
. -pc_sor_local_forward   - Activates local forward version
. -pc_sor_local_symmetric - Activates local symmetric version
- -pc_sor_local_backward  - Activates local backward version

  Note:
  To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
  which can be chosen with the option
.  -pc_type eisenstat - Activates Eisenstat trick

  Level: intermediate

.seealso: [](ch_ksp), `PCSOR`, `PCEisenstatSetOmega()`, `PCSORSetIterations()`, `PCSORSetOmega()`
@*/
PetscErrorCode PCSORSetSymmetric(PC pc, MatSORType flag)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscValidLogicalCollectiveEnum(pc, flag, 2);
  PetscTryMethod(pc, "PCSORSetSymmetric_C", (PC, MatSORType), (pc, flag));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCSORSetOmega - Sets the SOR relaxation coefficient, omega
  (where omega = 1.0 by default).

  Logically Collective

  Input Parameters:
+ pc    - the preconditioner context
- omega - relaxation coefficient (0 < omega < 2).

  Options Database Key:
. -pc_sor_omega <omega> - Sets omega

  Level: intermediate

  Note:
  If omega != 1, you will need to set the `MAT_USE_INODE`S option to `PETSC_FALSE` on the matrix.

.seealso: [](ch_ksp), `PCSOR`, `PCSORSetSymmetric()`, `PCSORSetIterations()`, `PCEisenstatSetOmega()`, `MatSetOption()`
@*/
PetscErrorCode PCSORSetOmega(PC pc, PetscReal omega)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscValidLogicalCollectiveReal(pc, omega, 2);
  PetscTryMethod(pc, "PCSORSetOmega_C", (PC, PetscReal), (pc, omega));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  PCSORSetIterations - Sets the number of inner iterations to
  be used by the SOR preconditioner. The default is 1.

  Logically Collective

  Input Parameters:
+ pc   - the preconditioner context
. lits - number of local iterations, smoothings over just variables on processor
- its  - number of parallel iterations to use; each parallel iteration has lits local iterations

  Options Database Keys:
+ -pc_sor_its <its>   - Sets number of iterations
- -pc_sor_lits <lits> - Sets number of local iterations

  Level: intermediate

  Note:
  When run on one processor the number of smoothings is lits*its

.seealso: [](ch_ksp), `PCSOR`, `PCSORSetOmega()`, `PCSORSetSymmetric()`
@*/
PetscErrorCode PCSORSetIterations(PC pc, PetscInt its, PetscInt lits)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscValidLogicalCollectiveInt(pc, its, 2);
  PetscTryMethod(pc, "PCSORSetIterations_C", (PC, PetscInt, PetscInt), (pc, its, lits));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*MC
     PCSOR - (S)SOR (successive over relaxation, Gauss-Seidel) preconditioning

   Options Database Keys:
+  -pc_sor_symmetric - Activates symmetric version
.  -pc_sor_backward - Activates backward version
.  -pc_sor_forward - Activates forward version
.  -pc_sor_local_forward - Activates local forward version
.  -pc_sor_local_symmetric - Activates local symmetric version  (default version)
.  -pc_sor_local_backward - Activates local backward version
.  -pc_sor_omega <omega> - Sets omega
.  -pc_sor_diagonal_shift <shift> - shift the diagonal entries; useful if the matrix has zeros on the diagonal
.  -pc_sor_its <its> - Sets number of iterations   (default 1)
-  -pc_sor_lits <lits> - Sets number of local iterations  (default 1)

   Level: beginner

   Notes:
   Only implemented for the `MATAIJ`  and `MATSEQBAIJ` matrix formats.

   Not a true parallel SOR, in parallel this implementation corresponds to block
   Jacobi with SOR on each block.

          For `MATAIJ` matrices if a diagonal entry is zero (and the diagonal shift is zero) then by default the inverse of that
          zero will be used and hence the `KSPSolve()` will terminate with `KSP_DIVERGED_NANORINF`. If the option
          `KSPSetErrorIfNotConverged()` or -ksp_error_if_not_converged the code will terminate as soon as it detects the
          zero pivot.

          For `MATSEQBAIJ` matrices this implements point-block SOR, but the omega, its, lits options are not supported.

          For `MATSEQBAIJ` the diagonal blocks are inverted using dense LU with partial pivoting. If a zero pivot is detected
          the computation is stopped with an error

          If used with `KSPRICHARDSON` and no monitors the convergence test is skipped to improve speed, thus it always iterates
          the maximum number of iterations you've selected for `KSP`. It is usually used in this mode as a smoother for multigrid.

          If omega != 1, you will need to set the `MAT_USE_INODES` option to `PETSC_FALSE` on the matrix.

.seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCJACOBI`,
          `PCSORSetIterations()`, `PCSORSetSymmetric()`, `PCSORSetOmega()`, `PCEISENSTAT`, `MatSetOption()`
M*/

PETSC_EXTERN PetscErrorCode PCCreate_SOR(PC pc)
{
  PC_SOR *jac;

  PetscFunctionBegin;
  PetscCall(PetscNew(&jac));

  pc->ops->apply           = PCApply_SOR;
  pc->ops->applytranspose  = PCApplyTranspose_SOR;
  pc->ops->applyrichardson = PCApplyRichardson_SOR;
  pc->ops->setfromoptions  = PCSetFromOptions_SOR;
  pc->ops->setup           = NULL;
  pc->ops->view            = PCView_SOR;
  pc->ops->destroy         = PCDestroy_SOR;
  pc->data                 = (void *)jac;
  jac->sym                 = SOR_LOCAL_SYMMETRIC_SWEEP;
  jac->omega               = 1.0;
  jac->fshift              = 0.0;
  jac->its                 = 1;
  jac->lits                = 1;

  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORSetSymmetric_C", PCSORSetSymmetric_SOR));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORSetOmega_C", PCSORSetOmega_SOR));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORSetIterations_C", PCSORSetIterations_SOR));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORGetSymmetric_C", PCSORGetSymmetric_SOR));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORGetOmega_C", PCSORGetOmega_SOR));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORGetIterations_C", PCSORGetIterations_SOR));
  PetscFunctionReturn(PETSC_SUCCESS);
}