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


#define MAX 5000

#define PRINTMAT(aaa,nnn,mmm) \
for(iii=0; iii<nnn; iii++)\
{ for(jjj=0; jjj<mmm; jjj++)\
{printf(" %11.4g \t",aaa[iii][jjj]);}printf("\n");}

#define IDENTITY(aaa,nnn) \
for(i=0; i<nnn; i++)\
{for(j=0; j<nnn; j++)\
{aaa[i][j]=0.;}aaa[i][i]=1.;}

#define SWITCHI(aaa,bbb); \
itemp=aaa;\
aaa=bbb;\
bbb=itemp;\

#define SWITCH(aaa,bbb); \
temp=aaa;\
aaa=bbb;\
bbb=temp;\


#define TRANSPOSE(aaa,bbb) \
for(i=0; i<n; i++)\
{ for(j=0; j<n; j++)\
{aaa[i][j]=bbb[j][i];}}

#define ZERO(aaa) \
for(i=0; i<n; i++)\
{ for(j=0; j<n; j++)\
{aaa[i][j]=0.;}}

#define ZEROG(aaa,mmm,nnn) \
for(i=0; i<mmm; i++)\
{ for(j=0; j<nnn; j++)\
{aaa[i][j]=0.;}}

#define ZEROV(aaa) \
for(i=0; i<n; i++)\
{aaa[i]=0.;}


#define MATMUL(ccc,aaa,bbb) \
for(i=0; i<n; i++)\
{for(j=0; j<n; j++)\
{for(k=0; k<n; k++)\
{ ccc[i][j]+=aaa[i][k]*bbb[k][j];	}}}

#define MATMULU(ccc,aaa,bbb,nnn,mmm) \
for(i=0; i<nnn; i++)\
{for(j=0; j<mmm; j++)\
{for(k=0; k<mmm; k++)\
{ ccc[i][j]+=aaa[i][k]*bbb[k][j];	}}}

#define MATMUL_TRANS(ccc,aaa,bbb) \
for(i=0; i<n; i++)\
{for(j=0; j<n; j++)\
{for(k=0; k<n; k++)\
{ ccc[i][j]+=aaa[k][i]*bbb[k][j];	}}}


#define MATEQUAL(aaa,bbb) \
for(iii=0; iii<n; iii++)\
{for(jjj=0; jjj<n; jjj++)\
{aaa[iii][jjj]=bbb[iii][jjj];}}\

#define TRANSPOSE(aaa,bbb) \
for(i=0; i<n; i++)\
{ for(j=0; j<n; j++)\
{aaa[i][j]=bbb[j][i];}}


void files(void);
void inverse(void);

void LDLT(void);
void LLT(void);
void Ktilda(void);
void eigen(void);
void PTKP(void);
void CalculateP(void);
void zerops(void);
void interactive(void);
void output(void);

const double pi =atan2(0.,-1.);

double tpi=2*pi;

float aabb,anum,ah,aq,s;
float temp,theta;

float a[MAX][MAX],aa[MAX][MAX],ac[MAX][MAX],b[MAX][MAX],bb[MAX][MAX];
float c[MAX][MAX],cc[MAX][MAX],ctm[MAX][MAX];

// float anm[MAX][MAX];
// float tempr[MAX][MAX];

// float cct[MAX][MAX];

float d[MAX][MAX];

long ngw[MAX];

float dt;


float fn[MAX];
float ktil[MAX][MAX];
float l[MAX][MAX],lt[MAX][MAX],li[MAX][MAX],lit[MAX][MAX];

long nt;
float sr;

float omegan[MAX];
float p[MAX][MAX],pt[MAX][MAX];
float stiff[MAX][MAX],x[MAX][MAX];

float prs[MAX][MAX];

float time_end;
float time_start;

long SWP;

long im[MAX];
long i,ic,iflag,ii,iii,ik,imne;
long j,jj,jjj,k,kk,kki,n,num;
long iv,jv;
long klast;

