/******************************************************************

  zero_cfrac.c generates zero points and associated residues
  of a continued fraction expansion terminated at 2n+2 level
  of the Ferm-Dirac function, which is derived from Gauss's
  hypergeometric function.

  This code is distributed under the constitution of GNU-GPL.

 (C) Taisuke Ozaki (AIST-RICS)
  
  Log of zero_cfrac.c:

     14/July/2005  Released by T.Ozaki

******************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <search.h>
#include <string.h>

#include "openmx_common.h"
#include "lapack_prototypes.h"

#ifdef nompi
#include "mimic_mpi.h"
#else
#include "mpi.h"
#endif


static void Eigen_DGGEVX( int n, double **a, double **s, double *eval,
                          double *resr, double *resi );


void zero_cfrac( int n, dcomplex *zp, dcomplex *Rp ) 
{
  static int i,j,N;
  static double **a,**s,*eval,*resr,*resi;

  /* check input parameters */

  if (n<=0){
    printf("\ncould not find the number of zeros\n\n");
    MPI_Finalize();
    exit(0);
  }

  /* the total number of zeros including minus value */

  N = 2*n + 1;

  /* allocation of arrays */

  a = (double**)malloc(sizeof(double*)*(N+2));
  for (i=0; i<(N+2); i++){
    a[i] = (double*)malloc(sizeof(double)*(N+2));
  }

  s = (double**)malloc(sizeof(double*)*(N+2));
  for (i=0; i<(N+2); i++){
    s[i] = (double*)malloc(sizeof(double)*(N+2));
  }

  eval = (double*)malloc(sizeof(double)*(n+3));
  resr = (double*)malloc(sizeof(double)*(n+3));
  resi = (double*)malloc(sizeof(double)*(n+3));

  /* initialize arrays */

  for (i=0; i<(N+2); i++){
    for (j=0; j<(N+2); j++){
      a[i][j] = 0.0;
      s[i][j] = 0.0;
    }
  }

  /* set matrix elements */

  s[2][1] =  1.0;
  s[2][2] = -0.5;

  for (i=3; i<=N; i++){
     s[i][i-1] =  -0.5;
     s[i-1][i] =   0.5;
  }

  a[1][1] = -2.0;
  a[1][2] =  1.0;
  a[2][2] = -1.0;

  for (i=3; i<=N; i++){
    a[i][i] = -(2.0*(double)i - 3.0);
  }

  /* diagonalization */

  Eigen_DGGEVX( N, a, s, eval, resr, resi );

  for (i=0; i<n; i++){
    zp[i].r = 0.0;
    zp[i].i = eval[i+1];
    Rp[i].r = resr[i+1];
    Rp[i].i = 0.0;
  }

  /* print result */
  /*
  for (i=1; i<=n; i++){
    printf("i=%4d  eval=%18.14f   resr=%18.15f resi=%18.15f\n",i,eval[i],resr[i],resi[i]);
  }
  */

  /* free of arrays */

  for (i=0; i<(N+2); i++){
    free(a[i]);
  }
  free(a);

  for (i=0; i<(N+2); i++){
    free(s[i]);
  }
  free(s);

  free(eval);
  free(resr);
  free(resi);
}













