static char help[] = "Standard symmetric eigenproblem for the Laplacian in 2d and 3d with finite elements.\n\
We form the Laplacian for a rectangular domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\
Use -real_eigenvector_view draw -draw_pause 3 to visualize the eigenvectors\n\n\n";

#include <slepceps.h>
#include <petscdmplex.h>
#include <petscsnes.h>
#include <petscds.h>

/*
  If we want to construct a multilevel eigensolver for this problem, we start with a traditional EPS solver on the
coarse level. Then on the finer grids, we relax on the homogeneous equation
\begin{align*}
  \left( A_l - \lambda_l B \right) u_l = 0
\end{align*}
and then correct the approximate eigenvalue using
\begin{align*}
  \lamnda_l = \frac{\HC{u_l} A u_l}{\HC{u_l} B u_l}
\end{align*}
which looks like the inverse Rayleigh quotient iteration, but with inversion reduced to relaxation. In order to
justify this, we take the solution $u_{l+1}$ from the coarse level,
\begin{align*}
  A_{l+1} u_{l+1} = \lambda_{l+1} B_{l+1} u_{l+1}
\end{align*}
where the operators are defined by the Galerkin projection
\begin{align*}
  A_{l+1} &= \HC{P^l_{l+1}} A_l P^l_{l+1} \\
          &= \HC{P_{l+1}} A P_{l+1}
\end{align*}
where we have defined the composite interpolation operator
\begin{align*}
  P_{l+1} = P^0_1 P^1_2 \cdots P^{l-1}_l P^l_{l+1}
\end{align*}
which interpolates from level $l+1$ directly to the finest space. That coarse eigenapproximation will satisfy
\begin{align*}
  A_{l+1} u_{l+1} &= \lambda_{l+1} B_{l+1} u_{l+1} \\
  \HC{P^l_{l+1}} A_l P^l_{l+1} u_{l+1} &= \lambda_{l+1} \HC{P^l_{l+1}} B_l P^l_{l+1} u_{l+1}
  \HC{P^l_{l+1}} A_l {u_l}' &= \lambda_{l+1} \HC{P^l_{l+1}} B_l {u_l}'
\end{align*}
making ${u_l}' = P^l_{l+1} u_{l+1}$ a good initial guess for the generalized eigenproblem on level $l$ since it
satifies the equation for all modes resolvable by the interpolator.
*/

typedef struct {
  /* Domain and mesh definition */
  PetscInt  dim;               /* The topological mesh dimension */
  PetscBool simplex;           /* Simplicial mesh */
  PetscInt  cells[3];          /* The initial domain division */
} AppCtx;

typedef struct {
  Mat A; /* The system operator */
  PC  M; /* The smoother */
  Vec w; /* Work vector */
} OpCtx;

static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
{
  *u = 0.0;
  return 0;
}

static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
{
  PetscInt d;
  *u = 0.0;
  for (d = 0; d < dim; ++d) *u += PetscSinReal(2.0*PETSC_PI*x[d]);
  return 0;
}

static void f0_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                      PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
{
  PetscInt d;
  for (d = 0; d < dim; ++d) f0[0] += -4.0*PetscSqr(PETSC_PI)*PetscSinReal(2.0*PETSC_PI*x[d]);
}

static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
{
  PetscInt d;
  for (d = 0; d < dim; ++d) f1[d] = u_x[d];
}

static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
                  const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
                  const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
                  PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
{
  PetscInt d;
  for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0;
}

static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
{
  PetscInt       n = 3;
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  options->dim      = 2;
  options->cells[0] = 1;
  options->cells[1] = 1;
  options->cells[2] = 1;
  options->simplex  = PETSC_TRUE;

  ierr = PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");CHKERRQ(ierr);
  ierr = PetscOptionsInt("-dim", "The topological mesh dimension", "ex13.c", options->dim, &options->dim, NULL);CHKERRQ(ierr);
  ierr = PetscOptionsIntArray("-cells", "The initial mesh division", "ex13.c", options->cells, &n, NULL);CHKERRQ(ierr);
  ierr = PetscOptionsBool("-simplex", "Simplicial (true) or tensor (false) mesh", "ex13.c", options->simplex, &options->simplex, NULL);CHKERRQ(ierr);
  ierr = PetscOptionsEnd();
  PetscFunctionReturn(0);
}