long mmm,nnn;

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

  float sum;
   
  long ijk;   

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

FILE *pFile[20];
char filename[20][FILENAME_MAX];


int main()
{

	printf("\n jacobi.cpp, ver 1.7, February 5, 2010");
	printf("\n\n by Tom Irvine ");
	printf("\n Email:  tomirvine@aol.com ");

    printf("\n\n This program uses the Jacobi method to solve the eigenvalue");
    printf("\n problem: \n\n      K - lambda M = 0");
	printf("\n\n where K and M are matrices.");

	printf("\n\n Restrictions: \n  1. Matrix M must be positive definite. \n  2. Matrices K and M must be symmetric. \n");

	     
//   Decomposition of a positive definite mass matrix into LDLT form and into LLT form.
//   Jacobi transformation method for eigenvalue and eigen vector analysis.

	    files();

		interactive();

  
	    printf("\n zerops ");
		zerops();
	
	    printf("\n LDLT ");
		LDLT();

	    printf("\n LLT ");
		LLT();
        
	    printf("\n Ktilda ");
		Ktilda();
		
	    printf("\n eigen ");
		eigen();
		
	    printf("\n output ");
	    output();
// 
  	    fclose(pFile[2]);
		fclose(pFile[3]);
		fclose(pFile[14]);
		fclose(pFile[15]); 
		
		int itm;
		
		printf("\n Multiply eigenvectors by transformation matrix? \n");
		printf("  1=yes  2=no \n");
		scanf("%d",&itm);
		
//		printf("\n ref 1 \n");
		
		if(itm==1)
		{
//  		    printf("\n ref 2 \n");
                          
            while(1)
            {
               printf("\n\n Enter the transformation matrix filename \n");
               scanf("%s",filename[2]);   
               pFile[2] = fopen(filename[2], "rb");
               
               if(pFile[2]>0){break;}                  
            }                 
 //     		printf("\n ref 3 \n");            
                  
            long nc;
            
            printf("\n");
            
            nc=0;
            for(i=0;i<MAX;i++)
            {
                if(fscanf(pFile[2],"%f",&anum)<=0)
                {
                   break;
                }               
                nc++;
            }
            rewind(pFile[2]);
   
            nc=long(nc/n);  
            
//            printf("\n  nc=%ld  n=%ld \n",nc,n);     
                  
            for(i=0; i< nc; i++)
            {
			   for(j=0; j<n; j++)
			   {
				   fscanf(pFile[2],"%f",&anum);
				   ctm[i][j]=anum;
			   }
		    }         
            fclose(pFile[2]);
            
            ZERO(a)
            
            MATMULU(a,ctm,cc,nc,n)  
            
            int iord;

		    printf("\n Reorder eigenvectors? \n");
		    printf("  1=yes  2=no \n");
		    scanf("%d",&iord);
            
            if(iord==1)
            {
               while(1)
               {
                  printf("\n\n Enter the reorder vector filename \n");
                  scanf("%s",filename[2]);   
                  pFile[2] = fopen(filename[2], "rb");
               
                  if(pFile[2]>0){break;}                  
               }      
               printf("\n");                   
              
               for(i=0; i< nc; i++)
               {
				   fscanf(pFile[2],"%ld",&iii);
				   ngw[i]=iii;
		       }        
               fclose(pFile[2]);         
               
               ZERO(b)
               
               for(i=0;i<nc;i++)
               {
                   for(j=0;j<n;j++)
                   {
                       iii=ngw[i];            
                       b[iii][j]=a[i][j];
                   }    
               }
               MATEQUAL(b,a)

            }  
            
            printf("\n\n Transformed eigenvectors \n\n");                      
            PRINTMAT(a,nc,n)
            
            while(1)
            {
                 printf("\n Enter the transformed eigenvector matrix filename \n");
                 scanf("%s",filename[4]);   
                 pFile[4] = fopen(filename[4], "w");     
                 
                 if(pFile[4]>0){break;}             
            }  
            
            for(i=0; i<nc; i++)
            {   
                for(j=0; j<n; j++)
                {              
                    fprintf(pFile[4],"%12.5e \t",a[i][j]);
                }
               fprintf(pFile[4],"\n");	        
            }
            fclose(pFile[4]);
                  
        }
                    
		printf("\n\n Calculation complete.  Press any key to exit.");
		getch();
		exit(1);

}

