static char help[] = "Solves a linear system in parallel with KSP.\n\
Input parameters include:\n\
	-random_exact_sol : use a random exact solution vector\n\
	-view_exact_sol : write exact solution vector to stdout\n\
	-m <mesh_x> : number of mesh points in x-direction\n\
	-n <mesh_n> : number of mesh points in y-direction\n\n";

#include <petscksp.h>

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc, char **argv)
{
  int m, n, p, q, nz, i0;
  int id, core, local_r,fd;
  Mat A, Ap, Ap1;                /* linear system matrix */
  IS  rowperm, colperm;     /* row and column permutations */
  PetscBool flg = PETSC_FALSE;
  PetscScalar	v;
  PetscInt i, j, Istart, Iend;
  MatOrderingType ordering = MATORDERINGRCM;
  PetscLogStage stage;
  char buf[PETSC_MAX_PATH_LEN];
  FILE *fp;
  PetscErrorCode ierr;

  ierr = PetscInitialize(&argc, &argv,(char*)0,help);CHKERRQ(ierr);

    // Load COO format data
    fp = fopen(argv[1], "r");
    do fgets(buf,PETSC_MAX_PATH_LEN-1,fp);
    while (buf[0] == '%');
    sscanf(buf,"%d %d %d\n",&m,&n,&nz);
    ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
    ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr);
    ierr = MatSetFromOptions(A);CHKERRQ(ierr);
    ierr = MatSeqAIJSetPreallocation(A, PetscMax(m/2, 4), NULL);CHKERRQ(ierr);
    ierr = PetscLogStageRegister("Assembly", &stage);CHKERRQ(ierr);
    ierr = PetscLogStagePush(stage);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "m = %d, n = %d, nz = %d\n", m, n, nz);CHKERRQ(ierr);
    for (i = 0; i < nz; i++) {
      fscanf(fp,"%d %d %le\n",&p,&q,&v);
      p -= 1;q -= 1;
      ierr = MatSetValues(A, 1, &p, 1, &q, &v,ADD_VALUES);CHKERRQ(ierr);
    }
    fclose(fp);
    ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
	ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
	ierr = PetscLogStagePop();CHKERRQ(ierr);
    
    ierr = MatGetOrdering(A,ordering,&rowperm,&colperm);CHKERRQ(ierr);
    ierr = MatPermute(A,rowperm,colperm,&Ap);CHKERRQ(ierr);
    ierr = MatViewFromOptions(Ap, NULL, "-mat_view");CHKERRQ(ierr);

    ierr = ISDestroy(&rowperm);CHKERRQ(ierr);
    ierr = ISDestroy(&colperm);CHKERRQ(ierr);
    ierr = MatDestroy(&A);CHKERRQ(ierr);
    ierr = MatDestroy(&Ap);CHKERRQ(ierr);
	PetscFinalize();
	return 0;
}