// rep09hinagata.c 
// LU 分解法 雛形プログラム
//
// 140行目付近の課題その１、150行目付近の課題その２を完成せよ

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>

// 行列のサイズは大域変数としてここで設定する。
int    size = 16;

#define N 128                                          // 行列の最大サイズ
#define RANDOM(n) (rand()%(2*n+1)-n)                   // 絶対値 |n| 以下の整数乱数 
#define REALRANDOM(x) (x*((2.0*rand()/RAND_MAX)-1))    // 絶対値 |x| 以下の実数乱数 

///// プロトタイプ宣言 /////

// 行列 A, ベクトル b の初期設定
void initialize(double a[][N], double b[]);    
// LU分解法の Step 1（LU分解）
void LUstep1(double a[][N], double ell[][N], double u[][N]);
// LU分解法の Step 2（前進代入）
void LUstep2(double ell[][N], double y[], double b[]);
// LU分解法の Step 3（後退代入）
void LUstep3(double u[][N], double x[], double y[]);
// 行列の積を計算する関数
void matmul(double a[][N], double b[][N], double c[][N]);
// 行列と列ベクトルの積を計算する関数
void matvecmul(double a[][N], double b[], double c[]);
// 行列を出力する関数
void matprint(double a[][N]);
// 列ベクトルを出力する関数
void vecprint(double b[]);

///// メイン関数 /////

int main()
{
    double a[N][N];     // 行列 A
    double ell[N][N];   // 行列 L
    double u[N][N];     // 行列 U
    double c[N][N];     // 検算のための行列
    double b[N];        // ベクトル b
    double x[N];        // 解ベクトル x
    double y[N];        // 前進代入で求めるベクトル y
    double z[N];        // 検算のためのベクトル
    double w;           // 総和を求めるための変数
    int i, j, k;        // for文の制御変数
    double e, f;        // 誤差測定のための変数

    ///// 初期設定の出力 /////
    initialize(a, b);
    printf("----- A と b の設定 -----\n\n");
    printf("A:\n"); matprint(a); printf("\n");
    printf("b:\n"); vecprint(b); printf("\n");

    ///// LU分解 /////
    LUstep1(a, ell, u);
    // L, U を出力
    printf("----- LU分解 -----\n\n");
    printf("L:\n"); matprint(ell); printf("\n");
    printf("U:\n"); matprint(u); printf("\n");

    // 検算のため積 C = LU を計算
    matmul(ell, u, c);
    printf("----- 検算-----\n\n");
    printf("LU:\n"); matprint(c); printf("\n");
    printf("A:\n"); matprint(a); printf("\n");

    ///// 前進代入:  Ly = b の解 y の計算 /////
    LUstep2(ell, y, b);
    printf("----- 前進代入 -----\n\n");
    printf("y:\n"); vecprint(y); printf("\n");

    ///// 後退代入:  Ux = y の解 x の計算 /////
    LUstep3(u, x, y);
    printf("----- 後退代入 -----\n\n");
    printf("x:\n"); vecprint(x); printf("\n");

    // 検算のため積 z = Ax を計算
    matvecmul(a, x, z);
    printf("----- 検算 -----\n\n");
    printf("Ax:\n"); vecprint(z); printf("\n");
    printf("b:\n"); vecprint(b); printf("\n");

    e = 0.0;
    for(i = 0; i < size; i++){
        f = fabs(z[i] - b[i]);
        if(f > e) e = f;
    }
    printf("最大誤差 = %12.10lf\n\n", e);
    return 0;
}

////////// 関数定義部 //////////

// 行列 A, ベクトル b の初期設定
void initialize(double a[][N], double b[N])
{
    int i, j;
    srand(time(NULL));    // 乱数の初期化

    ///// 課題 /////
    for(i = 0; i < size; i++) for(j = 0; j < size; j++) a[i][j] = 1.0 / (i + j + 1);
    for(i = 0; i < size; i++) b[i] = 1.0;

    /*
    ///// ランダムな行列・ベクトル /////
    for(i = 0; i < size; i++) for(j = 0; j < size; j++) a[i][j] = REALRANDOM(10);
    for(i = 0; i < size; i++) b[i] = REALRANDOM(10);

    ///// ランダムな整数行列・ベクトル /////
    for(i = 0; i < size; i++) for(j = 0; j < size; j++) a[i][j] = RANDOM(10);
    for(i = 0; i < size; i++) b[i] = RANDOM(10);
    */

    return;
}

// LU分解法の Step 1（LU 分解）
// A = LU の形に分解する
void LUstep1(double a[][N], double ell[][N], double u[][N])
{
    int i, j;

    // L, U の初期値は単位行列
    for(i = 0; i < size; i++){
        for(j = 0; j < size; j++){
            ell[i][j] = 0.0;
            u[i][j] = 0.0;
        }
        ell[i][i] = 1.0;
        u[i][i] = 1.0;
    }

    // LU 分解

    //////////////////////////////////////////////
    /////     課題その１：ここを完成せよ     /////
    //////////////////////////////////////////////

    return;
}

// LU分解法の Step 2（前進代入）
// Ly = b の解 y の計算
void LUstep2(double ell[][N], double y[], double b[])
{
    //////////////////////////////////////////////
    /////     課題その２：ここを完成せよ     /////
    //////////////////////////////////////////////

    return;
}

// LU分解法の Step 3（後退代入）
// Ux = y の解 x の計算
void LUstep3(double u[][N], double x[], double y[])
{
    int i, j;
    double w;
    for(i = size - 1; i >= 0; i--){
        w = 0.0;
        for(j = i + 1; j < size; j++) w += u[i][j] * x[j];
        x[i] = y[i] - w;
    }
    return;
}

// 行列の積を計算する関数
// C = AB
void matmul(double a[][N], double b[][N], double c[][N])
{
    int i, j, k;
    double x;
    for(i = 0; i < size; i++){
        for(j = 0; j < size; j++){
            x = 0.0;
            for(k = 0; k < size; k++) x += a[i][k] * b[k][j];
            c[i][j] = x;
        }
    }
    return;
}

// 行列と列ベクトルの積を計算する関数
// c = Ab
void matvecmul(double a[][N], double b[], double c[])
{
    int i, j;
    double w;
    for(i = 0; i < size; i++){
        w = 0.0;
        for(j = 0; j < size; j++) w += a[i][j] * b[j];
        c[i] = w;
    }
    return;
}

// 行列を出力する関数
void matprint(double a[][N])
{
    int i, j;
    for(i = 0; i < size; i++){
        for(j = 0; j < size; j++)
            printf("%12.7lf ", a[i][j]);
        printf("\n");
    }
    return;
}

// 列ベクトルを出力する関数
void vecprint(double b[])
{
    int i;
    for(i = 0; i < size; i++)
        printf("%15.10lf\n", b[i]);
    return;
}