static char help[] = "Test MatMult() for Hermitian matrix.\n\n";

#include <petscmat.h>

int main(int argc, char **args)
{
  Mat         A, As, B;
  Vec         l, r;
  PetscBool   flg;
  PetscMPIInt size;
  PetscInt    i, j;
  PetscScalar v, sigma2;
  PetscReal   h2, sigma1 = 100.0;
  PetscInt    dim, Ii, J, n = 3, rstart, rend;

  PetscFunctionBeginUser;
  PetscCall(PetscInitialize(&argc, &args, NULL, help));
  PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
  PetscCall(PetscOptionsGetReal(NULL, NULL, "-sigma1", &sigma1, NULL));
  PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
  dim = n * n;

  PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
  PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, dim, dim));
  PetscCall(MatSetType(A, MATAIJ));
  PetscCall(MatSetFromOptions(A));
  PetscCall(MatSetUp(A));

  sigma2 = 10.0 * PETSC_i;
  h2     = 1.0 / ((n + 1) * (n + 1));

  PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
  for (Ii = rstart; Ii < rend; Ii++) {
    v = -1.0;
    i = Ii / n;
    j = Ii - i * n;
    if (i > 0) {
      J = Ii - n;
      PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
    }
    if (i < n - 1) {
      J = Ii + n;
      PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
    }
    if (j > 0) {
      J = Ii - 1;
      PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
    }
    if (j < n - 1) {
      J = Ii + 1;
      PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
    }
    v = 4.0 - sigma1 * h2;
    PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, ADD_VALUES));
  }
  PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
  PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

  /* Check whether A is symmetric */
  PetscCall(PetscOptionsHasName(NULL, NULL, "-check_symmetric", &flg));
  if (flg) {
    PetscCall(MatIsSymmetric(A, 0.0, &flg));
    PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "A is not symmetric");
  }
  PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));

  /* make A complex Hermitian */
  Ii = 0;
  J  = dim - 1;
  if (Ii >= rstart && Ii < rend) {
    v = sigma2 * h2; /* RealPart(v) = 0.0 */
    PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
    v = -sigma2 * h2;
    PetscCall(MatSetValues(A, 1, &J, 1, &Ii, &v, ADD_VALUES));
  }

  Ii = dim - 2;
  J  = dim - 1;
  if (Ii >= rstart && Ii < rend) {
    v = sigma2 * h2; /* RealPart(v) = 0.0 */
    PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
    v = -sigma2 * h2;
    PetscCall(MatSetValues(A, 1, &J, 1, &Ii, &v, ADD_VALUES));
  }

  PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
  PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
  PetscCall(MatViewFromOptions(A, NULL, "-disp_mat"));

  if (flg) {
    PetscCall(MatIsSymmetric(A, 0.0, &flg));
    PetscCheck(!flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "A is symmetric");
  }
  PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_FALSE));
  /* Check whether A is Hermitian, then set A->hermitian flag */
  PetscCall(PetscOptionsHasName(NULL, NULL, "-check_Hermitian", &flg));
  if (flg && size == 1) {
    PetscCall(MatIsHermitian(A, 0.0, &flg));
    PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "A is not Hermitian");
  }
  PetscCall(MatSetOption(A, MAT_HERMITIAN, PETSC_TRUE));

#if defined(PETSC_HAVE_SUPERLU_DIST)
  /* Test Cholesky factorization */
  PetscCall(PetscOptionsHasName(NULL, NULL, "-test_choleskyfactor", &flg));
  if (flg) {
    Mat           F;
    IS            perm, iperm;
    MatFactorInfo info;
    PetscInt      nneg, nzero, npos;

    PetscCall(MatGetFactor(A, MATSOLVERSUPERLU_DIST, MAT_FACTOR_CHOLESKY, &F));
    PetscCall(MatGetOrdering(A, MATORDERINGND, &perm, &iperm));
    PetscCall(MatCholeskyFactorSymbolic(F, A, perm, &info));
    PetscCall(MatCholeskyFactorNumeric(F, A, &info));

    PetscCall(MatGetInertia(F, &nneg, &nzero, &npos));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n", nneg, nzero, npos));
    PetscCall(MatDestroy(&F));
    PetscCall(ISDestroy(&perm));
    PetscCall(ISDestroy(&iperm));
  }
#endif

  /* Create a Hermitian matrix As in sbaij format */
  PetscCall(MatConvert(A, MATSBAIJ, MAT_INITIAL_MATRIX, &As));
  PetscCall(MatViewFromOptions(As, NULL, "-disp_mat"));

  /* Test MatMult */
  PetscCall(MatMultEqual(A, As, 10, &flg));
  PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "MatMult not equal");
  PetscCall(MatMultAddEqual(A, As, 10, &flg));
  PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "MatMultAdd not equal");

  PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
  PetscCall(MatRealPart(B));
  PetscCall(MatSetOption(B, MAT_SYMMETRIC, PETSC_TRUE));

  PetscCall(MatCreateVecs(A, &r, &l));
  PetscCall(VecSetRandom(r, NULL));
  PetscCall(VecCopy(r, l));
  PetscCall(MatDiagonalScale(A, r, l));
  PetscCall(PetscOptionsHasName(NULL, NULL, "-check_Hermitian", &flg));
  if (flg && size == 1) {
    PetscCall(MatIsHermitian(A, 0.0, &flg));
    PetscCheck(!flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "A is Hermitian");
  }
  PetscCall(MatSetOption(A, MAT_HERMITIAN, PETSC_FALSE));

  /* Free spaces */
  PetscCall(MatDestroy(&A));
  PetscCall(MatDestroy(&As));
  PetscCall(MatDestroy(&B));
  PetscCall(VecDestroy(&r));
  PetscCall(VecDestroy(&l));
  PetscCall(PetscFinalize());
  return 0;
}

/*TEST

   build:
      requires: complex

   test:
      args: -n 1000 -check_symmetric -check_Hermitian
      output_file: output/empty.out

   test:
      suffix: 2
      nsize: 3
      args: -n 1000
      output_file: output/empty.out

   test:
      suffix: superlu_dist
      nsize: 3
      requires: superlu_dist
      args: -test_choleskyfactor -mat_superlu_dist_rowperm NOROWPERM
      output_file: output/ex127_superlu_dist.out
TEST*/