#define _GNU_SOURCE
#include <features.h>
#include <sys/types.h>
#include <stdlib.h>
#include <string.h>
#include <sys/types.h>
#include <stdio.h>
#include <math.h>
#include <ctype.h>
#include <unistd.h>
#include <sys/times.h>
#include <liboct_parser.h>
#include <mpi.h>
#include <slepceps.h>
#include <petscmat.h>

typedef struct _eigliststruct
{
  int index;
  double proj;
} eigliststruct;


int main(int argc, char *argv[])
{
  int ierr = 0;
  int myrank,nprocs;
  int findnextnonblankline(FILE *inputstream, size_t *inputsize, 
			   char **inputstring, int *offset);
  ierr = MPI_Init(&argc, &argv);

  static char slepchelp[] = "Dummy string to make slepc happy";
  SlepcInitialize(&argc, &argv, (char *) NULL, slepchelp);
  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
  MPI_Comm_size(MPI_COMM_WORLD, &nprocs);

  FILE *inputstream = NULL;
  ssize_t linelength = 0;
  size_t inputsize = 0;
  char *inputstring = NULL;
  int offset = 0;

  int nhamlines = 0;
  int maxi = -1;
  int *ind1 = NULL;
  int *ind2 = NULL;
  double *hval = NULL;
  int number = 0;
  int *targetset = NULL;
  int i,j,k;
  
  /* do the reading on proc 0 */
  if (myrank == 0) {
    inputstream = fopen("slepctestham.dat", "r");
    if (inputstream==NULL) {
      puts("error opening slepctestham.dat");
      exit(1);
    }
   linelength = findnextnonblankline(inputstream, &inputsize, &inputstring, 
				      &offset);
    /* read line by line */
    while (linelength>0) {
      ind1 = realloc(ind1, (nhamlines+1)*sizeof(int));
      ind2 = realloc(ind2, (nhamlines+1)*sizeof(int));
      hval = realloc(hval, (nhamlines+1)*sizeof(double));
      int nread = sscanf(inputstring, "HAM %d %d %lf ", ind1+nhamlines, ind2+nhamlines,
		     hval+nhamlines);
      if (nread==3) {
	if (ind1[nhamlines] > maxi) maxi = ind1[nhamlines];
	if (ind2[nhamlines] > maxi) maxi = ind2[nhamlines];
	nhamlines++;
      }
      linelength = findnextnonblankline(inputstream, &inputsize, &inputstring, 
					&offset);
    }
    /* I finished reading in the dummy hamiltonian */
    fclose(inputstream);
    inputstream = NULL;

    inputstream = fopen("slepctesttargets.dat", "r");
    if (inputstream==NULL) {
      puts("error opening slepctesttargets.dat");
      exit(1);
    }
    number = 1;
    targetset = calloc(1, sizeof(int));
    linelength = findnextnonblankline(inputstream, &inputsize, &inputstring, 
				      &offset);
    /* read line by line */
    while (linelength>0) {
      targetset = realloc(targetset, (number+1)*sizeof(int));
      int nread = sscanf(inputstring, " %d ", targetset+number);
      if (nread==1) {
	number++;
      }
      linelength = findnextnonblankline(inputstream, &inputsize, &inputstring, 
					&offset);
    }
    /* I finished reading in the dummy target list */
    fclose(inputstream);
    inputstream = NULL;
    free(inputstring);
    inputstring = NULL;
    targetset[0] = number-1;

    /* Now that I read the data, propagate them to other procs */
    /* Broadcast the number of vib modes to all nodes */
    ierr = MPI_Bcast(&nhamlines, 1, MPI_INT, 0, MPI_COMM_WORLD);
    ierr = MPI_Bcast(&maxi, 1, MPI_INT, 0, MPI_COMM_WORLD);
    ierr = MPI_Bcast(ind1, nhamlines, MPI_INT, 0, MPI_COMM_WORLD);
    ierr = MPI_Bcast(ind2, nhamlines, MPI_INT, 0, MPI_COMM_WORLD);
    ierr = MPI_Bcast(hval, nhamlines, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    ierr = MPI_Bcast(&number, 1, MPI_INT, 0, MPI_COMM_WORLD);
    ierr = MPI_Bcast(targetset, number, MPI_INT, 0, MPI_COMM_WORLD);
  }
  else {
    /* receive data from node 0 */
    ierr = MPI_Bcast(&nhamlines, 1, MPI_INT, 0, MPI_COMM_WORLD);
    ind1 = calloc(nhamlines, sizeof(int));
    ind2 = calloc(nhamlines, sizeof(int));
    hval = calloc(nhamlines, sizeof(double));
    ierr = MPI_Bcast(&maxi, 1, MPI_INT, 0, MPI_COMM_WORLD);
    ierr = MPI_Bcast(ind1, nhamlines, MPI_INT, 0, MPI_COMM_WORLD);
    ierr = MPI_Bcast(ind2, nhamlines, MPI_INT, 0, MPI_COMM_WORLD);
    ierr = MPI_Bcast(hval, nhamlines, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    ierr = MPI_Bcast(&number, 1, MPI_INT, 0, MPI_COMM_WORLD);
    targetset = calloc(number, sizeof(int));
    ierr = MPI_Bcast(targetset, number, MPI_INT, 0, MPI_COMM_WORLD);
  }   
    
  /* now I can create and build the dummy H matrix */
  Mat H;
  EPS eps;
  int nelems = maxi+1;
  ierr = MatCreate(MPI_COMM_WORLD, &H); CHKERRQ(ierr);
  ierr = MatSetSizes(H, PETSC_DECIDE, PETSC_DECIDE, nelems, nelems);
  CHKERRQ(ierr);
  ierr = MatSetFromOptions(H);CHKERRQ(ierr);
  ierr = MatSetUp(H);CHKERRQ(ierr);
  /* fill the dummy H matrix with the values I read */
  int first_row, last_row;
  MatGetOwnershipRange(H, &first_row, &last_row);
  last_row--;
  for (i=0; i<=nhamlines-1; i++) {
    if ((ind1[i]>=first_row)&&(ind1[i]<=last_row)) {
      MatSetValue(H, ind1[i], ind2[i], hval[i], ADD_VALUES);
    }
  }
  /* Assemble the hamiltonian H */
  ierr = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  if (ind1!=NULL) free(ind1);
  if (ind2!=NULL) free(ind2);
  if (hval!=NULL) free(hval);
		   
  PetscInt nconverged=0;
  int goodconverged=0;
  PetscScalar lambdar, lambdai;
  PetscReal tol=1000*PETSC_MACHINE_EPSILON;
  IS is;
  Vec localx, xr, xi;
  VecScatter scatter;
  PetscScalar *veccomps=NULL;
  int newtotconverged=0;
  
  double vectolerance = 0.1;
  int totconverged = 0;
  int maxmpd = 1000;
  
  ierr = MatCreateVecs(H, &xr, PETSC_NULL);CHKERRQ(ierr);
  ierr = MatCreateVecs(H, &xi, PETSC_NULL);CHKERRQ(ierr);
  ierr = ISCreateStride(MPI_COMM_WORLD,nelems,0,1,&is);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF, nelems, &localx);CHKERRQ(ierr);
  ierr = VecScatterCreate(xr,is,localx,is,&scatter); CHKERRQ(ierr);
  ierr = EPSCreate(MPI_COMM_WORLD,&eps);CHKERRQ(ierr);
  ierr = EPSSetOperators(eps,H,PETSC_NULL);CHKERRQ(ierr);
  ierr = EPSSetProblemType(eps,EPS_HEP);CHKERRQ(ierr);
  ierr = EPSSetTolerances(eps, tol, PETSC_DECIDE);
  ierr = EPSSetFromOptions(eps); CHKERRQ(ierr);
  ierr = PetscOptionsSetFromOptions(); CHKERRQ(ierr);
  const char common_options[] = "-st_ksp_type preonly \
                                       -st_pc_type cholesky";
  ierr = PetscOptionsInsertString(common_options); CHKERRQ(ierr);
  if (nprocs>1) {
    const char mumps_options[] = "-st_pc_factor_mat_solver_package mumps \
                                        -mat_mumps_icntl_13 1";
    ierr = PetscOptionsInsertString(mumps_options); CHKERRQ(ierr);
  }
  ierr = PetscOptionsInsert(&argc,&argv,PETSC_NULL); CHKERRQ(ierr);
  ierr = EPSSetFromOptions(eps); CHKERRQ(ierr);
  ierr = PetscOptionsSetFromOptions(); CHKERRQ(ierr);
  /* every time I solve, I want to find (at least) one eigenvector, and it must be
     the one with the largest component along the target state space */
  ierr = EPSSetWhichEigenpairs(eps, EPS_LARGEST_MAGNITUDE);
  PetscErrorCode computeprojection(PetscScalar er, PetscScalar ei, 
				   Vec xr, Vec xi, PetscScalar *rr,
				   PetscScalar *ri, void *ctx);
  ierr = EPSSetArbitrarySelection(eps, computeprojection, 
				  (void *) targetset);
  Vec *Cv = NULL;
  /* since it is apparently impossible to realloc with petsc, I allocate
     Cv as large as it can possibly need to be */
  PetscMalloc((nelems)*sizeof(Vec *), &Cv);
  Vec Iv;
  ierr = VecDuplicate(xr, &Iv);
  Vec yr;
  ierr = VecDuplicate(xr, &yr);
  double minconvergedcomponent = 0;
  int newminconvergedindex = 0;
  int minconvergedindex = -1;
  double *convergedcomponents = calloc(targetset[0], sizeof(double));
  for (i=0; i<=targetset[0]-1; i++) {
    convergedcomponents[i] = 0;
  }
  PetscInt ncv = PETSC_DECIDE;
  PetscInt mpd = PETSC_DECIDE;
  PetscInt nev = 1;
  ierr = EPSSetDimensions(eps, nev, ncv, mpd); CHKERRQ(ierr);
  int stopflag = 0;
  while (((1-minconvergedcomponent) >= vectolerance)&&(stopflag==0)) {
    if (newminconvergedindex != minconvergedindex) {
      minconvergedindex = newminconvergedindex;
      /* the initial guess coincides with the target vector indexed by minconvergedindex */
      ierr = VecZeroEntries(Iv);
      ierr = VecSetValue(Iv, (PetscInt) targetset[minconvergedindex+1], (PetscScalar) 1, 
			 INSERT_VALUES);
      /* I tell the solver to start from the target space, set before */
      ierr = VecAssemblyBegin(Iv);
      ierr = VecAssemblyEnd(Iv);
      ierr = EPSSetInitialSpace(eps, 1, &Iv);
    }
    if (totconverged>0) {
      /* I tell the solver to avoid the eigenstates already found before */
      ierr = EPSSetDeflationSpace(eps, totconverged, Cv);
    }

    nconverged = 0;
    while ((nconverged==0) && (ncv<=(nelems)) && (mpd<=maxmpd) &&
	   (mpd<=(nelems))) {
      ierr = EPSSolve(eps);CHKERRQ(ierr);
      /* How many eigenpairs did I get? */
      ierr = EPSGetConverged(eps, &nconverged); CHKERRQ(ierr);
      if (nconverged==0) {
	/* no solutions found, try harder */
	int newnev;
	if ((newnev = nev*4/3)> nev) nev = newnev;
	else nev++;
	printf("setting EPSSetDimensions nev=%d, ncv=%d, mpd=%d\n", nev, ncv, mpd);
	ierr = EPSSetDimensions(eps, nev, ncv, mpd);
	CHKERRQ(ierr);
      }
    }
    /* the eigenvectors should be sorted in decreasing component along
       the target state (that's what I asked at least, but it does not appear to quite work) */
    if (nconverged == 0) {
      /* EPSSolve did not get any solution even with the largest acceptable EPS size 
	 parameters, stop it here and exit the while loop */
      stopflag = 1;
    }
    else {
      goodconverged = nconverged;
      for (i=0; (i<=nconverged-1)&&((1-minconvergedcomponent)>=vectolerance); i++) {
	EPSGetEigenpair(eps, i, &lambdar, &lambdai, xr, xi);
	/* Check the component of this eigenvector on the target subspace */
	PetscScalar projr = 0;
	PetscScalar proji = 0;
	ierr = computeprojection(lambdar, lambdai, xr, xi, &projr, &proji, (void *) targetset);

	/* Check the component of this eigenvector on the deflation space 
	   (should be null by definition!) */
	int isdup = 0;
	PetscScalar defproj = 0;
	PetscScalar maxproj = 0;
	int maxprojindex;
	for (j=0; j<=totconverged-1; j++) {
	  PetscScalar dotprod = 0;
	  PetscScalar myproj = 0;
	  ierr = VecDot(xr, Cv[j], &dotprod);
	  defproj += (myproj = dotprod*dotprod);
	  if (myproj>maxproj) {
	    maxproj = myproj;
	    maxprojindex = j;
	  }
	}
	if (defproj > 1e-8) {
	  /* complain loudly, this should never happen */
	  printf("large projection %f of eigenstate %d on eigenstate %d\n",
		 maxproj, totconverged, maxprojindex);
	  isdup = 1;
	}
	if (isdup ==0) {
	  /* this is not a duplicated eigenvector, add it to the set */
  	  ierr = VecScatterBegin(scatter,xr,localx,INSERT_VALUES,
				 SCATTER_FORWARD);
	  ierr = VecScatterEnd(scatter,xr,localx,INSERT_VALUES,
			       SCATTER_FORWARD);
	  ierr = VecGetArray(localx,&veccomps);
	  
	  for (j=0; j<=targetset[0]-1; j++) {
	    convergedcomponents[j] += veccomps[targetset[j+1]]
	      *veccomps[targetset[j+1]];
	  }
	  minconvergedcomponent = convergedcomponents[0];
	  newminconvergedindex = 0;
	  for (j=1; j<=targetset[0]-1; j++) {
	    if (minconvergedcomponent>convergedcomponents[j]) {
	      minconvergedcomponent = convergedcomponents[j];
	      newminconvergedindex = j;
	    }
	  }
	  newtotconverged = totconverged+1;
	  ierr = VecRestoreArray(localx, &veccomps);
	  /* now add this eigenvector to the deflation space for next 
	     iteration */
	  ierr = VecDuplicate(xr, Cv+totconverged);
	  ierr = VecCopy(xr, Cv[totconverged]);
	  totconverged = newtotconverged;
	}
	else {
	  /* this was already in the deflation space to some extent, count it out */
	  goodconverged --;
	}
      }
      if (goodconverged==0) {
	/* All new eigenvectors were already in the deflation space! 
	   Try harder to get new ones*/
	int newnev;
	if ((newnev = nev*4/3)>nev) nev = newnev;
	else nev++;
	printf("setting EPSSetDimensions nev=%d, ncv=%d, mpd=%d\n", nev, ncv, mpd);
      	ierr = EPSSetDimensions(eps, nev, ncv, mpd);
      	CHKERRQ(ierr);
      }
      else {
	if (nconverged>1) { 
	  /* I found more eigenvectors than really necessary, reduce nev */
	  if ((nev=nev/2)<1) nev=1;
	  printf("setting EPSSetDimensions nev=%d, ncv=%d, mpd=%d\n", nev, ncv, mpd);
	  ierr = EPSSetDimensions(eps, nev, ncv, mpd);
	  CHKERRQ(ierr);
	}
      }
    }
  }

  /* Ok, now I can do some checks on whatever I may want to check before freeing everithing*/
  if (stopflag!=0) {
    printf("The eigensolver could not find enough solutions to cover the target space\n"); 
  }
  printf("final values of EPSSetDimensions nev=%d, ncv=%d, mpd=%d\n", nev, ncv, mpd);

  /* sort found eigenstates for projection on target space */
  eigliststruct *kseiglist = calloc(totconverged, sizeof(eigliststruct));
  for (i=0; i<=totconverged-1; i++) {
    kseiglist[i].index = i;
    PetscScalar projr = 0;
    PetscScalar proji = 0;
    ierr = computeprojection(0, 0, Cv[i], xi, &projr, &proji, (void *) targetset);
    kseiglist[i].proj = projr;
  }
  int eiglistcomp();
  qsort(kseiglist, totconverged, sizeof(eigliststruct), eiglistcomp);
  /* see how far down I have to go in the list to get the minimum coverage of target space 
     I want */
  for (j=0; j<=targetset[0]-1; j++) {
    convergedcomponents[j] = 0;
  }
  minconvergedcomponent = 0;
  for (i=0; (i<=totconverged-1) && ((1-minconvergedcomponent) >= vectolerance); i++) {
    minconvergedcomponent = 1;
    for (j=0; j<=targetset[0]-1; j++) {
      PetscScalar dotprod = 0;
      ierr = VecZeroEntries(yr);
      ierr = VecSetValue(yr, (PetscInt) targetset[j+1], (PetscScalar) 1,
			 INSERT_VALUES);
      ierr = VecDot(yr, Cv[kseiglist[i].index], &dotprod);
      convergedcomponents[j] += dotprod*dotprod;	
      if (minconvergedcomponent>convergedcomponents[j]) {
	minconvergedcomponent = convergedcomponents[j];
	minconvergedindex = j;
      }
    }
  }

  int needconverged = i;
  printf("Krylov-Schur run results:");
  printf("%d eigenstates out of %d determined, only %d actually needed\n",
	 totconverged, nelems, needconverged);

  /* Now redo everything with Lapack */
  /* Save K-S results I want to compare later */
  int kstotconverged = totconverged;
  Vec *ksCv = Cv;
  PetscMalloc((nelems)*sizeof(Vec *), &Cv);
  /* destroy the old eps, to recreate it anew */
  ierr = EPSDestroy(&eps);CHKERRQ(ierr);
  ierr = EPSCreate(MPI_COMM_WORLD,&eps);CHKERRQ(ierr);
  ierr = EPSSetOperators(eps,H,PETSC_NULL);CHKERRQ(ierr);
  ierr = EPSSetProblemType(eps,EPS_HEP);CHKERRQ(ierr);
  ierr = EPSSetTolerances(eps, tol, PETSC_DECIDE);
  ierr = EPSSetFromOptions(eps); CHKERRQ(ierr);
  ierr = PetscOptionsSetFromOptions(); CHKERRQ(ierr);
  ierr = EPSSetType(eps, EPSLAPACK); CHKERRQ(ierr);
  ierr = EPSSetFromOptions(eps); CHKERRQ(ierr);
  ierr = PetscOptionsSetFromOptions(); CHKERRQ(ierr);
  totconverged = 0;
  minconvergedcomponent = 0;
  newminconvergedindex = 0;
  minconvergedindex = -1;
  for (i=0; i<=targetset[0]-1; i++) {
    convergedcomponents[i] = 0;
  }
  nev = 1;
  ierr = EPSSolve(eps);CHKERRQ(ierr);
  /* How many eigenpairs did I get? */
  ierr = EPSGetConverged(eps, &nconverged); CHKERRQ(ierr);
  /* the eigenvectors should be sorted in decreasing component along
     the target state (that's what I asked at least, but it does not appear to quite work) */
  goodconverged = nconverged;
  for (i=0; (i<=nconverged-1)&&((1-minconvergedcomponent)>=vectolerance); i++) {
    EPSGetEigenpair(eps, i, &lambdar, &lambdai, xr, xi);
    /* Check the component of this eigenvector on the target subspace */
    PetscScalar projr = 0;
    PetscScalar proji = 0;
    ierr = computeprojection(lambdar, lambdai, xr, xi, &projr, &proji, (void *) targetset);

    ierr = VecScatterBegin(scatter,xr,localx,INSERT_VALUES,
			   SCATTER_FORWARD);
    ierr = VecScatterEnd(scatter,xr,localx,INSERT_VALUES,
			 SCATTER_FORWARD);
    ierr = VecGetArray(localx,&veccomps);

    for (j=0; j<=targetset[0]-1; j++) {
      convergedcomponents[j] += veccomps[targetset[j+1]]
	*veccomps[targetset[j+1]];
    }
    minconvergedcomponent = convergedcomponents[0];
    newminconvergedindex = 0;
    for (j=1; j<=targetset[0]-1; j++) {
      if (minconvergedcomponent>convergedcomponents[j]) {
	minconvergedcomponent = convergedcomponents[j];
	newminconvergedindex = j;
      }
    }
    newtotconverged = totconverged+1;
    ierr = VecRestoreArray(localx, &veccomps);
    /* now add this eigenvector to the deflation space for next 
       iteration */
    ierr = VecDuplicate(xr, Cv+totconverged);
    ierr = VecCopy(xr, Cv[totconverged]);
    totconverged = newtotconverged;
  }

  eigliststruct *lpeiglist = calloc(totconverged, sizeof(eigliststruct));
  for (i=0; i<=totconverged-1; i++) {
    lpeiglist[i].index = i;
    PetscScalar projr = 0;
    PetscScalar proji = 0;
    ierr = computeprojection(0, 0, Cv[i], xi, &projr, &proji, (void *) targetset);
    lpeiglist[i].proj = projr;
  }
  qsort(lpeiglist, totconverged, sizeof(eigliststruct), eiglistcomp);
  /* see how far down I have to go in the list to get the minimum coverage of target space 
     I want */
  for (j=0; j<=targetset[0]-1; j++) {
    convergedcomponents[j] = 0;
  }
  minconvergedcomponent = 0;
  for (i=0; (i<=totconverged-1) && ((1-minconvergedcomponent) >= vectolerance); i++) {
    minconvergedcomponent = 1;
    for (j=0; j<=targetset[0]-1; j++) {
      PetscScalar dotprod = 0;
      ierr = VecZeroEntries(yr);
      ierr = VecSetValue(yr, (PetscInt) targetset[j+1], (PetscScalar) 1,
			 INSERT_VALUES);
      ierr = VecDot(yr, Cv[lpeiglist[i].index], &dotprod);
      convergedcomponents[j] += dotprod*dotprod;	
      if (minconvergedcomponent>convergedcomponents[j]) {
	minconvergedcomponent = convergedcomponents[j];
	minconvergedindex = j;
      }
    }
  }

  needconverged = i;
  printf("Lapack run results:");
  printf("%d eigenstates out of %d determined, only %d actually needed\n",
	 totconverged, nelems, needconverged);

  /* Compare K-S and Lapack eigenvectors */
  for (i=0; i<=kstotconverged-1; i++) {
    PetscScalar dotprod = 0;
    double projnorm = 0;
    double maxprojnorm = 0;
    int maxprojindex = -1;
    for (j=0; j<=totconverged-1; j++) {
      ierr = VecDot(ksCv[kseiglist[i].index], Cv[lpeiglist[j].index], &dotprod);
      if ((projnorm=dotprod*dotprod)>maxprojnorm) {
	maxprojnorm = projnorm;
	maxprojindex = j;
      }
    }
    printf("KS eigenstate %d max match with Lapack eigenstate %d, proj^2=%lf\n",
	   i, maxprojindex, maxprojnorm);
  }
  
  free(kseiglist);
  free(lpeiglist);
  free(targetset);
  free(convergedcomponents);
  ierr = VecDestroy(&Iv);
  ierr = VecDestroy(&yr);
  ierr = ISDestroy(&is);
  ierr = VecScatterDestroy(&scatter); CHKERRQ(ierr);
  for (k=0; k<=kstotconverged-1; k++) {
    ierr = VecDestroy(ksCv+k);
  }
  ierr = PetscFree(ksCv);
  for (k=0; k<=totconverged-1; k++) {
    ierr = VecDestroy(Cv+k);
  }
  ierr = PetscFree(Cv);
  ierr = EPSDestroy(&eps);CHKERRQ(ierr);
  /* I can destroy the matrix H and relatives now */
  ierr = VecDestroy(&xr); CHKERRQ(ierr);
  ierr = VecDestroy(&xi); CHKERRQ(ierr);
  ierr = VecDestroy(&localx); CHKERRQ(ierr);
  ierr = MatDestroy(&H); CHKERRQ(ierr);



  ierr = SlepcFinalize(); CHKERRQ(ierr);
  MPI_Finalize();

  return 0;
}


PetscErrorCode computeprojection(PetscScalar er, PetscScalar ei, 
				 Vec xr, Vec xi, PetscScalar *rr,
				 PetscScalar *ri, void *ctx)
{

  int *targetset = (int *) ctx;
  int targetsize = targetset[0];
  int mytarget;
  PetscScalar val = 0.0, gval = 0.0;
  const PetscScalar *a;
  int j = 0;
  PetscInt low, high;

  VecGetOwnershipRange(xr, &low, &high);
  *rr = 0;
  *ri = 0;

  for (j=1; j<=targetsize; j++) {
    mytarget = targetset[j];
    gval = val = 0;
    if ((mytarget >= low) && (mytarget < high)) {
      VecGetArrayRead(xr,&a);
      val = a[mytarget-low];
      VecRestoreArrayRead(xr, &a);
    }
    MPI_Allreduce(&val, &gval, 1, MPIU_SCALAR, MPI_SUM, PETSC_COMM_WORLD);
    *rr += gval*gval;
  }

  return 0;

}


PetscErrorCode computeprojection1(PetscScalar er, PetscScalar ei, 
				 Vec xr, Vec xi, PetscScalar *rr,
				 PetscScalar *ri, void *ctx)
{

  int *target = (int *) ctx;
  int mytarget = *target;
  PetscScalar val = 0.0, gval = 0.0;
  const PetscScalar *a;
  int j = 0;
  PetscInt low, high;

  VecGetOwnershipRange(xr, &low, &high);
  *rr = 0;
  *ri = 0;

  gval = val = 0;
  if ((mytarget >= low) && (mytarget < high)) {
    VecGetArrayRead(xr,&a);
    val = a[mytarget-low];
    VecRestoreArrayRead(xr, &a);
  }
  MPI_Allreduce(&val, &gval, 1, MPIU_SCALAR, MPI_SUM, PETSC_COMM_WORLD);
  *rr += gval*gval;

  return 0;

}


int findnextnonblankline(FILE *inputstream, size_t *inputsize, 
			 char **inputstring, int *offset) {
    
  ssize_t linelength=0;

  do {
    linelength = getline(inputstring, inputsize, inputstream);
  } while (linelength==0);
    
  if (linelength>0) {
    *offset = 0;
    while (isspace(*(*inputstring+*offset))&&(*offset<linelength)) 
      (*offset)++;
  }

  return linelength;

}


int eiglistcomp(eigliststruct *state1, eigliststruct *state2)
{

  if (state1->proj < state2->proj) return 1;
  else if (state1->proj > state2->proj) return -1;
  else return 0;

}