レポート課題

必須課題１）
LU分解モジュールの未実装関数（lu_solve()）を実装し，LU分解を用いて連立方程式を解くアプリケーションを完成させる．

必須課題２）

オプション課題１）
LとUとxのために領域を確保しなくても良い場合もある．
そこで，係数行列AのLU分解後，Aの領域をLUで上書きし，それを用いて，連立方程式を解き，右辺ベクトルbの領域を解ベクトルxを上書きするようにして，メモリ領域を節約し，かつ冗長な計算を省くことが可能である．
必須課題１）で作成したLU分解モジュール（lu.h，lu.c）のlu_decomp()とlu_solve()を以下のlu_mod.hとlu_mod.cを参考にlu_decomp_mod()とlu_solve_mod()として作成する
テストのためのmain関数としてlu_mod_test.cを利用すること．

ソリューション・エクスプローラー構成


lu_mod.h

#ifndef LU_H
#define	LU_H

#include "vector.h"
#include "matrix.h"

/**
 * @brief LU分解
 * @param[in,out] A LU分解対象行列，LU分解後の行列
 * @return 行列式
 * @attention LU分解後，Aの非対角成分の下三角成分はLで，Aの非対角成分の上三角成分はUで，Aの対角成分はLの対角成分で上書きされる
 */
double lu_decomp_mod(matrix* A);

/**
 * @brief LU分解を用いて連立方程式を解く
 * @param[in,out] b 右辺ベクトル，解ベクトル
 * @param[in] A lu_decomp_mod後のLU分解行列
 * @attention bは解ベクトルで上書きされる
 */
void lu_solve_mod(vector* b, const matrix* A);

#endif	/* LU_H */

lu_mod.c

#include "lu_mod.h"

double lu_decomp_mod(matrix* A) {
    /* ここを実装せよ */
}

void lu_solve_mod(vector* b, const matrix* A) {
    /* ここを実装せよ */
}

lu_mod_test.c

#include <stdlib.h>
#include "vector.h"
#include "matrix.h"
#include "lu_mod.h"

/**
 *     / 8 16 24  32 \
 * A = | 2  7 12  17 |
 *     | 6 17 32  59 |
 *     \ 7 22 46 105 /
 * 
 *     / 160 \
 * b = |  70 |
 *     | 198 |
 *     \ 291 /
 * @param[out] A 係数行列
 * @param[out] b 右辺ベクトル
 */
void set_equation1(matrix* A, vector* b) {
    /* 行列サイズ指定 */
    int N = 4;
    /* 領域確保 */
    matrix_new(N, N, A);
    vector_new(N, b);
    /* 係数行列要素指定 */
    A->val[0][0] = 8;
    A->val[0][1] = 16;
    A->val[0][2] = 24;
    A->val[0][3] = 32;
    A->val[1][0] = 2;
    A->val[1][1] = 7;
    A->val[1][2] = 12;
    A->val[1][3] = 17;
    A->val[2][0] = 6;
    A->val[2][1] = 17;
    A->val[2][2] = 32;
    A->val[2][3] = 59;
    A->val[3][0] = 7;
    A->val[3][1] = 22;
    A->val[3][2] = 46;
    A->val[3][3] = 105;
    /* 右辺ベクトル要素指定 */
    b->val[0] = 160;
    b->val[1] = 70;
    b->val[2] = 198;
    b->val[3] = 291;
}

/**
 *     / 3 2 1 \
 * A = | 2 4 3 |
 *     \ 1 3 5 /
 * 
 *     /  6   \
 * b = | 10   |
 *     \ 12.5 /
 * @param[out] A 係数行列
 * @param[out] b 右辺ベクトル
 */
void set_equation2(matrix* A, vector* b) {
    /* 行列サイズ指定 */
    int N = 3;
    /* 領域確保 */
    matrix_new(N, N, A);
    vector_new(N, b);
    /* 係数行列要素指定 */
    A->val[0][0] = 3;
    A->val[0][1] = 2;
    A->val[0][2] = 1;
    A->val[1][0] = 2;
    A->val[1][1] = 4;
    A->val[1][2] = 3;
    A->val[2][0] = 1;
    A->val[2][1] = 3;
    A->val[2][2] = 5;
    /* 右辺ベクトル要素指定 */
    b->val[0] = 6;
    b->val[1] = 10;
    b->val[2] = 12.5;
}

/**
 *     / 1 2 2 \
 * A = | 3 6 2 |
 *     \ 1 1 1 /
 * 
 *     / 1 \
 * b = | 7 |
 *     \ 5 /
 * @param[out] A 係数行列
 * @param[out] b 右辺ベクトル
 */
void set_equation3(matrix* A, vector* b) {
    /* 行列サイズ指定 */
    int N = 3;
    /* 領域確保 */
    matrix_new(N, N, A);
    vector_new(N, b);
    /* 係数行列要素指定 */
    A->val[0][0] = 1;
    A->val[0][1] = 2;
    A->val[0][2] = 2;
    A->val[1][0] = 3;
    A->val[1][1] = 6;
    A->val[1][2] = 2;
    A->val[2][0] = 1;
    A->val[2][1] = 1;
    A->val[2][2] = 1;
    /* 右辺ベクトル要素指定 */
    b->val[0] = 1;
    b->val[1] = 7;
    b->val[2] = 5;
}

int main(void) {
    double det; /* 行列式 */
    matrix A; /* 係数行列 */
    vector b; /* 右辺ベクトル */
    char sw[2]; /* 方程式切替スイッチ */

    printf("input 1,2,3 ? ");
    fgets(sw, sizeof (sw), stdin);
    switch (atoi(sw)) {
        case 1:
            set_equation1(&A, &b);
            break;
        case 2:
            set_equation2(&A, &b);
            break;
        case 3:
            set_equation3(&A, &b);
            break;
        default:
            return EXIT_FAILURE;
    }
    /* 係数行列表示 */
    printf("A =\n");
    matrix_print(&A, "%f");
    /* 右辺ベクトル表示 */
    printf("b = ");
    vector_print(&b, "%f");
    /* LU分解 */
    det = lu_decomp_mod(&A);
    /* LU分解後の係数行列を表示 */
    printf("A =\n");
    matrix_print(&A, "%f");
    /* 行列式表示 */
    printf("det = %3.0f\n", det);
    /* LU分解を用いて連立方程式を解く */
    lu_solve_mod(&b, &A);
    /* 解ベクトル表示 */
    printf("b = ");
    vector_print(&b, "%f");
    /* 領域解放 */
    matrix_delete(&A);
    vector_delete(&b);

    return EXIT_SUCCESS;
}


input 1,2,3 ? 1
A =
[8.000000, 16.000000, 24.000000, 32.000000
2.000000, 7.000000, 12.000000, 17.000000
6.000000, 17.000000, 32.000000, 59.000000
7.000000, 22.000000, 46.000000, 105.000000]
b = [160.000000, 70.000000, 198.000000, 291.000000]
A =
[8.000000, 2.000000, 3.000000, 4.000000
2.000000, 3.000000, 2.000000, 3.000000
6.000000, 5.000000, 4.000000, 5.000000
7.000000, 8.000000, 9.000000, 8.000000]
det = 768
b = [4.000000, 3.000000, 2.000000, 1.000000]

レポート提出