void files(void)
{   
    
   while(pFile[4] == NULL )
   {
       printf("\n Enter the mass matrix filename \n");
       scanf("%s",filename[4]);   
       pFile[4] = fopen(filename[4], "rb");                  
   }  
   
   while(pFile[5] == NULL )
   {
       printf("\n Enter the stiffness matrix filename \n");
       scanf("%s",filename[5]);   
       pFile[5] = fopen(filename[5], "rb");                  
   }  


// Bathe, p 578

	strcpy(filename[1],"LDLT.out");
	strcpy(filename[2],"LLT.out");

//   printf("\n\n The input file must..... ");

//   printf("\n\n Enter input filename\n");
//   scanf("%s",&filename[0]);
//   pFile[0]=fopen(filename[0], "rb");

//   printf("\n Enter output filename for LDLT decomposition.\n");
//   scanf("%s",&filename[1]);
   pFile[1]=fopen(filename[1], "w");

//   printf("\n Enter output filename for LLT decomposition.\n");
//   scanf("%s",&filename[2]);
   pFile[2]=fopen(filename[2], "w");

   printf("\n Enter output filename for full text output. \n");
   scanf("%s",&filename[3]);
   pFile[3]=fopen(filename[3], "w"); 
 
   
   printf("\n Enter output filename mass normalized eigenvectors. \n");
   scanf("%s",&filename[14]);
   pFile[14]=fopen(filename[14], "w");
   
   printf("\n Enter output filename for natural frequencies (Hz). \n");
   scanf("%s",&filename[15]);
   pFile[15]=fopen(filename[15], "w");       
      
/*   
   for(i=0; i<=2; i++)
   {
      if(pFile[i] == NULL )
      {
         printf(" Failed to open file: %s \n", filename[i]);
         fclose(pFile[i]);
      }
      else
      {
         printf(" File: %s opened. \n", filename[i]);
      }
   }
   printf("\n");
*/
}

void inverse(void)
{
	IDENTITY(x,MAX)

	for(i=0; i<n; i++)  // x column
	{
	
		MATEQUAL(aa,l)

        for(j=0; j< n; j++)
		{
			x[j][i]/=aa[j][j];
			aa[j][j]=1.;

			for(k=j+1; k<n; k++)
			{
			    x[k][i]+=-aa[k][j]*x[j][i];
			}
		}
	}

	MATEQUAL(li,x)

}

