blob: aca65e8ff329e0f8042d1e6292f2a730a98b0d6e (
plain) (
blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
|
#include "sle.h"
#include "stdio.h"
#include <stdlib.h>
#include <math.h>
/*
* Gaussian elimination
*/
int sle_solve(double *a, size_t n, double *x) {
int i, j, k, index;
double max, t;
const int col = n + 1;
const float eps = 0.00001;
if (n == 0) {
return - 1;
}
for (k = 0; k < n; k++) {
max = fabs(a[col*k + k]);
index = k;
for (i = k + 1; i < n; i++) {
t = fabs(a[col * i + k]);
if (t > max) {
max = t;
index = i;
}
}
if (max < eps) {
return -1;
}
if (index != k) {
for (i = 0; i < col; i++) {
t = a[col * k + i];
a[col * k + i] = a[col * index + i];
a[col * index + i] = t;
}
}
for (i = k; i < n; i++) {
t = a[col * i + k];
if (fabs(t) < eps)
continue;
for (j = 0; j < col; j++) {
a[i * col + j] /= t;
}
if (i != k) {
for (j = 0; j < col; j++) {
a[i * col + j] -= a[k * col + j];
}
}
}
}
for (k = n - 1; k >= 0; k--) {
x[k] = a[col * k + n];
for (i = 0; i < k; i++) {
a[col * i + n] -= a[col * i + k] * x[k];
}
}
return 0;
}
|
