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 <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <petscksp.h>
#include "petscmat.h" 
#include <omp.h>
#include "mpi.h"
#include "mkl.h"

typedef double FLOAT_TYPE;

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **argv)
{
    int m, n, p, q, nz, i0;
    int id, core, local_r,fd;
	Vec x,b,u;                /* approx solution, RHS, exact solution */
	Mat A, Ap;                /* linear system matrix */
	KSP	ksp;                  /* linear solver context */
	PetscRandom	rctx;         /* random number generator context */
	PetscReal norm;           /* norm of solution error */
    IS  rowperm, colperm;     /* row and column permutations */
	PetscErrorCode ierr;
	PetscBool flg = PETSC_FALSE;
	PetscScalar	v;
    PetscInt i, j, Istart, Iend;
    MatOrderingType  ordering = MATORDERINGRCM;
    PetscViewer    view_out,view_in;
    #if defined(PETSC_USE_LOG)
    PetscLogStage stage;
    #endif
    char buf[PETSC_MAX_PATH_LEN];
    FILE *fp;
    //FLOAT_TYPE value;
	PetscInitialize(&argc, &argv,(char*)0,help);

    fp = fopen(argv[1], "r");
    // Process header with comments
    do fgets(buf,PETSC_MAX_PATH_LEN-1,fp);
        while (buf[0] == '%');
    // Get the initial parameters
    sscanf(buf,"%d %d %d\n",&m,&n,&nz);

    // Build Matrix
    MatCreate(PETSC_COMM_WORLD,&A);
    MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n);
    MatSetFromOptions(A);
    // Type in number of nonzeros per row (same for all rows)
    MatSeqAIJSetPreallocation(A,MAX(m/2,4),NULL);
    MatGetOwnershipRange(A,&Istart,&Iend);
    PetscLogStageRegister("Assembly", &stage);
    PetscLogStagePush(stage);
    printf("m = %d, n = %d, nz = %d\n", m, n, nz);
    // load COO format data
    for (i = 0; i < nz; i++) {
        fscanf(fp,"%d %d %le\n",&p,&q,&v);
        p -= 1;q -= 1;
        MatSetValues(A, 1, &p, 1, &q, &v,ADD_VALUES);
    }
    fclose(fp);
	
    MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
	MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
	PetscLogStagePop();
    
    // Start Reverse Cuthill Mckee
    MatGetOrdering(A,ordering,&rowperm,&colperm);

    MatPermute(A,rowperm,colperm,&Ap);
    //MatView(A,PETSC_VIEWER_STDOUT_SELF);
    //MatView(Ap, PETSC_VIEWER_STDOUT_WORLD);
    //ISView(rowperm,PETSC_VIEWER_STDOUT_SELF);
    //ISView(colperm,PETSC_VIEWER_STDOUT_SELF);

    // Open viewer for binary output
    PetscViewerBinaryOpen(PETSC_COMM_SELF,"./input.dat",FILE_MODE_WRITE,&view_out);
    PetscViewerBinaryGetDescriptor(view_out,&fd);
    // !!! Caution, this can't be annotated
    MatView(Ap,view_out);
    
    // Free work space Destroy the output viewer and work array
    PetscViewerDestroy(&view_out);
    MatDestroy(&A);
    MatDestroy(&Ap);
    
	PetscFinalize();
	return 0;
}