void LDLT(void)
{
//		printf("\n LDLT ref 1");
		for(i=0; i<n; i++)
		{
			for(j=0; j<n; j++)
			{
				if( fabs(a[i][j]) > 1.0e-80)
				{
					im[i]=j;
					break;
				}
			}
		}

//		printf("\n LDLT ref 2");
		
		ZERO(d)
		ZERO(l)

		d[0][0]=a[0][0];


//  i = row number 
      
//		printf("\n LDLT ref 3");
		for(i=1; i<n; i++)
		{

//  j = column number

			l[i][0]=a[i][0]/d[0][0];
        
			if(i>=1)
			{
				for( j=im[i]; j<= (i-1); j++)
				{
					s=0.;

					for(k=0; k<j; k++)
					{
						s+= d[k][k]*l[i][k]*l[j][k];
					}
					l[i][j]=(a[i][j]-s)/(d[j][j]);
				}
			}

//  find d coefficients
         
			s=0.;

			for(kk=0; kk <= (i-1); kk++)
			{
				s+=pow( l[i][kk],2. )*d[kk][kk];
			}
			d[i][i]=a[i][i]-s;

			if(fabs(d[i][i]) < 1.0e-25)
			{
				printf("\n Zero on diagonal. \n"); 
				break;
			}
			l[i][i]=1.;

		}
		l[0][0]=1.;

//		printf("\n LDLT ref 4");
		fprintf(pFile[1],"\n\n\n L matrix \n\n");
		for(i=0; i<n; i++)
		{  
				for(j=0; j<n; j++)
				{              
					    fprintf(pFile[1],"%12.5e \t",l[i][j]);
				}
                fprintf(pFile[1],"\n");		
		}

//		printf("\n LDLT ref 5");
		fprintf(pFile[1],"\n\n\n D matrix \n\n");
		for(i=0; i<n; i++)
		{
				for(j=0; j<n; j++)
				{              
					    fprintf(pFile[1],"%12.5e \t",d[i][j]);
				}
                fprintf(pFile[1],"\n");			
		}
		fprintf(pFile[1],"\n\n\n LT matrix \n\n");
		
//		printf("\n LDLT ref 6");
		for(i=0; i<n; i++)
		{  
				for(j=0; j<n; j++)
				{              
					    fprintf(pFile[1],"%12.5e \t",l[j][i]);
				}
                fprintf(pFile[1],"\n");			
		}
//		printf("\n LDLT ref 7");

		
}
void LLT(void)
{
	
		fprintf(pFile[2],"\n\n\n L matrix \n\n");
				
		for(i=0; i<n; i++)
		{
			if(d[i][i] < 0. )
			{
				printf("\n LLT error. \n Press any key to exit.");
		        getch();
	
				exit(1);
			}

			d[i][i] =sqrt(d[i][i]);
		}


		ZERO(b)

		MATMUL(b,l,d)

		MATEQUAL(l,b)
		
		for(i=0; i<n; i++)
		{  
				for(j=0; j<n; j++)
				{              
				 	 
			         fprintf(pFile[2],"%12.5e \t",l[i][j]);
					 
					 aa[i][j]=l[i][j];
				}
                fprintf(pFile[2],"\n");			
		}
		fprintf(pFile[2],"\n\n\n LT matrix \n\n");


	    ZERO(lt)

        TRANSPOSE(lt,l)

		for(i=0; i<n; i++)
		{  
				for(j=0; j<n; j++)
				{              
					    fprintf(pFile[2],"%12.5e \t",lt[i][j]);
				}
                fprintf(pFile[2],"\n");			
		}

		

		inverse();
		

		fprintf(pFile[2],"\n\n\n L inverse matrix \n\n");
		for(i=0; i<n; i++)
		{  
				for(j=0; j<n; j++)
				{              
					    fprintf(pFile[2],"%12.5e \t",li[i][j]);
				}	
                fprintf(pFile[2],"\n");				
		}

		ZERO(lit)

		for(i=0; i<n; i++)
		{  
				for(j=0; j<n; j++)
				{     
					lit[i][j]=li[j][i];
				}
		}
		
		fprintf(pFile[2],"\n\n\n L inverse T matrix \n\n");
		for(i=0; i<n; i++)
		{  
				for(j=0; j<n; j++)
				{              
					    fprintf(pFile[2],"%12.5e \t",lit[i][j]);
				}	
                fprintf(pFile[2],"\n");						
		}

		ZERO(b)
		MATMUL(b,l,lt)

		fprintf(pFile[2],"\n\n\n L LT matrix \n\n");
		for(i=0; i<n; i++)
		{  
				for(j=0; j<n; j++)
				{              
					    fprintf(pFile[2],"%12.5e \t",b[i][j]);	
				}			
                fprintf(pFile[2],"\n",b[i][j]);   				
		}

		
		ZERO(b)
		MATMUL(b,l,li)

		fprintf(pFile[2],"\n\n\n L LI matrix \n\n");
		for(i=0; i<n; i++)
		{  
				for(j=0; j<n; j++)
				{              
					    fprintf(pFile[2],"%12.5e \t",b[i][j]);
				}			
                fprintf(pFile[2],"\n");     				
		}

}

