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

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

#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 ADD_MATRICES(aaa,bbb,ccc,nnn) \
for(iii=0; iii<nnn; iii++)\
{for(jjj=0; jjj<nnn; jjj++)\
{aaa[iii][jjj]=bbb[iii][jjj]+ccc[iii][jjj];}}\

#define MATEQUAL_TIMES_CONSTANT(aaa,bbb,ccc,nnn) \
for(iii=0; iii<nnn; iii++)\
{for(jjj=0; jjj<nnn; jjj++)\
{aaa[iii][jjj]=ccc*bbb[iii][jjj];}}\


#define VECTMULUT(ccc,aaa,nnn,mmm) \
for(i=0; i<nnn; i++)\
{for(j=0; j<mmm; j++)\
{ ccc[i][j]=aaa[i]*aaa[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 IDENTITY(aaa,nnn) \
for(i=0; i<nnn; i++)\
{ for(j=0; j<nnn; j++)\
{aaa[i][j]=0.;}aaa[i][i]=1.;}

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

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

#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];}}

#define MAX 2000
#define MAX2 MAX*MAX

void data_entry(void);
void QR(void);
void output(void);

double nvk2;

long i,j,k;
long n,nv;
long m;

long iii,jjj;

long iv;

double anum;
double norm;

double A[MAX][MAX];
double Q[MAX][MAX];
double Qk[MAX][MAX];
double ident[MAX][MAX];
double Hk[MAX][MAX];
double vvt[MAX][MAX];
double temp[MAX][MAX];
double v[MAX];
double ak[MAX];
double ek[MAX];
double vk[MAX];

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

int main()
{

    printf("\n\n QR_Householder.m   ver 1.1   March 16, 2010 \n");
    printf(" By Tom Irvine  Email: tomirvine@aol.com \n");
    
    printf("\n The input matrix must be square. \n");
    n=MAX;
    
    ZERO(A,n,n);

    data_entry();

    QR();
  
    output();  
  
    printf("\n Calculation complete. \n Press any key to exit. \n\n");
    getch();
   
}
void output(void)
{
    for(i=0;i<n;i++)
    {
       for(j=0;j<n;j++)
       {
           if(fabs(A[i][j])<1.0e-10)
           {
               A[i][j]=0.;                     
           }
           if(fabs(Q[i][j])<1.0e-10)
           {
               Q[i][j]=0.;                     
           }                       
       }                  
    }
    
    
    for(i=0; i<n; i++)
    { 
        for(j=0; j<n; j++)
        {
             fprintf(pFile[1]," %12.6g \t",Q[i][j]);
        }
        fprintf(pFile[1],"\n");
    }    
    
    for(i=0; i<n; i++)
    { 
        for(j=0; j<n; j++)
        {
             fprintf(pFile[2]," %12.6g \t",A[i][j]);  // R=A
        }
        fprintf(pFile[2],"\n");
    }  
   
    
    printf("\n\n Q \n");
    PRINTMAT(Q,n,n)    
    
    printf("\n\n R \n");
    PRINTMAT(A,n,n)
    
    fclose(pFile[1]);
    fclose(pFile[2]); 
    
    printf("\n\n Output files \n");
    
    printf(" %s \n",filename[1]);    
    printf(" %s \n",filename[2]);  
        
}
void QR(void)
{
     double sign_ak;
     double norm_ak;
     
     m=n;  
     
     long kv;   
    
     IDENTITY(Q,m)
     
     IDENTITY(ident,m)
     
     long LAK=n;
  
     for(kv=0;kv<(n-1);kv++)
     {  
        
//        printf(" kv=%ld  \n",kv);
        
//        PRINTMAT(A,n,n)   
        
        for(i=0;i<LAK;i++)
        {
           ak[i]=A[kv+i][kv];                       
        }       
        
//        printf("\n ak \n");
//        PRINTVECT(ak,LAK)
        
        norm_ak=0;          

        for(i=0;i<LAK;i++)
        {
           norm_ak+=pow(ak[i],2);
        }
        
        norm_ak=sqrt(norm_ak);
        
 //       printf(" norm_ak = %8.4g \n",norm_ak);
        
        sign_ak=ak[0]/fabs(ak[0]);
      
            
        ZEROV(ek,LAK)
      
        ek[0]=1;

        nvk2=0;
        
        for(i=0;i<LAK;i++)
        {
            vk[i]=ak[i] + sign_ak*norm_ak*ek[i];
            
            nvk2+=pow(vk[i],2);
        }
        
 //       printf("\n vk \n");
 //       PRINTVECT(vk,LAK);  
        
        
        MATEQUAL(Qk,ident)  
        
       
    
//  construct reflector
    
//       Hk=ident - 2*vk*vk'/(norm(vk)^2);    

        for(i=0;i<LAK;i++)
        {    
            for(j=0;j<LAK;j++)
            {
                 Hk[i][j]=vk[i]*vk[j];                              
            }
        }


        for(i=0;i<LAK;i++)
        {    
            for(j=0;j<LAK;j++)
            {
                 Hk[i][j]=ident[i][j]-2*Hk[i][j]/nvk2;        
            }
        } 
//                printf("\n\n Hk \n");
//                   PRINTMAT(Hk,LAK,LAK) 
//        Qk=[eye(k-1)  zeros(k-1,m-k+1);  zeros(m-k+1,k-1)  Hk];

        for(i=0;i<LAK;i++)
        {    
            for(j=0;j<LAK;j++)
            {
                 Qk[kv+i][kv+j]=Hk[i][j];        
            }
        }
//              printf("\n\n Qk \n");  
//           PRINTMAT(Qk,n,n) 
//        A=Qk*A;

        ZERO(temp,n,n)
        MATMULU(temp,Qk,A,n,n)   
        MATEQUAL(A,temp)   
        
            
//         printf("\n\n A \n");
//                PRINTMAT(A,n,n)   
        
//       Q=Q*Qk;
        
        ZERO(temp,n,n)
        MATMULU(temp,Q,Qk,n,n)   
        MATEQUAL(Q,temp) 
        
 //        printf("\n\n Q \n");
 //               PRINTMAT(Q,n,n)       
        
        LAK--; 
    } 
        
}
void data_entry(void)
{
    while(pFile[0] == NULL )
    {
       printf("\n\n Enter the input matrix filename \n");
       scanf("%s",filename[0]);   
       pFile[0] = fopen(filename[0], "rb");                  
    }     

    n=0;
    while(1)
    {
      if(fscanf(pFile[0],"%lf",&anum)<=0)
      {
         break;
      }               
      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);
    }

    printf("\n reading input marix \n");

    for(i=0; i< n; i++)
    {
			for(j=0; j<n; j++)
			{
				fscanf(pFile[0],"%lf",&anum);
				
				A[i][j]=anum;
			}
    }
    fclose(pFile[0]);
    
    strcpy(filename[1],"Q.out");
    strcpy(filename[2],"R.out");  
    
    pFile[1] = fopen(filename[1], "w");  
    pFile[2] = fopen(filename[2], "w");         
}