static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
{
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  /* Create box mesh */
  ierr = DMPlexCreateBoxMesh(comm, user->dim, user->simplex, user->cells, NULL, NULL, NULL, PETSC_TRUE, dm);CHKERRQ(ierr);
  /* Distribute mesh over processes */
  {
    DM               dmDist = NULL;
    PetscPartitioner part;

    ierr = DMPlexGetPartitioner(*dm, &part);CHKERRQ(ierr);
    ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr);
    ierr = DMPlexDistribute(*dm, 0, NULL, &dmDist);CHKERRQ(ierr);
    if (dmDist) {
      ierr = DMDestroy(dm);CHKERRQ(ierr);
      *dm  = dmDist;
    }
  }
  /* TODO: This should be pulled into the library */
  {
    char      convType[256];
    PetscBool flg;

    ierr = PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");CHKERRQ(ierr);
    ierr = PetscOptionsFList("-dm_plex_convert_type","Convert DMPlex to another format","ex12",DMList,DMPLEX,convType,256,&flg);CHKERRQ(ierr);
    ierr = PetscOptionsEnd();
    if (flg) {
      DM dmConv;

      ierr = DMConvert(*dm,convType,&dmConv);CHKERRQ(ierr);
      if (dmConv) {
        ierr = DMDestroy(dm);CHKERRQ(ierr);
        *dm  = dmConv;
      }
    }
  }
  /* TODO: This should be pulled into the library */
  ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr);

  ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr);
  ierr = DMSetApplicationContext(*dm, user);CHKERRQ(ierr);
  ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);
  ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
{
  PetscDS        prob;
  const PetscInt id = 1;
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
  ierr = PetscDSSetResidual(prob, 0, f0_trig_u, f1_u);CHKERRQ(ierr);
  ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr);
  ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "wall", "marker", 0, 0, NULL, (void (*)(void)) trig_u, 1, &id, user);CHKERRQ(ierr);
  ierr = PetscDSSetExactSolution(prob, 0, trig_u);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user)
{
  DM             cdm = dm;
  PetscFE        fe;
  char           prefix[PETSC_MAX_PATH_LEN];
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  /* Create finite element */
  ierr = PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name);CHKERRQ(ierr);
  ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) dm), user->dim, 1, user->simplex, name ? prefix : NULL, -1, &fe);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) fe, name);CHKERRQ(ierr);
  /* Set discretization and boundary conditions for each mesh */
  ierr = DMSetField(dm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr);
  ierr = DMCreateDS(dm);CHKERRQ(ierr);
  ierr = (*setup)(dm, user);CHKERRQ(ierr);
  while (cdm) {
    ierr = DMCopyDisc(dm, cdm);CHKERRQ(ierr);
    /* TODO: Check whether the boundary of coarse meshes is marked */
    ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr);
  }
  ierr = PetscFEDestroy(&fe);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