void Eigen_DGGEVX( int n, double **a, double **s, double *eval, double *resr, double *resi )
{
  static int i,j,k,l,num;

  static char balanc = 'N';
  static char jobvl = 'V';
  static char jobvr = 'V';
  static char sense = 'B';
  static double *A;
  static double *b;
  static double *alphar;
  static double *alphai;
  static double *beta;
  static double *vl;
  static double *vr;
  static double *lscale;
  static double *rscale;
  static double abnrm;
  static double bbnrm;
  static double *rconde;
  static double *rcondv;
  static double *work;
  static double *tmpvecr,*tmpveci;
  static INTEGER *iwork;
  static INTEGER lda,ldb,ldvl,ldvr,ilo,ihi;
  static INTEGER lwork,info;
  static logical *bwork; 
  static double sumr,sumi,tmpr,tmpi;
  static double *sortnum;

  lda = n;
  ldb = n;
  ldvl = n;
  ldvr = n;

  A = (double*)malloc(sizeof(double)*n*n);
  b = (double*)malloc(sizeof(double)*n*n);
  alphar = (double*)malloc(sizeof(double)*n);
  alphai = (double*)malloc(sizeof(double)*n);
  beta = (double*)malloc(sizeof(double)*n);

  vl = (double*)malloc(sizeof(double)*n*n);
  vr = (double*)malloc(sizeof(double)*n*n);

  lscale = (double*)malloc(sizeof(double)*n);
  rscale = (double*)malloc(sizeof(double)*n);

  rconde = (double*)malloc(sizeof(double)*n);
  rcondv = (double*)malloc(sizeof(double)*n);

  lwork = 2*n*n + 12*n + 16;
  work = (double*)malloc(sizeof(double)*lwork);

  iwork = (INTEGER*)malloc(sizeof(INTEGER)*(n+6));
  bwork = (logical*)malloc(sizeof(logical)*n);

  tmpvecr = (double*)malloc(sizeof(double)*(n+2));
  tmpveci = (double*)malloc(sizeof(double)*(n+2));

  sortnum = (double*)malloc(sizeof(double)*(n+2));

  /* convert two dimensional arrays to one-dimensional arrays */

  for (i=0; i<n; i++) {
    for (j=0; j<n; j++) {
       A[j*n+i]= a[i+1][j+1];
       b[j*n+i]= s[i+1][j+1];
    }
  }

  /* call dggevx_() */

  F77_NAME(dggevx,DGGEVX)(
           &balanc, &jobvl, & jobvr, &sense, &n, A, &lda, b, &ldb,
           alphar, alphai, beta, vl, &ldvl, vr, &ldvr, &ilo, &ihi,
           lscale, rscale, &abnrm, &bbnrm, rconde, rcondv, work,
           &lwork, iwork, bwork, &info );

  if (info!=0){
    printf("Errors in dggevx_() info=%2d\n",info);
  }

  /*
  for (i=0; i<n; i++){
    printf("i=%4d  %18.13f %18.13f %18.13f\n",i,alphar[i],alphai[i],beta[i]);
  }
  printf("\n");
  */

  num = 0;
  for (i=0; i<n; i++){

    if ( 1.0e-13<fabs(beta[i]) && 0.0<alphai[i]/beta[i] ){

      /* normalize the eigenvector */

      for (j=0; j<n; j++) {

        sumr = 0.0;
        sumi = 0.0;

        for (k=0; k<n; k++) {
          sumr += s[j+1][k+1]*vr[i*n    +k];
          sumi += s[j+1][k+1]*vr[(i+1)*n+k];
        }
        
        tmpvecr[j] = sumr;
        tmpveci[j] = sumi;
      }

      sumr = 0.0;
      sumi = 0.0;

      for (k=0; k<n; k++) {
        sumr += vl[i*n+k]*tmpvecr[k] + vl[(i+1)*n+k]*tmpveci[k];
        sumi += vl[i*n+k]*tmpveci[k] - vl[(i+1)*n+k]*tmpvecr[k];
      }

      /* calculate zero point and residue */

      eval[num+1] = alphai[i]/beta[i];
      tmpr =  vr[i*n]*vl[i*n] + vr[(i+1)*n]*vl[(i+1)*n];
      tmpi = -vr[i*n]*vl[(i+1)*n] + vr[(i+1)*n]*vl[i*n];
      resr[num+1] =  tmpi/sumi;
      resi[num+1] = -tmpr/sumi;

      num++;
    }
    else{
      /*
      printf("i=%4d  %18.13f %18.13f %18.13f\n",i+1,alphar[i],alphai[i],beta[i]);
      */
    }
  }

  /* check round-off error */

  for (i=1; i<=num; i++){
    if (1.0e-8<fabs(resi[i])){
      printf("Could not calculate zero points and residues due to round-off error\n");
      MPI_Finalize();
      exit(0);
    }
  }

  /* sorting */

  qsort_double(num,eval,resr);

  /* free arraies */

  free(A);
  free(b);
  free(alphar);
  free(alphai);
  free(beta);

  free(vl);
  free(vr);

  free(lscale);
  free(rscale);

  free(rconde);
  free(rcondv);

  free(work);

  free(iwork);
  free(bwork);

  free(tmpvecr);
  free(tmpveci);
  free(sortnum);
}