//#include "lapack.h"
#include<math.h>
#include <stdio.h>
#include <stdlib.h>


void calcula(int **m,int k,double *mf,double *mfo,double *mfe,double *S,double *c,double *cg,double *conng,int flag, FILE *fco, FILE *fcoconn, FILE *fmf, char **nodos)
{
  int i,j,kk,aristas,sum;
  double * A, *V2, *fb, *C, *EDL, *vp, *EL;


  //  char *JOBZ, *UPLO;
  long N, M, L, LDA, LDB, LDC, LWORK, INFO, incx, incy, one;
  double *W, *WORK,*mfv,*mfov,*mfev;
  int *grado;
  double alpha,beta,media();
  void nomemo();
  

  N=k;
  M=k;
  L=k;
  LDA=k;
  LDB=k;
  LDC=k;
  LWORK=3*k-1;
  incx=1;
  incy=1;
  alpha=1;
  beta=0;
  one=1;

  mfv= (double *) malloc (k * sizeof(double)); if (mfv==NULL) nomemo();
  mfov= (double *) malloc (k * sizeof(double)); if (mfov==NULL) nomemo();
  mfev= (double *) malloc (k * sizeof(double)); if (mfev==NULL) nomemo();
  A= (double *) malloc (k*k * sizeof(double)); if (A==NULL) nomemo();
  V2= (double *) malloc (k*k * sizeof(double)); if (V2==NULL) nomemo();
  C= (double *) malloc (k*k * sizeof(double)); if (C==NULL) nomemo();
  EL= (double *) malloc (k*k * sizeof(double));if (EL==NULL) nomemo(); 
  W= (double *) malloc (k * sizeof(double)); if (W==NULL) nomemo();
  //  fb= (double *) malloc (k * sizeof(double)); if (fb==NULL) nomemo();
  EDL= (double *) malloc (k * sizeof(double));if (EDL==NULL) nomemo(); 
  vp= (double *) malloc (k * sizeof(double)); if (vp==NULL) nomemo();
  WORK= (double *) malloc ((3*k-1) * sizeof(double)); if (WORK==NULL) nomemo();
  grado = (int *) malloc(k * sizeof(int)); if (grado==NULL) nomemo();
 
  //  printf("Inicio calcula\n");

  for(i=0;i<k;i++)
    {
      sum=0;
      for(j=0;j<k;j++)
	{
	  A[i*k+j]=(double) m[i][j];
	  sum+=m[i][j];
	}
      grado[i]=sum;
    }

  /***** AUTOVALORES Y AUTOVECTORES *****/

  N=k;
  LDA=k;
  LWORK=3*k-1;

  dsyev ("V","U",&N,A,&LDA,W,WORK,&LWORK,&INFO,sizeof(char),sizeof(char)); 

  //    printf("\nINFO=%d \nautovectores\n",INFO);

  for (i=0;i<k*k;i++)
    {
      V2[i]=A[i]*A[i];
    }

  for (i=0;i<k;i++)
    {
      EDL[i]=exp(W[i]);
    }

  dgemv("N",&M,&N,&alpha,V2,&LDA,EDL,&incx,&beta,mfv,&incy,sizeof(char));

  //  mf=fb;
  //  *mf=(double) media(fb,k);


  for (i=0;i<k;i++)
    {
      EDL[i]=sinh(W[i]);
    }

  dgemv("N",&M,&N,&alpha,V2,&LDA,EDL,&incx,&beta,mfov,&incy,sizeof(char));

  //  mfo=fb;
  //  *mfo=(double) media(fb,k);


  for (i=0;i<k;i++)
    {
      EDL[i]=cosh(W[i]);
    }

  dgemv("N",&M,&N,&alpha,V2,&LDA,EDL,&incx,&beta,mfev,&incy,sizeof(char));

  //  mfe=fb;
  //  *mfe=(double) media(fb,k);


  /***** Calculos de Excel *****/

  for (i=0;i<k;i++)
    {
      vp[i]=fabs(A[(k-1)*k+i]);               /* autovector principal */
      EDL[i]=log(vp[i]);
    }

  //  printf("flag=%d,%d",flag,k);
  if (flag==1)
    {
      for (i=0;i<k;i++)
	{
	  fprintf(fmf,"%s\t%g\t%g\t%g\t%g\n",nodos[i],mfv[i],mfov[i],mfev[i],vp[i]);
	  //	  printf("%s\t%g\t%g\t%g\t%g\n",nodos[i],mf[i],mfo[i],mfe[i],vp[i]);
	}
    }

  *mf=(double) media(mfv,k);
  *mfo=(double) media(mfov,k);
  *mfe=(double) media(mfev,k);
  //  printf("\tmedia de fb = %g\n\n",*mf);
  //  printf("\tmedia centralidades odd = %g\n\n",*mfo);
  //  printf("\tmedia centralidades even = %g\n\n",*mfe);

  alpha=-1/log(2);
  dgemv("T",&N,&one,&alpha,vp,&LDA,EDL,&incx,&beta,S,&incy,sizeof(char));
  //  printf("\tS = %g\n\n",*S);
       
  /***** Communicability *****/

  for(i=0;i<k*k;i++)
    {
      EL[i]=0;
    }

  for (i=0;i<k;i++)
    {
      EL[i*k+i]=exp(W[i]);
    }

 alpha=1;
 dgemm("N","N",&M,&N,&L,&alpha,A,&LDA,EL,&LDB,&beta,C,&LDC,sizeof(char),sizeof(char));
 dgemm("N","T",&M,&N,&L,&alpha,C,&LDA,A,&LDB,&beta,V2,&LDC,sizeof(char),sizeof(char));

  *cg=0;
  for(i=0;i<k-1;i++)
    {
      for(j=i+1;j<k;j++)
	{
	  *cg+=fabs(V2[i*k+j]);  // poner en valor absoluto
	}
    }
  *cg=*cg*2/(k*k-k);
  //  printf("\tcg = %g,%d\n\n",*cg,k);

  if(flag==1)
    {
      for(i=0;i<k-1;i++)
	{
	  for(j=i+1;j<k;j++)
	    {
	      fprintf(fco,"%s\t%s\t%d\t%d\t%g\n",nodos[i],nodos[j],grado[i],grado[j],fabs(V2[i*k+j]));
	    }
	}
    }

  *conng=0;
  aristas=0;
  for(i=0;i<k-1;i++)
    {
      for(j=i+1;j<k;j++)
	{
	  if (m[i][j]==1)
	    {
	      aristas++;
	      *conng+=fabs(V2[i*k+j]);
	    }
	}
    }
  *conng=*conng/aristas;
  //  printf("\tconng = %g\n\n",*conng);

  if(flag==1)
    {
      for(i=0;i<k-1;i++)
	{
	  for(j=i+1;j<k;j++)
	    {
	      if (m[i][j]==1)
		{
		  fprintf(fcoconn,"%s\t%s\t%d\t%d\t%g\n",nodos[i],nodos[j],grado[i],grado[j],fabs(V2[i*k+j]));
		}
	    }
	}
    }

  /*  
  for(i=0;i<k;i++)
    {
      for(j=0;j<k;j++)
	{
	  C[i*k+j]=0;
	  for(kk=0;kk<k;kk++)
	    {
	      C[i*k+j]+=(A[kk*k+i]*EL[kk*k+kk])*A[kk*k+j];
	    }
	  C[i*k+j]=fabs(C[i*k+j]);
	}
    }

  *cg=0;
  for(i=0;i<k-1;i++)
    {
      for(j=i+1;j<k;j++)
	{
	  *cg+=C[i*k+j];  // poner en valor absoluto
	}
    }
  *cg=*cg*2/(k*k-k);
    printf("\tcg = %g\n\n",*cg);

  *conng=0;
  aristas=0;
  for(i=0;i<k-1;i++)
    {
      for(j=i+1;j<k;j++)
	{
	  if (m[i][j]==1)
	    {
	      aristas++;
	      *conng+=C[i*k+j];
	    }
	}
    }
  *conng=*conng/aristas;
    printf("\tconng = %g\n\n",*conng);
  */

  free(mfv);
  free(mfov);
  free(mfev);
  free(A);
  free(V2);
  free(C);
  free(EL);
  free(W);
  //  free(fb);
  free(EDL);
  free(vp);
  free(WORK);
  
}


double media(double * vector, int k)
{
  int i;
  double sum=0;
  for(i=0;i<k;i++)
    {
      sum+=vector[i];
    }
  return(sum/k);
}
