#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 MATMUL(ccc,aaa,bbb,nra,nca,ncb) \
for(i=0; i<nra; i++)\
{for(j=0; j<ncb; j++)\
{for(k=0; k<nca; k++)\
{ ccc[i][j]+=aaa[i][k]*bbb[k][j];	}}}


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


void read_input();


long i,ijk,j,k,n;

long iii,jjj;

double temp[MAX][MAX];
double P[MAX][MAX];
double K[MAX][MAX];
double V[MAX];

double anum;

double alpha;
double sum;
double ccc;
double IM;

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


int main()
{
    
     printf("\n\n Householder_tridiagonal.cpp   ver 1.2   March 25, 2010 \n"); 
     printf(" By Tom Irvine  Email: tomirvine@aol.com \n");
    
//
//  Reference: Meirovitch p156
//

    read_input();

    PRINTMAT(K,n,n)

    for(ijk=0;ijk<=(n-2);ijk++)
    {
        sum=0.; 
        
        for(j=0;j<n;j++)
        {      
            V[j]=0.;  
        }  
                     
        for(j=ijk+1;j<n;j++)
        {
           sum+=pow(K[ijk][j],2);
           
        }                  
        alpha=sqrt(sum);
        
 //         printf("\n alpha=%8.4g \n",alpha);                
        
        ccc=1+ (K[ijk][ijk+1]/alpha);    

        V[ijk+1]=sqrt(0.5*ccc);

        for(j=ijk+2;j<n;j++)
        {
           V[j]=K[ijk][j]/(2*alpha*V[ijk+1]);

        }

        for(i=0;i<n;i++)
        {
                     
 //            printf("  V[%ld]=%8.4g  \n",i,V[i]);             
                        
           for(j=0;j<n;j++)
           {         
              IM=0.;
              
              if(i==j)
              {
                 IM=1.;     
              }       
                                        
              P[i][j]=IM-2*(V[i]*V[j]);
           } 
        }

//        K=P*K*P
       
        ZERO(temp,n,n)
       
        MATMUL(temp,K,P,n,n,n) 
        
        ZERO(K,n,n)
        
        MATMUL(K,P,temp,n,n,n) 

    }
    
    for(i=0;i<n;i++)
    {
        for(j=0;j<n;j++)
        {
            if(fabs(K[i][j])< 1.0e-9)
            {
                   K[i][j]=0.;              
            } 
            fprintf(pFile[1]," %12.8e \t",K[i][j]);           
        }     
        fprintf(pFile[1],"\n",K[i][j]);                   
    }
   
    printf("\n\n Tridiagonal Matrix \n");
    
    PRINTMAT(K,n,n)

    fclose(pFile[1]);

    printf("\n Press any key to exit.\n");
    getch();
    exit(1);

}
void read_input()
{
     printf("\n The input matrix must be real, symmetric. \n");
     
     printf("\n Enter the input matrix filename. \n");
     
     
     scanf("%s",filename[0]);   
     pFile[0] = fopen(filename[0], "rb");           

     
	 while(pFile[0] == NULL )
	 {

		printf("\n Failed to open file: %s \n", filename[0]);
		printf("\n Please enter the input filename: \n");

		scanf("%s",filename[0]);

		pFile[0] = fopen(filename[0], "rb");

	 }

	 printf(" File: %s opened. \n", filename[0]);     
     
     
     printf("\n\n Enter the output matrix filename \n");
     
     scanf("%s",filename[1]);          
     pFile[1] = fopen(filename[1], "w");


     n=0;
     while(1)
     {
        if(fscanf(pFile[0],"%lf",&anum)<=0)
        {
           break;
        }    
        else
        {
 //          printf(" n=%ld  anum=%8.4g \n",n,anum); 
        }           
        n++;
    }
    rewind(pFile[0]);
   
    n=long(sqrt(n));

    if( n<2 || n>MAX)
	{
		printf("\n\n Maximum number of dofs exceeded. \n");
		printf("\n Press any key to exit.");
		getch();
		exit(1);
	}  
	else
	{
        printf("\n n=%ld \n",n);
    }
  
  
    for(i=0; i< n; i++)
    {
		for(j=0; j<n; j++)
		{
		
            fscanf(pFile[0],"%lf",&anum);
	
	        K[i][j]=anum;
				
	    }
	}
  
    fclose(pFile[0]);

}