/*  pdnmesh - a 2D finite element solver
    Copyright (C) 2001-2005 Sarod Yatawatta <sarod@users.sf.net>  
  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2 of the License, or
  (at your option) any later version.
  
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
  
  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  $Id: eig_lapack.c,v 1.19 2005/02/17 03:43:16 sarod Exp $
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "types.h"

#ifdef USE_MKL 
#include <mkl.h>
#endif /* !USE_MKL */

typedef struct eigenabs_ {
  MY_DOUBLE absval;
  MY_INT loc;
} eigenabs; /* used to sort eigenvalues by magnitude */

static int
my_dspgvx(int ITYPE, char JOBZ, char RANGE, char UPLO, int N,
    double *A,  double *B,  double VL, double VU,
    int IL, int IU, double ABSTOL, int G, double *W, double *Z,
    int LDZ, double *WORK,  int *IWORK, int *IFAIL) {

 int info=0;
#ifdef USE_LAPACK
 extern void dspgvx_( int *ITYPE_, char *JOBZ_, char *RANGE_, char *UPLO_, int *N_,
     double *AP_,  double *BP_,  double *VL_, double *VU_,
     int *IL_, int *IU_, double *ABSTOL_, int *M_, double *W_, double *Z_,
     int *LDZ_, double *WORK_,  int *IWORK_, int *IFAIL_, int *INFO_);
 
 dspgvx_(&ITYPE, &JOBZ, &RANGE, &UPLO, &N, A, B, &VL, &VU,
   &IL, &IU, &ABSTOL, &G, W, Z, &LDZ, WORK, IWORK, IFAIL, &info);
#endif/* !USE_LAPACK */ 
#ifdef USE_MKL
 extern void dspgvx( int *ITYPE_, char *JOBZ_, char *RANGE_, char *UPLO_, int *N_,
     double *AP_,  double *BP_,  double *VL_, double *VU_,
     int *IL_, int *IU_, double *ABSTOL_, int *M_, double *W_, double *Z_,
     int *LDZ_, double *WORK_,  int *IWORK_, int *IFAIL_, int *INFO_);
 
 dspgvx(&ITYPE, &JOBZ, &RANGE, &UPLO, &N, A, B, &VL, &VU,
   &IL, &IU, &ABSTOL, &G, W, Z, &LDZ, WORK, IWORK, IFAIL, &info);
#endif /* !USE_MKL */
 return(info);
}

static double
my_dlamch(char CMACH)
{
#ifdef USE_LAPACK
  extern  double  dlamch_ (char *CMACHp);
  return  dlamch_ (&CMACH);
#endif /* !USE_LAPACK */
#ifdef USE_MKL
  extern  double  dlamch (char *CMACHp);
  return  dlamch(&CMACH);
#endif /* !USE_MKL */
  return(0);
}


