static char help[] = "Tests Cholesky factorization for a SBAIJ matrix, (bs=2).\n";
/*
  This code is modified from the code contributed by JUNWANG@uwm.edu on Apr 13, 2007
*/

#include <petscmat.h>

int main(int argc, char **args)
{
  Mat           mat, fact, B;
  PetscInt      ind1[2], ind2[2];
  PetscScalar   temp[4];
  PetscInt      nnz[3];
  IS            perm, colp;
  MatFactorInfo info;
  PetscMPIInt   size;

  PetscFunctionBeginUser;
  PetscCall(PetscInitialize(&argc, &args, 0, help));
  PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
  PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");

  nnz[0] = 2;
  nnz[1] = 1;
  nnz[2] = 1;
  PetscCall(MatCreateSeqSBAIJ(PETSC_COMM_SELF, 2, 6, 6, 0, nnz, &mat));

  ind1[0] = 0;
  ind1[1] = 1;
  temp[0] = 3;
  temp[1] = 2;
  temp[2] = 0;
  temp[3] = 3;
  PetscCall(MatSetValues(mat, 2, ind1, 2, ind1, temp, INSERT_VALUES));
  ind2[0] = 4;
  ind2[1] = 5;
  temp[0] = 1;
  temp[1] = 1;
  temp[2] = 2;
  temp[3] = 1;
  PetscCall(MatSetValues(mat, 2, ind1, 2, ind2, temp, INSERT_VALUES));
  ind1[0] = 2;
  ind1[1] = 3;
  temp[0] = 4;
  temp[1] = 1;
  temp[2] = 1;
  temp[3] = 5;
  PetscCall(MatSetValues(mat, 2, ind1, 2, ind1, temp, INSERT_VALUES));
  ind1[0] = 4;
  ind1[1] = 5;
  temp[0] = 5;
  temp[1] = 1;
  temp[2] = 1;
  temp[3] = 6;
  PetscCall(MatSetValues(mat, 2, ind1, 2, ind1, temp, INSERT_VALUES));

  PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
  PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));

  PetscCall(MatDuplicate(mat, MAT_SHARE_NONZERO_PATTERN, &B));
  ind1[0] = 0;
  ind1[1] = 1;
  temp[0] = 3;
  temp[1] = 2;
  temp[2] = 0;
  temp[3] = 3;
  PetscCall(MatSetValues(mat, 2, ind1, 2, ind1, temp, INSERT_VALUES));
  ind2[0] = 4;
  ind2[1] = 5;
  temp[0] = 1;
  temp[1] = 1;
  temp[2] = 2;
  temp[3] = 1;
  PetscCall(MatSetValues(mat, 2, ind1, 2, ind2, temp, INSERT_VALUES));
  ind1[0] = 2;
  ind1[1] = 3;
  temp[0] = 4;
  temp[1] = 1;
  temp[2] = 1;
  temp[3] = 5;
  PetscCall(MatSetValues(mat, 2, ind1, 2, ind1, temp, INSERT_VALUES));
  ind1[0] = 4;
  ind1[1] = 5;
  temp[0] = 5;
  temp[1] = 1;
  temp[2] = 1;
  temp[3] = 6;
  PetscCall(MatSetValues(mat, 2, ind1, 2, ind1, temp, INSERT_VALUES));

  PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
  PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));

  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "mat: \n"));
  PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_SELF));

  /* begin cholesky factorization */
  PetscCall(MatGetOrdering(mat, MATORDERINGNATURAL, &perm, &colp));
  PetscCall(ISDestroy(&colp));

  info.fill = 1.0;
  PetscCall(MatGetFactor(mat, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &fact));
  PetscCall(MatCholeskyFactorSymbolic(fact, mat, perm, &info));
  PetscCall(MatCholeskyFactorNumeric(fact, mat, &info));

  PetscCall(ISDestroy(&perm));
  PetscCall(MatDestroy(&mat));
  PetscCall(MatDestroy(&fact));
  PetscCall(MatDestroy(&B));
  PetscCall(PetscFinalize());
  return 0;
}