/*
 * testgauss1.c --- gauss1.c の中の関数gauss1() の使い方
 *
 * コンパイルは
 *   gcc -c matutil.c
 *   gcc -c gauss1.c
 *   gcc -o testgauss1 testgauss1.c gauss1.o matutil.o
 *
 *   2 x + 3 y -   z = 5
 *   4 x + 4 y - 3 z = 3
 *  -2 x + 3 y -   z = 1
 *
 * の解が(x,y,z)=(1,2,3) であることで確認。
 */

#include <stdio.h>
#include "matutil.h"

/* 次の二つの関数は testmatmul1.c の中と同じ*/
double dotprod1(int n, vector x, vector y)
{
   int i;
   double s = 0;
   for (i = 1; i <= n; i++)
     s += x[i] * y[i];
   return s;
}

/* m 行 n 列の行列 A と n 次元ベクトル b の積を c に代入する*/
void mul_matrix_vector1(int m, int n, vector c, matrix a, vector b)
{
  int i;
  /* c の第 i 成分は A の第 i 行ベクトルと b の内積*/
  for (i = 1; i <= m; i++)
    c[i] = dotprod1(n, a[i], b);
}

/* ベクトル x を表示する*/
void show_vector1(int n, vector x)
{
  int i;
  for (i = 1; i <= n; i++)
    printf("%g ", x[i]);
  printf("\n");
}

/* 行列 A を表示する*/
void show_matrix1(int m, int n, matrix A)
{
  int i;
  for (i = 1; i <= m; i++)
    show_vector1(n, A[i]);
}

main()
{
  int i, j, n;
  matrix a;
  vector x, b;
  n = 3;
  /* 行列とベクトルの準備 */
  a = new_matrix(n+1, n+1);
  x = new_vector(n+1);
  b = new_vector(n+1);
  /* 真の解を先に決める*/
  x[1] = 1; x[2] = 2; x[3] = 3;
  /* 係数行列を決める */
  a[1][1] =  2; a[1][2] = 3; a[1][3] = -1;
  a[2][1] =  4; a[2][2] = 4; a[2][3] = -3;
  a[3][1] = -2; a[3][2] = 3; a[3][3] = -1;
  /* 右辺の b を計算する */
  mul_matrix_vector1(n, n, b, a, x);
  /* 行列 A と右辺 b を計算する */
  printf("A=\n");
  show_matrix1(n, n, a);
  printf("b=\n");
  show_vector1(n, b);
  /* Gauss の消去法で解く */
  gauss1(n, a, b);
  /* 真の解と計算で求めた解を表示する */
  printf("真の解=\n");
  show_vector1(n, x);
  printf("計算解=\n");
  show_vector1(n, b);
}