void Ktilda(void)
{

	  ZERO(b)
	  ZERO(ktil)
	
	  MATMUL(b,stiff,lit)

	  MATMUL(ktil,li,b)


      fprintf(pFile[2],"\n\n\n K tilda \n\n");

	  for(i=0; i<n; i++)
	  {  
				for(j=0; j<n; j++)
				{              
					    fprintf(pFile[2],"%12.5e \t",ktil[i][j]);
				}
                fprintf(pFile[2],"\n");		
		}
}
void eigen(void)
{
//    printf("\n ref 1 "); 
    float ee;
    
    float max_error;

    float cbefore[MAX];
    
    int iflag=0;
    int jflag=0;

	ZERO(ac)

	MATEQUAL(c,ktil)
	
//	printf("\n Enter the number of sweeps.  Suggest >= 8 \n");
//    scanf("%ld",&SWP);

    SWP=200;
  //  printf("\n ref 2 "); 
	for(imne=1; imne <= SWP; imne++)   // sweeps
	{
		fprintf(pFile[2],"\n\n SWP = %ld \n",imne);
//		 printf("\n SWP = %ld \n",imne);

		kki=0;
	
		for( iv=0; iv<n; iv++)
		{
			for( jv=iv+1; jv<n; jv++)
			{
				iflag=0;
				CalculateP();

				if(iflag==9){break;}

				PTKP();

				kki++;
			}
			if(iflag==9){break;}
		}

		fprintf(pFile[2],"\n\n\n Transformed K tilda Matrix\n\n");

		for(i=0; i<n; i++)
		{  
				for(j=0; j<n; j++)
				{              
                   fprintf(pFile[2],"%12.5e \t",c[i][j]);
				}
                fprintf(pFile[2],"\n");		
		}


		fprintf(pFile[2],"\n\n\n Eigenvalues \n\n");
	
		for(i=0; i<n; i++)
		{
			fprintf(pFile[2],"\n %12.5e ",c[i][i]);
		}
		
		max_error=0.;
		
		if(imne>2)
		{
           if( cbefore[i]>1.0e-10 && c[i][i]>1.0e-10 )
           {
		        for(i=0; i<n; i++)
		        {               
                    ee=log10(cbefore[i]/c[i][i]);
                    ee=fabs(ee);
                
                    if( ee>max_error)
                    {
                       max_error=ee;
                    }
                }
           }
           if(max_error<0.1)  //  0.1 db
           {
               jflag=1;
           }   
        }
		
		
		for(i=0; i<n; i++)
		{
			cbefore[i]=c[i][i];
		}


		fprintf(pFile[2],"\n\n\n Eigenvector Matrix \n\n");

		for(i=0; i<n; i++)
		{  
				for(j=0; j<n; j++)
				{              
				 	 if( j <(n-1) )
					 {
					    fprintf(pFile[2],"%12.5e \t",prs[i][j]);
					 }
					 else
					 {
			            fprintf(pFile[2],"%12.5e \n",prs[i][j]);
					 }
				}		
		}


		fprintf(pFile[2],"\n\n\n Transformed Eigenvector Matrix \n\n");

        ZERO(cc)

		MATMUL(cc,lit,prs)

		float an[MAX];
		
		for(i=0; i<n; i++)
		{
				an[i]=0.;
		}

		for(i=0; i<n; i++)
		{ 
				for(j=0; j<n; j++)
				{              
				 	 if( j <(n-1) )
					 {
					    fprintf(pFile[2],"%12.5e \t",cc[i][j]);
					 }
					 else
					 {
			            fprintf(pFile[2],"%12.5e \n",cc[i][j]);
					 }
					 an[i]+=pow( cc[j][i],2. );
				}		
		}

//   printf("\n ref 3 "); 
		fprintf(pFile[2],"\n\n\n Normalized Transformed Eigenvector Matrix \n\n");
	
		for(i=0; i<n; i++)
		{  
				for(j=0; j<n; j++)
				{              
                   fprintf(pFile[2],"%12.5e \t",cc[i][j]/sqrt(an[j]));
				}
                fprintf(pFile[2],"\n");		
		}
		
		if(jflag==1)
		{
           printf("\n %ld sweeps \n",imne);         
                    
           break;
        }

	}
//    printf("\n ref 4 "); 

		if(jflag==0)
		{
           printf("\n Warning: maximum error limit exceeded.  %ld sweeps \n",imne);         
                    
        }
}