static PetscErrorCode identity(Mat A, Vec x, Vec y)
{
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  ierr = VecCopy(x, y);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

static PetscErrorCode lufactor(Mat A, IS row, IS col, MatFactorInfo *info)
{
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  ierr = MatSetFactorType(A, MAT_FACTOR_LU);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

static PetscErrorCode factor(Mat A, MatFactorType ftype, Mat *F)
{
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  ierr = PetscObjectReference((PetscObject) A);CHKERRQ(ierr);
  *F   = A;
  PetscFunctionReturn(0);
}

static PetscErrorCode MinvA(Mat A, Vec x, Vec y)
{
  OpCtx         *ctx;
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  ierr = MatShellGetContext(A, &ctx);CHKERRQ(ierr);
  ierr = MatMult(ctx->A, x, ctx->w);CHKERRQ(ierr);
  ierr = PCApply(ctx->M, ctx->w, y);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

int main(int argc, char **argv)
{
  DM             dm;   /* Problem specification */
  EPS            eps;  /* Eigensolver */
  Mat            A, S; /* Operator and smoother */
  Vec            u;
  OpCtx          ctx;  /* Context for shell matmult */
  AppCtx         user; /* User-defined work context */
  PetscErrorCode ierr;

  ierr = SlepcInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
  ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
  ierr = EPSCreate(PETSC_COMM_WORLD, &eps);CHKERRQ(ierr);
  ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr);
  ierr = SetupDiscretization(dm, "potential", SetupPrimalProblem, &user);CHKERRQ(ierr);
  ierr = DMCreateMatrix(dm, &A);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) A, "Laplacian");CHKERRQ(ierr);
  ierr = DMGetGlobalVector(dm, &u);CHKERRQ(ierr);
  ierr = DMPlexSNESComputeJacobianFEM(dm, u, A, A, &user);CHKERRQ(ierr);
  ierr = DMRestoreGlobalVector(dm, &u);CHKERRQ(ierr);
  ctx.A = NULL;
  ctx.M = NULL;
  ctx.w = NULL;
#if 0
  /* Make a shell matrix for the smoother
       -eps_gen_hermitian -eps_smallest_magnitude -st_pc_factor_mat_ordering_type natural -st_pc_factor_in_place
  */
  {
    PetscInt m, n, M, N;

    ierr = MatCreate(PETSC_COMM_WORLD, &S);CHKERRQ(ierr);
    ierr = MatGetLocalSize(A, &m, &n);CHKERRQ(ierr);
    ierr = MatGetSize(A, &M, &N);CHKERRQ(ierr);
    ierr = MatSetSizes(S, m, n, M, N);CHKERRQ(ierr);
    ierr = MatSetType(S, MATSHELL);CHKERRQ(ierr);
    ierr = MatShellSetOperation(S, MATOP_MULT,     (void(*)(void)) identity);CHKERRQ(ierr);
    ierr = MatShellSetOperation(S, MATOP_LUFACTOR, (void(*)(void)) lufactor);CHKERRQ(ierr);
    ierr = MatShellSetOperation(S, MATOP_SOLVE,    (void(*)(void)) identity);CHKERRQ(ierr);
    ierr = MatSetUp(S);CHKERRQ(ierr);

    ierr = MatSolverTypeRegister("petsc", MATSHELL, MAT_FACTOR_LU, factor);CHKERRQ(ierr);
  }
  ierr = EPSSetOperators(eps, A, S);CHKERRQ(ierr);
#endif
#if 1
  /* Make shell matrix for the preconditioned operator
       -eps_non_hermitian -eps_smallest_magnitude -minva_pc_type sor -minva_pc_sor_its 4
  */
  {
    PetscInt m, n, M, N;

    ierr = PCCreate(PETSC_COMM_WORLD, &ctx.M);CHKERRQ(ierr);
    ierr = PCSetOptionsPrefix(ctx.M, "minva_");CHKERRQ(ierr);
    ierr = PCSetOperators(ctx.M, A, A);CHKERRQ(ierr);
    ierr = PCSetFromOptions(ctx.M);CHKERRQ(ierr);
    ierr = PCSetUp(ctx.M);CHKERRQ(ierr);
    ierr = MatCreateVecs(A, NULL, &ctx.w);CHKERRQ(ierr);
    ierr = PetscObjectReference((PetscObject) A);CHKERRQ(ierr);
    ctx.A = A;

    ierr = MatCreate(PETSC_COMM_WORLD, &S);CHKERRQ(ierr);
    ierr = MatGetLocalSize(A, &m, &n);CHKERRQ(ierr);
    ierr = MatGetSize(A, &M, &N);CHKERRQ(ierr);
    ierr = MatSetSizes(S, m, n, M, N);CHKERRQ(ierr);
    ierr = MatSetType(S, MATSHELL);CHKERRQ(ierr);
    ierr = MatShellSetOperation(S, MATOP_MULT, (void(*)(void)) MinvA);CHKERRQ(ierr);
    ierr = MatShellSetContext(S, &ctx);CHKERRQ(ierr);
    ierr = MatSetUp(S);CHKERRQ(ierr);

    ierr = MatSolverTypeRegister("petsc", MATSHELL, MAT_FACTOR_LU, factor);CHKERRQ(ierr);
  }
  ierr = EPSSetOperators(eps, S, NULL);CHKERRQ(ierr);
#endif
  ierr = EPSSetFromOptions(eps);CHKERRQ(ierr);
  ierr = EPSSolve(eps);CHKERRQ(ierr);
  {
    PetscDS   prob;
    Vec       xr, xi;
    PetscReal kr, ki;
    PetscInt  ncv, n;

    /* Change boundary conditions to 0 for the eigenvectors */
    ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
    ierr = PetscDSUpdateBoundary(prob, 0, DM_BC_ESSENTIAL, NULL, NULL, PETSC_DETERMINE, PETSC_DETERMINE, NULL, (void (*)(void)) zero, PETSC_DETERMINE, NULL, NULL);CHKERRQ(ierr);

    ierr = EPSErrorView(eps, EPS_ERROR_RELATIVE, NULL);CHKERRQ(ierr);
    ierr = EPSGetConverged(eps, &ncv);CHKERRQ(ierr);
    ierr = DMGetGlobalVector(dm, &xr);CHKERRQ(ierr);
    ierr = DMGetGlobalVector(dm, &xi);CHKERRQ(ierr);
    for (n = 0; n < ncv; ++n) {
      ierr = PetscObjectSetName((PetscObject) xr, "Real Eigenvector");CHKERRQ(ierr);
      ierr = PetscObjectSetName((PetscObject) xi, "Imaginary Eigenvector");CHKERRQ(ierr);
      ierr = EPSGetEigenpair(eps, n, &kr, &ki, xr, xi);CHKERRQ(ierr);
      ierr = PetscPrintf(PETSC_COMM_WORLD, "Eigenvector %D for value %.4g + %.4g i\n", n, (double) kr, (double) ki);CHKERRQ(ierr);
      ierr = VecViewFromOptions(xr, NULL, "-real_eigenvector_view");CHKERRQ(ierr);
      ierr = VecViewFromOptions(xi, NULL, "-imag_eigenvector_view");CHKERRQ(ierr);
    }
    ierr = DMRestoreGlobalVector(dm, &xr);CHKERRQ(ierr);
    ierr = DMRestoreGlobalVector(dm, &xi);CHKERRQ(ierr);
  }

  /* Cleanup */
  ierr = MatDestroy(&ctx.A);CHKERRQ(ierr);
  ierr = PCDestroy(&ctx.M);CHKERRQ(ierr);
  ierr = VecDestroy(&ctx.w);CHKERRQ(ierr);
  ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = MatDestroy(&S);CHKERRQ(ierr);
  ierr = EPSDestroy(&eps);CHKERRQ(ierr);
  ierr = DMDestroy(&dm);CHKERRQ(ierr);
  ierr = SlepcFinalize();
  return ierr;
}

/*TEST

  test:
    suffix: 2d_p1_0
    requires: triangle
    args: -potential_petscspace_degree 1 -dm_refine 2 -eps_hermitian -eps_smallest_magnitude

TEST*/