text: /*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !!!!!!! FEAST Driver banded example !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !!!!!!! solving Ax=eBx with A real and B spd --- A and B banded matrix!!!
    !!!!!!! by Eric Polizzi- 2009-2012!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/
  
  #include <stdio.h> 
  #include <stdlib.h> 
  #include <sys/time.h>
  #include <mpi.h>
  
  #include "feast.h"
  #include "feast_banded.h"
  int main(int argc, char **argv) { 
    /*!!!!!!!!!!!!!!!!! Feast declaration variable */
    int  feastparam[64]; 
    float epsout;
    int loop,LDA,LDB;
    char UPLO='F'; // ! 'L' or 'U' also fine
  
    /*!!!!!!!!!!!!!!!!! Matrix declaration variable */
    FILE *fp;
    char name[]="../../system1";
    int  N,nnz,kl,ku;
    float *A,*B;
  
    /*!!!!!!!!!!!!!!!!! Others */
    struct timeval t1, t2;
    int  i,j,k,n2,err;
    int  M0,M,info;
    float Emin,Emax,trace,dum;
    float *X; //! eigenvectors
    float *E,*res; //! eigenvalue+residual
  
  /*********** MPI *****************************/
  int rank,numprocs;
  MPI_Init(&argc,&argv); 
  MPI_Comm_size(MPI_COMM_WORLD,&numprocs); 
  MPI_Comm_rank(MPI_COMM_WORLD,&rank); 
  /*********************************************/
  
    /*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !!!!!!!!!!!!!!! read input file in coordinate format!!!!!!!
      !!!!!!!!!!!!!!!!transform  in banded format directly !!!!!!
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/
  
    // !! find kl,ku (kl==ku here)
    kl=0;
    fp = fopen (name, "r");
    err=fscanf (fp, "%d%d\n",&N,&nnz);
    for (k=0;k<=nnz-1;k++){
      err=fscanf(fp,"%d%d%f%f\n",&i,&j,&dum,&dum);
      if (abs(i-j)>kl){ kl=abs(i-j); }
    }
    fclose(fp);
    ku=kl;
    // !! form the banded matrices A and B
    LDA=(kl+ku+1);
    LDB=LDA;
    n2=N*(kl+ku+1);
    A=calloc(n2,sizeof(float));
    B=calloc(n2,sizeof(float));
    for (i=0;i<=n2-1;i++){
      *(A+i)=(float) 0.0;
      *(B+i)=(float) 0.0;
    };
  
    fp = fopen (name, "r");
    err=fscanf (fp, "%d%d\n",&N,&nnz); 
    for (k=0;k<=nnz-1;k++){
      err=fscanf(fp,"%d%d",&i,&j);
      err=fscanf(fp,"%f%f\n",A+(j-1)*LDA+ku+(i-j),B+(j-1)*LDB+ku+(i-j));
    };
    fclose(fp);
  
      
    /*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !!!!!!!!!!!!!!!!!!!!!! INFORMATION ABOUT MATRIX !!!!!!!!!!!!!!!
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/
  if (rank==0) {
    printf("dense matrix -system1- size %.d\n",N);
    printf("bandwidth %d \n",kl+ku+1);
  }
    gettimeofday(&t1,NULL);
    /*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !!!!!!!!!!!!!!!!!!! FEAST in banded format !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/
  
    /*!!! search interval [Emin,Emax] including M eigenpairs*/
    Emin=(float) 0.18;
    Emax=(float) 1.0;
    M0=25; // !! M0>=M
  
    /*!!!!!!!!!!!!! ALLOCATE VARIABLE */
    E=calloc(M0,sizeof(float));  // eigenvalues
    res=calloc(M0,sizeof(float));// eigenvectors 
    X=calloc(N*M0,sizeof(float));// residual 
  
  
    /*!!!!!!!!!!!!  FEAST */
    feastinit(feastparam);
    feastparam[0]=1;  /*change from default value */
    sfeast_sbgv(&UPLO,&N,&kl,A,&LDA,&kl,B,&LDB,feastparam,&epsout,&loop,&Emin,&Emax,&M0,E,X,&M,res,&info);
  
    gettimeofday(&t2,NULL);
    /*!!!!!!!!!! REPORT !!!!!!!!!*/
  if (rank==0) {
    printf("FEAST OUTPUT INFO %d\n",info);
    if (info==0) {
      printf("*************************************************\n");
      printf("************** REPORT ***************************\n");
      printf("*************************************************\n");
      printf("# of processors %d \n",numprocs);
      printf("SIMULATION TIME %f\n",(t2.tv_sec-t1.tv_sec)*1.0+(t2.tv_usec-t1.tv_usec)*0.000001);
      printf("# Search interval [Emin,Emax] %f %f\n",Emin,Emax);
      printf("# mode found/subspace %d %d \n",M,M0);
      printf("# iterations %d \n",loop);
      trace=(float) 0.0;
      for (i=0;i<=M-1;i=i+1){
        trace=trace+*(E+i);
      }	  
      printf("TRACE %f\n", trace);
      printf("Relative error on the Trace %f\n",epsout );
      printf("Eigenvalues/Residuals\n");
      for (i=0;i<=M-1;i=i+1){
        printf("   %d %f %f\n",i,*(E+i),*(res+i));
      }
    }
  }
  MPI_Finalize(); /************ MPI ***************/
    return 0;
  }
  
  
  
  
  
  
  