void CalculateP(void)
{
//	float cf,con1,con2,theta,cf_thresh,sf_thresh;


//	cf_thresh = pow( 10.,-2.*(2.*imne) );
//	sf_thresh = pow( 10.,-24);
 
 // coupling factor to determine whether sweep is necessary
//	cf= sqrt( pow(c[iv][jv],2.) / ( c[iv][iv]*c[jv][jv] ));


//	iflag=9;
//	for(ic=0; ic<n; ic++)
//	{
//		con1=fabs( c[ic][ic] - ac[ic][ic] )/c[ic][ic];

//		if(con1 > sf_thresh) // convergence not yet achieved
//		{
//			iflag=0;
//		}
//	}


//    if( cf > cf_thresh && iflag == 0)

// perform cyclic Jacobi procedure rather than threshold Jacobi method

	iflag=0;

	if(iflag == 0)
	{

		if( fabs( c[iv][iv]-c[jv][jv] ) > 1.0e-20)
		{
			theta = 0.5*atan( (2.*c[iv][jv]) / ( c[iv][iv]-c[jv][jv] )  ); 
		}
		else
		{
			theta = pi/4.;
		}

//	fprintf(pFile[2],"\n\n\n theta= %12.4g \n",theta);

		ZERO(p)

		for( i=0; i<n; i++)
		{
			p[i][i]=1.;		
		}

		p[iv][iv]=cos(theta);
		p[iv][jv]=-sin(theta);
		p[jv][iv]=sin(theta);
		p[jv][jv]=cos(theta);

		for( i=0; i<n; i++)
		{
			for( j=0; j<n; j++)
			{
				pt[i][j]=p[j][i];
			}
		}


		ZERO(cc)

		MATMUL(cc,prs,p)

		MATEQUAL(prs,cc)

	}
}

