// JGS_hinagata_utf8.c
// 連立一次方程式の反復解法: 
// ヤコビ法, ガウス-ザイデル法
// 雛形プログラム

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

////////// マクロと大域変数の宣言 //////////////////////////////////////////////

#define N 128    // 行列のサイズの上限
int     size;    // 行列のサイズ

////////// ヤコビ法･ガウス-ザイデル法の関数定義 ////////////////////////////////

// ヤコビ法で近似ベクトル x を更新する関数
void Jacobi(double a[][N], double b[], double x[])
{

    //////////////////////////////////////////////////////
    ////// 課題：                                    /////
    ////// 次のガウス-ザイデル法をまねてここを埋めよ /////
    //////////////////////////////////////////////////////

    return;
}

// ガウス-ザイデル法で近似ベクトル x を更新する関数
void GaussSeidel(double a[][N], double b[], double x[])
{
    int i, j;
    double w[N];     // w = 新x
    double z;
    for(i = 0; i < size; i++){
        z = b[i];
        for(j = 0; j < i; j++) z -= a[i][j] * w[j];
        for(j = i + 1; j < size; j++) z -= a[i][j] * x[j];
        w[i] = z / a[i][i];
    }
    for(i = 0; i < size; i++) x[i] = w[i];    // w を x に上書き
    return;
}

////////// その他の関数のプロトタイプ宣言 //////////////////////////////////////

// 行列 a を出力する関数
void matprint(double a[][N]);
// ベクトル v を出力する関数
void vecprint(double v[]);        
// 行列 a とベクトル v の積をベクトル w に格納する関数
void matvecmul(double a[][N], double v[], double w[]);
// ベクトル u とベクトル v の差 u - v をベクトル w に格納する関数
void vecsub(double u[], double v[], double w[]);
// ベクトル v の長さを返す関数
double norm(double v[]);

////////// メインルーティン ////////////////////////////////////////////////////

int main()
{
    double a[N][N];     // 行列 A を表す２次元配列
    double b[N];        // ベクトル b を表す配列
    double x[N];        // 解ベクトル x を表す配列
    double ax[N];       // ベクトル Ax を格納する配列
    double e[N];        // 誤差ベクトル Ax - b を格納する配列
    double enorm;       // 誤差ベクトル Ax - b の長さを格納する変数
    int maxrec = 128;   // 最大反復回数
    double adm = 1e-10; // 許容誤差
    int i, j, step, mode, trace;

    ////////// 方程式の設定 //////////

    printf("行列を設定してください：\n");
    printf("1 サンプルデータ\n");
    printf("2 ランダムな行列\n");
    printf("3 課題の行列\n");
    scanf("%d", &mode);

    switch(mode){
        case 1:
            // サンプルデータ
            size = 2;
            a[0][0] = 3; a[0][1] = -2; a[1][0] = 1; a[1][1] = 3;
            b[0] = 1; b[1] = 4;
            break;

        case 2:
            // ランダムな行列
            srand(time(NULL));    // 乱数の初期化
            printf("ランダムな行列を生成します。\n");
            printf("行列のサイズを入力してください。\n");
            scanf("%d", &size);
            for(i = 0; i < size; i++){
                for(j = 0; j < size; j++)
                    if(i == j) a[i][i] = 20 * ((1.0 * rand()) / RAND_MAX - 0.5);
                    else  a[i][j] = (20.0 / size) * ((1.0 * rand()) / RAND_MAX - 0.5);    // 非対角要素は小さく
                b[i] = 1.0;
            }
            break;

        case 3:
            // 課題のヒルベルト行列
            printf("課題の方程式で計算します。\n");
            printf("行列のサイズを入力してください。\n");
            scanf("%d", &size);
            for(i = 0; i < size; i++){
                for(j = 0; j < size; j++) a[i][j] = 0.;
                if(i > 0) a[i][i - 1] =  1.;
                a[i][i] = 3.;
                if(i < size - 1) a[i][i + 1] = 1.;
                b[i] = 1.;
            }
            break;
    }

    printf("\n----- A と b の設定 -----\n\n");
    printf("A:\n"); matprint(a); printf("\n");
    printf("b:\n"); vecprint(b); printf("\n");

    printf("途中経過を出力しますか ( Yes:1, No:0 )\n");
    scanf("%d", &trace);

    ////////// Jacobi 法 実行部 //////////

    printf("\n----- Jacobi 法 -----\n\n");
    for(i = 0; i < size; i++) x[i] = 0.;
    step = 0;
    enorm = 1.;    // 誤差の初期値はダミー
    while(enorm > adm && step < maxrec){
        if(trace){ printf("%d-th step :\n\n", step); vecprint(x); printf("\n"); }
        Jacobi(a, b, x);        // 近似ベクトル x の更新
        matvecmul(a, x, ax);    // ax = Ax
        vecsub(ax, b, e);        // e  = Ax - b
        enorm = norm(e);        // enorm = (e の長さ)
        step++;
    }
    printf("%d-th step :\n\n", step); vecprint(x); printf("\n");
    printf("検算 Ax - b :\n\n");
    vecprint(e); printf("\n");
    if(enorm <= adm)
        printf("ヤコビ法では %d ステップ掛かりました。\n\n", step);
    else
        printf("ヤコビ法では収束しませんでした。\n\n", step);
    getchar(); printf("リターンキーを押してください\n"); getchar();
    
    ////////// Gauss-Seidel 法 実行部 //////////

    printf("----- Gauss-Seidel 法 -----\n\n");
    for(i = 0; i < size; i++) x[i] = 0.;
    step = 0;
    enorm = 1.;    // 誤差の初期値はダミー
    while(enorm > adm && step < maxrec){
        if(trace){ printf("%d-th step :\n\n", step); vecprint(x); printf("\n"); }
        GaussSeidel(a, b, x);    // 近似ベクトル x の更新
        matvecmul(a, x, ax);    // ax = Ax
        vecsub(ax, b, e);        // e  = Ax - b
        enorm = norm(e);        // enorm = (e の長さ)
        step++;
    }
    printf("%d-th step :\n\n", step); vecprint(x); printf("\n");
    printf("検算 Ax - b :\n\n");
    vecprint(e); printf("\n");
    if(enorm <= adm)
        printf("ガウス・ザイデル法では %d ステップ掛かりました。\n\n", step);
    else
        printf("ガウス・ザイデル法では収束しませんでした。\n\n", step);

    return 0;
}

////////// 関数定義（以下は触らなくてよい） ////////////////////////////////////

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

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

// 行列 a とベクトル v の積をベクトル w に格納する関数
void matvecmul(double a[][N], double v[], double w[])
{
    int i, j;
    double z;
    for(i = 0; i < size; i++){
        z = 0;
        for(j = 0; j < size; j++)
            z += a[i][j] * v[j];
        w[i] = z;
    }
    return;
}

// ベクトル u とベクトル v の差 u - v をベクトル w に格納する関数
void vecsub(double u[], double v[], double w[])
{
    int i;
    for(i = 0; i < size; i++) w[i] = u[i] - v[i];
    return;
}

// ベクトル v の長さを返す関数
double norm(double v[])
{
    int i;
    double z = 0.;
    for(i = 0; i < size; i++) z += v[i] * v[i];
    return sqrt(z);
}