/*
 * check_periodic_lu1.c
 */

#include <stdio.h>
#include <matrix.h>

main()
{
    int i, j, n;
    vector x, b;
    vector al, ad, au, ab, ar;
    matrix a, u, L, Lu;
    int debug = 1;

    printf("n="); scanf("%d", &n);

    /* ベクトルの準備 */
    x = new_vector(n + 1);
    b = new_vector(n + 1);

    /* 行列の準備 */
    al = new_vector(n + 1); ad = new_vector(n + 1); au = new_vector(n + 1);
    ab = new_vector(n + 1); ar = new_vector(n + 1);

    a = new_matrix(n + 1, n + 1);
    L = new_matrix(n + 1, n + 1);
    u = new_matrix(n + 1, n + 1);
    Lu = new_matrix(n + 1, n + 1);

    clear_vector(n, al); clear_vector(n, ad); clear_vector(n, au);
    clear_vector(n, ab); clear_vector(n, ar);

    clear_matrix(n, L);
    clear_matrix(n, u);

    /* A を作る */
    for (i = 1; i <= n; i++) {
	ad[i] = 3.0;
	if (i != 1) al[i] = -1.0;
	if (i != n) au[i] = -1.0;
    }
    ab[1] = -1.0; ab[n - 1] = -1.0;
    ar[1] = -1.0; ar[n - 1] = -1.0;

    /* A を表示する */
    if (debug) {
	printf("A=\n");
	show_vector(n, al); show_vector(n, ad); show_vector(n, au);
	show_vector(n, ab); show_vector(n, ar);
    }
    printf("A=\n");
    normal_matrix(n, a, al, ad, au, ab, ar);
    show_matrix(n, a);

    /* 連立1次方程式の右辺を作る */
    for (i = 1; i <= n; i++)
        x[i] = i;
    multiply_mv(n, b, a, x);

    /* A を LU 分解する */
    ptrilu1(n, al, ad, au, ab, ar);

    /* A を LU 分解した結果を表示する */
    printf("L & U=\n");
    normal_matrix(n, a, al, ad, au, ab, ar);
    show_matrix(n, a);

    /* L, U を二つの matrix に分ける */
    divide_to_lu(n, a, L, u);

    /* L を表示する */
    printf("\nL=\n"); show_matrix(n, L);
    /* U を表示する */
    printf("\nU=\n"); show_matrix(n, u);

    /* L と U をかける */
    multiply_matrix(n, Lu, L, u);

    /* L U を表示する (この結果が A と等しくなればよい) */
    printf("\nL U=\n"); show_matrix(n, Lu);

    /* 連立1次方程式を解く */
    ptrisol1(n, al, ad, au, ab, ar, b);
    /* 解を表示する */
    printf("x(解いたもの)=\n");
    show_vector(n, b);
    printf("x(本当の答え)=\n");
    show_vector(n, x);
 
    return 0;
}

clear_matrix(int n, matrix a)
{
    int i, j;
    for (i = 1; i <= n; i++)
	for (j = 1; j <= n; j++)
	    a[i][j] = 0.0;
}

clear_vector(int n, vector a)
{
    int i;
    for (i = 1; i <= n; i++)
	a[i] = 0.0;
}

show_matrix(int n, matrix a)
{
    int i, j;
    for (i = 1; i <= n; i++) {
	for (j = 1; j <= n; j++)
	    printf("%5.2f ", a[i][j]);
	printf("\n");
    }
}

show_vector(int n, vector a)
{
    int i;
    for (i = 1; i <= n; i++) {
	printf("%5.2f ", a[i]);
    }
    printf("\n");
}

multiply_mv(int n, vector ab, matrix a, vector b)
{
    int i, j;
    double s;
    for (i = 1; i <= n; i++) {
        s = 0.0;
        for (j = 1; j <= n; j++)
            s += a[i][j] * b[j];
        ab[i] = s;
    }
}

multiply_matrix(int n, matrix ab, matrix a, matrix b)
{
    int i, j, k;
    double s;
    for (i = 1; i <= n; i++)
	for (j = 1; j <= n; j++) {
	    s = 0.0;
	    for (k = 1; k <= n; k++)
		s += a[i][k] * b[k][j];
	    ab[i][j] = s;
	}
}

normal_matrix(int n, matrix a,
	      vector al, vector ad, vector au, vector b, vector r)
{
    int i, nm1;
    nm1 = n - 1;
    clear_matrix(n, a);
    for (i = 1; i <= n; i++) {
	a[i][i] = ad[i];
	if (i != 1)
	    a[i][i - 1] = al[i];
	if (i != n)
	    a[i][i + 1] = au[i];
    }
    for (i = 1; i < nm1; i++) {
	a[n][i] = b[i];
	a[i][n] = r[i];
    }
}

divide_to_lu(int n, matrix a, matrix L, matrix u)
{
    int i, j;
    clear_matrix(n, L);
    clear_matrix(n, u);
    for (i = 1; i <= n; i++) {
	for (j = 1; j <= n; j++)
	    if (i < j)
		u[i][j] = a[i][j];
	    else if (i > j)
		L[i][j] = a[i][j];
	    else {
		u[i][j] = a[i][j];
		L[i][j] = 1.0;
	    }
    }
}