void
solve_helmholtz_lapack(MY_DOUBLE **P,MY_DOUBLE **T, MY_DOUBLE **x,  MY_DOUBLE *ev, MY_INT  N, MY_INT G)
{
 /* solves the symmetrix, generalized eigenvalue problem 
    * (P-lambda. T)x= 0
    * P, T = N by N matrices, only lower triangle
    * x = eigenvectors to be found size (N x M)
    * using LAPACK 
    */

 int i,j;
 double *Af,  *Bf, *WORK, *W, *Z;
 int *IWORK,*IFAIL;
 int info;
 MY_DOUBLE ab=0; /* normalization */

 update_status_bar("Start solving Eigenproblem.");
 if((Af= (double*) calloc(N*(N+1)/2, sizeof(double)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 if((Bf= (double*) calloc(N*(N+1)/2, sizeof(double)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 /* copy A and B to Af and Bf */
 for (j=1;j<=N;j++) {
  for (i=j;i<=N;i++) {
    Af[i+(j-1)*(2*N-j)/2-1]=P[i-1][j-1];
    Bf[i+(j-1)*(2*N-j)/2-1]=T[i-1][j-1];
  }
 }

 if((W= (double*) calloc(N, sizeof(double)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 /* find first M eigenvalue, eigenvector pairs */ 
 if((Z= (double*) calloc(G*N, sizeof(double)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 /* workspace */
 if((WORK= (double*) calloc(8*N, sizeof(double)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 if((IWORK= (int*) calloc(5*N, sizeof(int)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 /* indices of eigenvectors faild to converge */
 if((IFAIL= (int*) calloc(N, sizeof(int)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

 info=my_dspgvx(1,'V','I','L',N,Af,Bf,0,0,1,G,2*my_dlamch('S'),G,W,Z,N,WORK,IWORK,IFAIL);
 if ( info ) {
   fprintf(stderr,"solve_helmholtz_lapack: something is wrong. flag=%d\n",info);
   if (info >= N) {
    fprintf(stderr,"solve_helmholtz_lapack: we have a singular matrix. bailing out\n");
    exit(1);
   }
 }

#ifdef DEBUG
 printf("solve_helmholtz_lapack: found %d eigenvalue/eigenvector pairs\n",G);
#endif
 for (i=0;i<G;i++) {
#ifdef DEBUG
  printf("%lf :[",W[i]);
#endif
   ev[i]=W[i]; /* copy eigenvalues */
  for(j=0;j<N;j++) {
#ifdef DEBUG
   printf("%lf ",Z[N*i+j]);
#endif
   x[j][i]=Z[N*i+j];
   ab+=x[j][i]*x[j][i];
  }
  /* normalization loop */
  ab=sqrt(ab);
  for(j=0;j<N;j++) {
   x[j][i]/=ab;
  }
#ifdef DEBUG
  printf("]\n");
#endif
 }

 free(Af);
 free(Bf);
 free(W);
 free(Z);
 free(WORK);
 free(IWORK);
 free(IFAIL);

 update_status_bar("Done solving Eigenproblem.");
}

void
solve_helmholtz_lapack_sparse(SPMat *P,SPMat *T, MY_DOUBLE **x,  MY_DOUBLE *ev, MY_INT  N, MY_INT G)
{
 /* solves the symmetrix, generalized eigenvalue problem 
    * (P-lambda. T)x= 0
    * P, T = N by N matrices, only lower triangle
    * x = eigenvectors to be found size (N x M)
    * using LAPACK 
    */

 int i,j;
 double *Af,  *Bf, *WORK, *W, *Z;
 int *IWORK,*IFAIL;
 int info;
 MY_DOUBLE ab=0; /* normalization */

 update_status_bar("Start solving Eigenproblem.");
 if((Af= (double*) calloc(N*(N+1)/2, sizeof(double)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 if((Bf= (double*) calloc(N*(N+1)/2, sizeof(double)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 /* copy A and B to Af and Bf */
 for (j=1;j<=N;j++) {
  for (i=j;i<=N;i++) {
    Af[i+(j-1)*(2*N-j)/2-1]=SPM_e(P,i-1,j-1);
    Bf[i+(j-1)*(2*N-j)/2-1]=SPM_e(T,i-1,j-1);
  }
 }

 if((W= (double*) calloc(N, sizeof(double)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 /* find first M eigenvalue, eigenvector pairs */ 
 if((Z= (double*) calloc(G*N, sizeof(double)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 /* workspace */
 if((WORK= (double*) calloc(8*N, sizeof(double)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 if((IWORK= (int*) calloc(5*N, sizeof(int)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 /* indices of eigenvectors faild to converge */
 if((IFAIL= (int*) calloc(N, sizeof(int)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

 info=my_dspgvx(1,'V','I','L',N,Af,Bf,0,0,1,G,2*my_dlamch('S'),G,W,Z,N,WORK,IWORK,IFAIL);
 if ( info ) {
   fprintf(stderr,"solve_helmholtz_lapack: something is wrong. flag=%d\n",info);
   if (info >= N) {
    fprintf(stderr,"solve_helmholtz_lapack: we have a singular matrix. bailing out\n");
    exit(1);
   }
 }

#ifdef DEBUG
 printf("solve_helmholtz_lapack: found %d eigenvalue/eigenvector pairs\n",G);
#endif
 for (i=0;i<G;i++) {
#ifdef DEBUG
  printf("%lf :[",W[i]);
#endif
   ev[i]=W[i]; /* copy eigenvalues */
  for(j=0;j<N;j++) {
#ifdef DEBUG
   printf("%lf ",Z[N*i+j]);
#endif
   x[j][i]=Z[N*i+j];
   ab+=x[j][i]*x[j][i];
  }
  /* normalization loop */
  ab=sqrt(ab);
  for(j=0;j<N;j++) {
   x[j][i]/=ab;
  }
#ifdef DEBUG
  printf("]\n");
#endif
 }

 free(Af);
 free(Bf);
 free(W);
 free(Z);
 free(WORK);
 free(IWORK);
 free(IFAIL);

 update_status_bar("Done solving Eigenproblem.");
}



/* wrapper for LAPACK non symmetric generalized eigenvalue routine */
static int
my_dggev(char JOBVL, char JOBVR, int N, double *A, int LDA, double *B, int LDB,
   double *ALPHAR, double *ALPHAI, double *BETA, double *VL, int LDVL,
   double *VR, int LDVR, double *WORK, int LWORK) {

 int info=0;
#ifdef USE_LAPACK
   extern void
   dggev_(char *JOBVL_, char *JOBVR_, int *N_, double *A_, int *LDA_,
   double *B_, int *LDB_, double *ALPHAR_, double *ALPHAI_, double *BETA_,
   double *VL_, int *LDVL_, double *VR_, int *LDVR_, double *WORK_, int *LWORK_, int *INFO);

 dggev_(&JOBVL,&JOBVR,&N,A,&LDA,B,&LDB,ALPHAR,ALPHAI,BETA,VL,&LDVL,VR,&LDVR,WORK, &LWORK, &info);
#endif /* !USE_LAPACK */
#ifdef USE_MKL
   extern void
   dggev(char *JOBVL_, char *JOBVR_, int *N_, double *A_, int *LDA_,
   double *B_, int *LDB_, double *ALPHAR_, double *ALPHAI_, double *BETA_,
   double *VL_, int *LDVL_, double *VR_, int *LDVR_, double *WORK_, int *LWORK_, int *INFO);

 dggev(&JOBVL,&JOBVR,&N,A,&LDA,B,&LDB,ALPHAR,ALPHAI,BETA,VL,&LDVL,VR,&LDVR,WORK, &LWORK, &info);
#endif /* !USE_MKL */

 return(info);

}

/* function used to compare abs values of eigenvalues for sorting */
static int
compare_abs_eigenvalues(const void *a, const void *b) {
 const eigenabs *da=(const eigenabs *)a;
 const eigenabs *db=(const eigenabs *)b;
  
 return((ABS(da->absval)>=ABS(db->absval)?1:-1));
}

/* function to solve generelized eigenvalue problem
   when the matrices are singular 
  Px = lambda Tx
  P, T = N by N matrices -full 
  x = N by N matrix of generelized eigenvectors 
  ev = generalized eigenvalues
*/
void
solve_generalized_eig_lapack(SPMat P,SPMat T, MY_DOUBLE **x, MY_DOUBLE *ev,  MY_INT  N, MY_INT G) {

 double *Af, *Bf, *al_R, *al_I, *bet;
 double *Vf,  *WORK;
 eigenabs *a_abs;
 int i,j,info,work_size;
 double dummy;

 if((Af= (double*) malloc(N*N*sizeof(double)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 if((Bf= (double*) malloc(N*N*sizeof(double)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 /* copy A and B to Af and Bf */
 for (i=0;i<N;i++) {
  for (j=0;j<N;j++) {
    Af[j+i*N]=SPM_e(&P,i,j);
    Bf[j+i*N]=SPM_e(&T,i,j);
  }
 }

 if((al_R= (double*) malloc(N*sizeof(double)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 if((al_I= (double*) malloc(N*sizeof(double)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 if((bet= (double*) malloc(N*sizeof(double)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 if((a_abs= (eigenabs*) malloc(N*sizeof(eigenabs)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 if((Vf= (double*) malloc(N*N*sizeof(double)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

 /* arrays for eigenvalues lambda=(al_R+j*al_I)/bet */
 /* compute only the right generelized eigenvectors */
  /* Vf=stores right eigenvectors NxN */ 
 if((WORK= (double*) malloc((size_t)(10*N)*sizeof(double)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 } 
 info=my_dggev('N','V',N,Af,N,Bf,N,al_R,al_I,bet,&dummy,1,Vf,N,WORK,-1);
#ifdef DEBUG
 if (!info  && (int)WORK[0]> 10*N) {
  printf("solve_generalized_eig_lapack: WORK size=%lf %d\n",WORK[0],(int)WORK[0]);
 } 
#endif

 work_size=(int)WORK[0];  
 if((WORK= (double*) realloc((void*)WORK,(size_t)(work_size)*sizeof(double)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 info=my_dggev('N','V',N,Af,N,Bf,N,al_R,al_I,bet,&dummy,1,Vf,N,WORK,10*N);
 if ( !info ) {
#ifdef DEBUG
   printf("optimal work area=%lf used %d\n",WORK[0],work_size);
   printf("solve_generalized_eig_lapack: generelized eigenvalues\n");
#endif
   /* calculate magnitudes of eigenvalues to sort them */
   for (i=0;i <N;i++) {

#ifdef DEBUG
    printf("%d/%d=(%lf+%lf)/%lf\n",i,N,al_R[i],al_I[i],bet[i]);
#endif
    if ( ABS(bet[i]) > TOL ) {
     a_abs[i].absval=(al_R[i]*al_R[i]+al_I[i]*al_I[i])/(bet[i]*bet[i]);
     /* eliminate zero eigenvalues too */
     if ( a_abs[i].absval <TOL ) {
      a_abs[i].absval=10000000.0;
     }
    } else {
     a_abs[i].absval=10000000.0;
    }
    a_abs[i].loc=i;
   }
   /* sort eigenvalues */
   qsort(a_abs,N,sizeof(eigenabs),compare_abs_eigenvalues);
#ifdef DEBUG
   /* print sorted values */
   for (i=0; i<N; i++) {
    printf("sort:%d=%lf\n",a_abs[i].loc,pow(a_abs[i].absval,0.25));
   }
#endif
 } 

 /* copy the G number of eigenvectors to v, according to the 
   order of smallest eigenvalues i.e. N-G,N-G+1,...,N */
  /* Vf[N^2-GN]..Vf[N^2-(G-1)N-1],....,Vf[N^2-N]...,Vf[N^2-1] */
  /* the locations are given be first G elements of a_abs[] array */
  for (j=0;j<G;j++) {
#ifdef DEBUG
   printf("column %d->%d\n",j,a_abs[j].loc);
#endif 
   /* copy eigenvalue assuming bet[] not zero*/
   ev[j]=al_R[a_abs[j].loc]/bet[a_abs[j].loc];
    for (i=0;i<N;i++) {
       x[i][j]=Vf[(a_abs[j].loc)*N+i];
#ifdef DEBUG
       printf("%d  %lf\n",i,x[i][j]);
#endif
   }
  } 
  free(Af);
  free(Bf);
  free(al_R);
  free(al_I);
  free(bet);
  free(a_abs);
  free(WORK); 
  free(Vf);

}