# hw1

###### tags: `tutorials`

#include <iostream>
#include <iomanip>

using namespace std;

void printMatrix(float arr[][4]);
void qr_decomposition(float arr[4][4], int temp);
float global[4][4];

int main()
{
    // input matrix
    float arr[4][4] = { {1,2,3,4},
                        {0,-1,3,2},
                        {0,0,3,3},
                        {0,0,0,2}};
    qr_decomposition(arr, 1);
}

// print the 4x4 matrix
void printMatrix(float arr[][4]) 
{
    for (int i = 0; i < 4; i++) {
        for (int j = 0; j < 4; j++)cout << fixed <<setprecision(5)  << arr[i][j] << " ";
        cout << endl;
    }cout << endl;
}

// assign one matrix to another
void copyMatrix(float a[4][4], float b[4][4]) {
    for (int i = 0; i < 4; i++) {
        for (int j = 0; j < 4; j++) {
            a[i][j] = b[i][j];
        }
    }
}

// multiply two matrix and assign the value result to global variable
void multiply(float a[4][4], float b[4][4]) {
    for (int i = 0; i < 4; i++) {
        float temp = 0;
        for (int j = 0; j < 4; j++) {
            global[i][j] = a[i][0] * b[0][j] 
                         + a[i][1] * b[1][j] 
                         + a[i][2] * b[2][j] 
                         + a[i][3] * b[3][j];
        }
    }
}

// caculate the value of two matrix
float multiply_1d(float a[4], float b[4]) {
    float temp = 0;
    for (int i = 0; i < 4; i++) {
        temp += (a[i] * b[i]);
    }
    return temp;
}

/* 
apply QR algorithm

@param arr : input matrix
@param temp : iterator time
*/
void qr_decomposition(float matrix[4][4], int temp) {

    if (temp > 15) {
        cout << "最後結果" << endl;
        printMatrix(matrix);
        return;
    }

    float q[4][4];
    float q_transpose[4][4];

    // initialize matrix q and matrix q_transpose
    copyMatrix(q, matrix);
    copyMatrix(q_transpose, matrix);

    // caculate Q1~Q4 and assign to matrix q

    for (int i = 0; i < 4; i++) {
        for (int a = 0; a < 4; a++) {
            for (int b = 0; b < 4; b++) {
                q_transpose[a][b] = q[b][a];
            }
        }
        
        /*
           iterate Q1~Q4
           the Qk needs k-1 times iteration
        */
        int iterator = i;
        while (iterator) {
            
            // do the prog algorithm
            for (int j = 0; j < 4; j++) {
                q[j][i] -= (multiply_1d(q_transpose[i], q_transpose[iterator -1]) / multiply_1d(q_transpose[iterator - 1], q_transpose[iterator - 1]))
                            * q[j][iterator - 1];
            }
            iterator--;
        } 
    }

    /*
    *
    * 
    // 將每個 q 矩陣列進行單位向量化 (e.g. A/|A|)
    for (int i = 0; i < 4; i++) {
        float base = 0;
        for (int j = 0; j < 4; j++) {
            base += pow(q[j][i],2);
        }

        for (int j = 0; j < 4; j++) {
            q[j][i] /= sqrt(base);
        }
    }

    // 計算出 R 矩陣 (R = q_transpose*A)
    float r[4][4];

    for (int a = 0; a < 4; a++) {
        for (int b = 0; b < 4; b++) {
            q_transpose[a][b] = q[b][a];
        }
    }

    multiply(q, matrix);

    memcpy(r, global, 4 * 4 * sizeof(float));

    multiply(r, q_transpose);

    qr_decomposition(global, temp + 1);
    
    */

    float matrix_temp[4][4]; // non change original matrix
    copyMatrix(matrix_temp, matrix);

    // 找出特徵向量，先遍歷四個特徵值
    for (int i = 3; i < 4; i++) {
        cout << "EigenValue : " << q[i][i] << endl;
        cout << endl;

        // 遍歷整個矩陣，計算(A+I)
        cout << "A-λI" << endl;
        copyMatrix(matrix, matrix_temp);
        for (int j = 0; j < 4; j++) {
            matrix[j][j] -= q[i][i];
        }
        printMatrix(matrix);
        
        // 以高斯消去法求 homogeneous 線性方程組
        // 先把不為零的首項都變一
        for (int row = 0; row < 4; row++) {
            int j = 0;
            while (j < 4) {
                if (matrix[row][j] != 0) {
                    for (int k = 3; k >= 0; k--) {
                        matrix[row][k] /= matrix[row][j];
                    }
                    break;
                }
                j++;
            }
        }

        cout << endl;
        cout <<"首項化為1"<< endl;
        printMatrix(matrix);
        //進行高斯消去，從底部開始
        for (int row = 3; row >= 0; row--) {
            if (!matrix[row][row])continue;
            // 剩下幾行
            for (int j = row - 1; j >= 0; j--) {
                float mag = matrix[j][row];
                // 幾個元素要處理
                for (int k = row ; k < 4; k++) {
                    matrix[j][k] -= matrix[row][k]*mag;
                }
            }
        }
        // 找特徵矩陣，先求最大交集是哪一個
        float eigenVector[5] = {0,0,0,0,0};

        // 去掉一定為零行(方便程式而已)
        for (int row = 0; row < 4; row++) {
            int count = 0;
            for (int col = 0; col < 4; col++) {
                if (matrix[row][col] == 0) {
                    count++;
                }
            }
            if (count == 3) {
                for (int j = 0; j < 4; j++) {
                    matrix[row][j] = 0;
                }
            }
        }

        cout << endl;
        cout << "去掉一定為零行(方便程式而已)" << endl;
        printMatrix(matrix);

        // 遍歷列，出現最多次非零值的其特徵向量為1
        int most = 0;
        int last_time = 0;
        for (int col = 0; col < 4; col++) {
            int count = 0;
            for (int j = 0; j < 4; j++) {
                if (matrix[j][col] != 0) {
                    count++;
                }
            }
            if (count == 0) continue;
            if (count == last_time and most==0) {
                most = col;
                last_time = count;
            }
            else if (count > last_time) {
                most = col;
                last_time = count;
            }
        }
        eigenVector[most+1] = 1;

        cout << endl;
        cout << "找到設為 x = 1 的列" << endl;
        cout << most << endl;
        
        // 遍歷行，若其主對角線值存在則計算其特徵向量值，不存在(或全零)則跳過
        for (int row = 0; row < 4; row++) {
            int count = 0;
            for (int col = 0; col < 4; col++) {
                if (matrix[row][col] == 0) {
                    count++;
                }
            }
            if (count == 4)continue;
            
            eigenVector[row + 1] = -(matrix[row][most]);
        }

        cout << endl;
        cout << "主對角線值不為零則計算其特徵向量值" << endl;
        printMatrix(matrix);

        // 打印特徵向量
        cout << "EigenVector : [";
        for (int j = 1; j < 5; j++) {
            if (j == 4) {
                cout << eigenVector[j] << "]" << endl;
                break;
            }
            cout << eigenVector[j] << ",";
        }
        cout << "================" << endl;
    }

}