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,vec,vec1;              /* initial vector */
  Mat            B,C, K,M;               /* operator matrix */
  EPS            eps;             /* eigenproblem solver context */
  const EPSType  type;
  PetscInt       N,m=15,nev,max_size,row,column,test_int;
  PetscErrorCode ierr;
  PetscScalar eigen_r,eigen_i;
  double value,test;
  int i,j, counter, counter2;
  FILE *myfile;

  SlepcInitialize(&argc,&argv,(char*)0,help);
  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Compute the operator matrix that defines the eigensystem, Ax=kx
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = MatCreate(PETSC_COMM_WORLD,&B);CHKERRQ(ierr);
  ierr = MatCreate(PETSC_COMM_WORLD,&C);CHKERRQ(ierr);
  ierr = MatCreate(PETSC_COMM_WORLD,&K);CHKERRQ(ierr);
  ierr = MatCreate(PETSC_COMM_WORLD,&M);CHKERRQ(ierr);

  ierr=  VecCreate(PETSC_COMM_WORLD,&vec); CHKERRQ(ierr);
  ierr=  VecCreate(PETSC_COMM_WORLD,&vec1); CHKERRQ(ierr);
// Initilize Matrices
  ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,48,48);CHKERRQ(ierr);
  ierr = MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,48,48);CHKERRQ(ierr);
  ierr = MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,48,48);CHKERRQ(ierr);
  ierr = MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,48,48);CHKERRQ(ierr);

  ierr = VecSetSizes(vec,PETSC_DECIDE,48);
  ierr = VecSetSizes(vec1,PETSC_DECIDE,48);
// Set up the matrices
  ierr = MatSetUp(B);CHKERRQ(ierr);
  ierr = MatSetUp(C);CHKERRQ(ierr);
  ierr = VecSetUp(vec); CHKERRQ(ierr);
  ierr = VecSetUp(vec1); CHKERRQ(ierr);

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

    counter = 0; counter2 = 0;
    myfile= fopen("mass.txt","r");
    if (myfile != NULL)
    {
        fscanf(myfile,"%d",&max_size);
        while ( fscanf(myfile,"%d",&row) !=EOF )
        {
          fscanf(myfile,"%d",&column);
          fscanf(myfile,"%lf",&value);
          MatSetValue(C,row-1,column-1,value,INSERT_VALUES);
//        MatSetValue(C,column-1,row-1,value,INSERT_VALUES);
        }
    }
  fclose(myfile);

  ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
  ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
  ierr = VecAssemblyBegin(vec);
  ierr = VecAssemblyBegin(vec1);

  ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
  ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
  ierr = VecAssemblyEnd(vec);
  ierr = VecAssemblyEnd(vec1);

  VecMax(vec,&i, &test);
  VecGetSize(vec,&i);
  ierr = PetscPrintf(PETSC_COMM_WORLD,": %f : %d",test,i);CHKERRQ(ierr);
  MatGetVecs(B,PETSC_NULL,&vec);
  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                Create the eigensolver and set various options
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  /*  Create eigensolver context  */
  ierr = EPSCreate(PETSC_COMM_WORLD,&eps);CHKERRQ(ierr);
  //ierr = EPSSetDimensions(eps,5,5,PETSC_NULL);
  //ierr = EPSSetProblemType(eps,EPS_NHEP);CHKERRQ(ierr);
  //ierr = EPSSetType(eps, EPSSUBSPACE);
  ierr = EPSSetOperators(eps,B,C);CHKERRQ(ierr);
  //EPSSetTolerances(eps,0.0001,300);
  EPSSetWhichEigenpairs(eps,EPS_SMALLEST_MAGNITUDE);
  //EPSSetInterval(eps,0,10000.0);
//  EPSSetLeftVectorsWanted(eps,PETSC_TRUE);
  EPSSetFromOptions(eps);
  /* - - - Solve the eigensystem - - - */
  ierr = EPSSolve(eps);CHKERRQ(ierr);

//POST PROCESSING
  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);
  //EPSGetEigenvalue(eps,0,&eigen_r,&eigen_i);
  EPSGetConverged(eps,&test_int);
  //EPSGetEigenpair(eps,5,&eigen_i,&eigen_r,vec,vec1);
  EPSGetEigenvector(eps,1,vec,vec1);
  PetscPrintf(MPI_COMM_WORLD,"%f",eigen_r);
  for (i=0;i<9;i++)
  {
    PetscPrintf(MPI_COMM_WORLD," %.8f:   %.8f:",&vec[i],vec1[i] );
  }

  ierr = EPSDestroy(&eps);CHKERRQ(ierr);
  ierr = VecDestroy(&vec);CHKERRQ(ierr);
  ierr = SlepcFinalize();CHKERRQ(ierr);
  return 0;
}