static char help[] = "Simple Hello World example program in SLEPc\n";

#include "slepceps.h"

#undef __FUNCT__
#define __FUNCT__ "main"
int main( int argc, char **argv )
{
    Vec            v0;              /* initial vector */
  Mat            A,B,C;               /* operator matrix */
  EPS            eps;             /* eigenproblem solver context */
  const EPSType  type;
  PetscInt       N,m=15,nev,max_size,row,column;
  PetscErrorCode ierr;
  double test,value;
  FILE *myfile;

  SlepcInitialize(&argc,&argv,(char*)0,help);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Compute the operator matrix that defines the eigensystem, Ax=kx
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatCreate(PETSC_COMM_WORLD,&B);CHKERRQ(ierr);
  ierr = MatCreate(PETSC_COMM_WORLD,&C);CHKERRQ(ierr);

  ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr);
  ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,96,96);CHKERRQ(ierr);
  ierr = MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,96,96);CHKERRQ(ierr);

//  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
//  ierr = MatSetUp(A);CHKERRQ(ierr);

  ierr = MatSetUp(B);CHKERRQ(ierr);
  MatSetValue(B,1,1,13.45,ADD_VALUES);

  ierr = MatSetUp(C);CHKERRQ(ierr);
  MatSetValue(C,1,1,13.45,ADD_VALUES);

    myfile= fopen("IN.txt","r");
    if (myfile != NULL)
    {
        fscanf(myfile,"%lf",&max_size);
        while ( fscanf(myfile,"%d",&row) !=EOF )
        {
          fscanf(myfile,"%d",&column);
          fscanf(myfile,"%lf",&value);
          MatSetValue(B,row-1,column-1,value,ADD_VALUES);
        }
    }
    fclose(myfile);

    myfile= fopen("mass.txt","r");
    if (myfile != NULL)
    {
        fscanf(myfile,"%lf",&max_size);
        while ( fscanf(myfile,"%d",&row) !=EOF )
        {
          fscanf(myfile,"%d",&column);
          fscanf(myfile,"%lf",&value);
          MatSetValue(C,row-1,column-1,value,ADD_VALUES);
        }
    }
    fclose(myfile);
  //ierr = MatMarkovModel(m,A);CHKERRQ(ierr);
  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                Create the eigensolver and set various options
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  /*  Create eigensolver context  */
  ierr = EPSCreate(PETSC_COMM_WORLD,&eps);CHKERRQ(ierr);
  /*
     Set operators. In this case, it is a standard eigenvalue problem
  */
  ierr = EPSSetType(eps,EPSSUBSPACE);
  ierr = EPSSetOperators(eps,B,C);CHKERRQ(ierr);
  ierr = EPSSetProblemType(eps,EPS_GNHEP);CHKERRQ(ierr);
  /*
     Set the initial vector. This is optional, if not done the initial
     vector is set to random values
  */
//  ierr = MatGetVecs(A,&v0,PETSC_NULL);CHKERRQ(ierr);
//  ierr = VecSet(v0,1.0);CHKERRQ(ierr);
//  ierr = EPSSetInitialSpace(eps,1,&v0);CHKERRQ(ierr);
  /* - - - Solve the eigensystem - - - */

  ierr = EPSSolve(eps);CHKERRQ(ierr);
  /*
     Optional: Get some information from the solver and display it
  */
  ierr = EPSGetType(eps,&type);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);CHKERRQ(ierr);
  ierr = EPSGetDimensions(eps,&nev,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);CHKERRQ(ierr);
  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                    Display solution and clean up
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = EPSPrintSolution(eps,PETSC_NULL);CHKERRQ(ierr);
  ierr = EPSDestroy(&eps);CHKERRQ(ierr);
  ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = VecDestroy(&v0);CHKERRQ(ierr);
  ierr = SlepcFinalize();CHKERRQ(ierr);
  return 0;
}