void PTKP(void)
{

	ZERO(b)

	MATEQUAL(ac,c)

	MATMUL(b,c,p)

	ZERO(c)

	MATMUL(c,pt,b);
	
}
void zerops(void)
{
        ZERO(prs)

		for(i=0; i<MAX; i++)
		{
			prs[i][i]=1.;
		}

}
void interactive(void)
{
     
    int idm;
     
    printf("\n divide mass by 386? \n  1=yes 2=no \n");
    scanf("%d",&idm);

// reading

   n=0;
   for(i=0;i<MAX;i++)
   {
      if(fscanf(pFile[4],"%f",&anum)<=0)
      {
         break;
      }               
      n++;
   }
   rewind(pFile[4]);
   
   n=long(sqrt(n));

    if( n<2 || n>MAX)
	{
		printf("\n\n Illegal n value. \n");
		printf("\n Press any key to exit.");
		getch();
		exit(1);
	}  
	else
	{
        printf("\n n=%ld \n",n);
    }

//	printf("\n\n Enter Mass Matrix M \n");

        printf("\n reading mass \n");

		for(i=0; i< n; i++)
		{
			for(j=0; j<n; j++)
			{
				fscanf(pFile[4],"%f",&anum);
				if(idm==1)
                {
                    anum/=386.;
                }
				a[i][j]=anum;
			}
		}

//  Let stiff[i][j] = stiffness matrix
		
//	printf("\n\n Enter Stiffness Matrix K \n");

        printf("\n reading stiffness \n");
        
		for(i=0; i< n; i++)
		{
			for(j=0; j<n; j++)
			{
				fscanf(pFile[5],"%f",&anum);
				stiff[i][j]=anum;
			}
		}
	
        fclose(pFile[4]);
        fclose(pFile[5]);  
        
}
void output(void)
{   
 //  printf("\n ref 1 \n");
  
	float temp;

	long ik;

        for(i=0;i<n-1;i++)
		{
              for(j=i+1; j<n; j++)
			  {
			      if(c[j][j] < c[i][i] )
				  {
				      SWITCH(c[j][j],c[i][i]);

					  for(ik=0; ik<n; ik++)
					  {
						  SWITCH(cc[ik][j],cc[ik][i]);
					  }
				  }
			  }
		  }
    //    printf("\n ref 2 \n");
   
     
		fprintf(pFile[3],"\n\n\n Eigenvalues \n\n");
	
		for(i=0; i<n; i++)
		{
			fprintf(pFile[3],"\n lambda %2d = %10.5g \t sqrt(lambda %d )= %10.5g",i+1,c[i][i],i, sqrt( c[i][i]) );
			
			omegan[i]=sqrt( c[i][i]);
			fn[i]=omegan[i]/tpi;
		}
		fprintf(pFile[3],"\n");

	    fprintf(pFile[3],"\n\n\n Natural Frequencies (Hz) \n\n");
	     printf("\n\n\n Natural Frequencies (Hz) \n\n");  
	
		for(i=0; i<n; i++)
		{
		     fprintf(pFile[3]," %ld  %10.5g \n",i+1,fn[i]);    
             fprintf(pFile[15]," %10.5g \n",fn[i]);    		     
		      printf(" %ld  %10.5g \n",i+1,fn[i]);                     
        }
	
		float an[MAX];
		
		for(i=0; i<n; i++)
		{
				an[i]=0.;
		}
 
		for(i=0; i<n; i++)
		{ 
				for(j=0; j<n; j++)
				{              
					 an[i]+=pow( cc[j][i],2. );
				}		
		}



		fprintf(pFile[3],"\n\n\n Mass Normalized Eigenvector Matrix (Column format) \n\n");

        ZERO(p)

        MATMUL(p,a,cc)		
	
//        ZERO(tempr)

//        MATMUL(tempr,a,cc)
        
//       TRANSPOSE(cct,cc)
	
//	    ZERO(anm)

        ZERO(b)
	
//	    MATMUL(anm,cc,tempr)
	    
//   	    MATMUL_TRANS(anm,cc,tempr)    

        MATMUL_TRANS(b,cc,p)   
 		
 	    for(i=0; i< n; i++)
		{
			fprintf(pFile[3],"%10d \t",i+1);
		}
		fprintf(pFile[3],"\n\n");
			
	    float aq[MAX];
		
		for(i=0; i<n; i++)
		{  
             aq[i]=sqrt(b[i][i]);           	
        }
		
        for(i=0; i<n; i++)
		{  
		    for(j=0; j<n; j++)
		    { 
                 cc[i][j]/=aq[j];         
            }
        } 


		for(i=0; i<n; i++)
		{  
				for(j=0; j<n; j++)
				{              
				 	 if( j <(n-1) )
					 {
					    fprintf(pFile[3],"%10.4e \t",cc[i][j]);
                        fprintf(pFile[14],"%12.5e \t",cc[i][j]);		    
					 }
					 else
					 {
			            fprintf(pFile[3],"%10.4e \n",cc[i][j]);
		                fprintf(pFile[14],"%12.5e \n",cc[i][j]);		            
					 }
				}		
		}
}