#include <iostream>
#include <stdlib.h>
#include <math.h>
#define MaxD 10

using namespace std;

float Matrix[MaxD][MaxD];
float Q[MaxD][MaxD];
float tempQ[MaxD][MaxD];
float R[MaxD][MaxD];
float H[MaxD][MaxD];   //QR分解時的H矩陣 
int N;

int sgn(float Number)
{
    int retVal = 1;
    if(Number<0) 
        retVal = -1;
    return retVal;    
}

float MultHM(int row, int col)
{
    float retVal = 0;
    for(int i=0; i<N; i++)
    {
        retVal += H[row][i]*Matrix[i][col];
    }
    return retVal;
}

float MultQH(int row, int col)
{
    float retVal = 0;
    for(int i=0; i<N; i++)
    {
        retVal += Q[row][i]*H[i][col];
    }
    return retVal;
}

void MultRQtoMatix()
{
    int i=0,j=0,k=0;
    float temp = 0.0;
    for(i=0; i<N; i++)
    {
        for(j=0; j<N; j++)
        {
            temp = 0.0;
            for(k=0; k<N; k++)
            {
                temp += Matrix[i][k]*Q[k][j];
            }
            R[i][j] = temp;
        }
    }
    
    for(i=0; i<N; i++)
    {
        for(j=0; j<N; j++)
        {
            Matrix[i][j] = R[i][j];
        }
    }
}

void QRAnalyze()
{
    int i=0,j=0,k=0;

    float V[MaxD];         //QR分解時的v向量 
    float a;               //a = sgn(a11)|a1| 
    float b;               //b = a(a+a11) 
    
    for(i=0; i<N-1; i++)   //用i來控制列 
    {
        a = 0;
        for(j=0; j<N; j++)
        {
            V[j] = Matrix[j][i];
        }
        
        for(j=i; j<N; j++)
        {
            a += V[j]*V[j];
        }
        
        a = sqrt(a);
        a *= sgn(V[i]);
        b = a*(a + V[i]);
        
        V[i] = V[i] + a;
        
        for(j=0; j<N; j++)
        {
            for(k=0; k<N; k++)
            {
                H[j][k] = 0.0;                
            }
        }
        
        for(j=i; j<N; j++)
        {
            for(k=i; k<N; k++)
            {
                H[j][k] = V[j]*V[k]*(-1/b);
            }
        }
        
        for(j=0; j<N; j++)
        {
            H[j][j] += 1;
        }

        //Ai = Hi * Ai-1
        for(j=0; j<N; j++)
        {
            for(k=0; k<N; k++)
            {
                R[j][k] = MultHM(j,k);
                tempQ[j][k] = MultQH(j,k);
            }
        }
        for(j=0; j<N; j++)
        {
            for(k=0; k<N; k++)
            {
                Matrix[j][k] = R[j][k];
                Q[j][k] = tempQ[j][k];
            }
        }
    }
}

int main(int argc, char *argv[])
{
    int i=0,j=0;
    do
    {
        printf("請輸入矩陣維數 不大于10\n");
        scanf("%d",&N);
    }while(N > 10);
    
    printf("請輸入矩陣元素:\n");
    for(i=0; i<N; i++)
    {
        for(j=0; j<N; j++)
        {
            printf("Matrix[%d][%d]=",i,j);
            scanf("%f",&Matrix[i][j]);
        }
    }
    
    /*  測試數據 
    Matrix[0][0]=1.0;
    Matrix[0][1]=2.0;
    Matrix[0][2]=3.0;
    Matrix[1][0]=2.0;
    Matrix[1][1]=3.0;
    Matrix[1][2]=4.0;
    Matrix[2][0]=2.0;
    Matrix[2][1]=1.0;
    Matrix[2][2]=4.0;
    */
    for(int num=0; num<10; num++)
    {
        for(i=0; i<N; i++)
        {
            for(j=0; j<N; j++)
            {
                Q[i][j] = 0;
                R[i][j] = 0;
            }
            Q[i][i] = 1;
        }
        QRAnalyze();
        //Matrix = R*Q;
        MultRQtoMatix();
    }
    printf("矩陣的特征值為：\n");
    for(i=0; i<N; i++)
    {
        printf("a[%d] = %5f\n",i,Matrix[i][i]);
    }

    system("PAUSE");	
    return 0;
}