// 가우스 소거법 입니다.
// 사용법은 gauss() 함수의 파라메터가 4개 있지요.
// 첫번째가 입력및 출력을 맏은 부분이구요. 두번째가
// 행의 수, 세번째가 열의 수, 네번째가 오차한계입니다.
// 출력은 Call by reference 로 첫번째 파라메터에 들어갑니다.
// 입력은 AX=B 에서 A B행렬을 붙여서 2차원 행렬로
// 넣으면 되고, 출력은 그 배열의 뒤쪽에 들어갑니다.
///////////////////////////////////////////////////////
#include "stdio.h"
#include "math.h"
int gauss(double a[], int row, int col, double eps);
void main()
{
int i, j, row = 3, col = 5;
static double a[3][5] = {
{ 5.0, 1.0,
1.0, 10.0, 18.0 },
{ 1.0, -7.0, 2.0,
-7.0, -9.0 },
{ 1.0, 1.0,
6.0, 21.0, 11.0 }};
double eps = 1.0e-6;
for(i = 0; i < row; i++) {
for(j
= 0; j < col; j++)
printf("%5.0f
|", a[i][j]);
printf("\n");
}
gauss((double *)a, row, col, eps);
printf(" Solution \n");
for(i = 0; i < row; i++) {
for(j
= row; j < col; j++)
printf("
%15.7e", a[i][j]);
printf("\n");
}
}
int gauss(double a[], int row, int col, double eps)
{
int i, j, k, s, t, u, v, iw, r, kp1,
mm1, work[500];
double w, wmax, pivot, api;
if(row < 2 || row > 500 ||
row >= col || eps <= 0.0)
return(999);
for(i = 0; i < row; i++) work[i] = i;
for(k = 0; k < row; k++) {
wmax
= 0.0;
for(i
= k; i < row; i++) {
w
= fabs(a[i*col + k]);
if(w
> wmax) {
wmax
= w;
r
= i;
}
}
u =
r * col;
v =
k * col;
s =
u + k;
pivot
= a[s];
api
= fabs(pivot);
if(api
<= eps) return(1);
if(r
!= k) {
iw
= work[k];
work[k]
= work[r];
work[r]
= iw;
for(j
= 0; j < col; j++) {
s
= v + j;
t
= u + j;
w
= a[s];
a[s]
= a[t];
a[t]
= w;
}
}
for(i
= k; i < col; i++) a[v+i] /= pivot;
kp1
= k + 1;
if(kp1
> row) break;
for(i
= kp1; i < row; i++) {
u
= i * col;
s
= u + k;
w
= a[s];
if(w
!= 0) {
for(j
= kp1; j < col; j++)
a[u+j]
-= w * a[v+j];
}
}
}
mm1 = row -1;
for(k = row; k < col; k++) {
for(i
= 0; i < mm1; i++) {
s
= row - i - 2;
t
= s + 1;
u
= row + t - 1;
r
= s * col;
w
= 0.0;
for(j
= t; j < row; j++) {
v
= u - j;
w
+= a[r+v] * a[v*col+k];
}
a[r+k]
-= w;
}
}
return(0);
}