#include < stdio.h >
#include < math.h >
#define N 3
#define FO "%-10.5g"
void multm(float * a, float * b, int m, int n, int k, float * c) {
int i,
j,
l,
u;
for (i = 0; i < m; i++) for (j = 0; j < k; j++) {
u = i * k + j;
c = 0.0;
for (l = 0; l < n; l++) c += a * b[l * k + j];
}
}
void scanfm(float * a, int m, int n, char aa) {
int i,
j;
if (n == 1) for (i = 0; i < m; i++) {
printf("%c[%d]=", aa, i + 1);
scanf("%f", a + i);
} else {
for (i = 0; i < m; i++) for (j = 0; j < n; j++) {
printf("%c[%d][%d]=", aa, i + 1, j + 1);
scanf("%f", a + i * n + j);
}
}
}
void printm(float * a, int m, int n, char aa) {
int i,
j,
k;
printf("%c=\n", aa);
for (i = 0; i < m; i++) {
for (j = 0; j < n; j++) {
k = i * n + j;
if (fabs(a[k]) <= 1e-6) a[k] = 0.0;
printf(FO, a[k]);
}
printf("\n");
}
}
void copym(float * a, float * b, int m, int n) {
int i,
j,
u;
for (i = 0; i < m; i++) for (j = 0; j < n; j++) {
u = i * n + j;
b = a;
}
}
float fanshu(float * a, int m, int n) {
int i,
j,
u;
float r = 0.0,
r1;
for (i = 0; i < m; i++) {
r1 = 0.0;
for (j = 0; j < n; j++) {
u = i * n + j;
r1 += fabs(a);
}
if (r1 > r) r = r1;
}
return r;
}
int Jacobi(float * a, int n, float * f, float * x, float eps, float * xk) {
float r,
r1,
r2;
int i,
j,
u;
for (i = 0; i < n; i++) {
u = i * n + i;
if (fabs(a) <= 1e-5) return 0;
for (j = 0; j < n; j++) if (j != i) a /= -a;
f /= a;
a = 0.0;
}
r1 = fanshu(a, n, n);
if (r1 >= 1) return 0;
r1 = r1 / (1.0 - r1);
r = 1.0;
for (i = 0; i < n; i++) xk = 0.0;
while (r > eps) {
for (i = 0; i < n; i++) {
r2 = 0.0;
for (j = 0; j < n; j++) if (j != i) r2 += a * xk[j];
x = r2 + f;
}
for (i = 0; i < n; i++) xk = x - xk;
r = r1 * fanshu(xk, n, 1);
for (i = 0; i < n; i++) xk = x;
}
return 1;
}
int JGS(float * a, int n, float * f, float * x, float eps, float * xk, float * L, float * U) {
float r,
r1,
r2;
int i,
j,
u;
for (i = 0; i < n; i++) {
u = i * n + i;
if (fabs(a) <= 1e-5) return 0;
for (j = 0; j < n; j++) if (j != i) a /= -a;
f /= a;
a = 0.0;
}
for (i = 0; i < n; i++) {
L = 1.0;
for (j = 0; j < n; j++) {
u = i * n + j;
if (j < i) L = a;
else if (j > i) U = a;
}
}
multm(L, U, n, n, n, a);
r1 = fanshu(a, n, n);
if (r1 >= 1) return 0;
r1 = r1 / (1.0 - r1);
r = 1.0;
for (i = 0; i < n; i++) xk = 0.0;
while (r > eps) {
for (i = 0; i < n; i++) {
r2 = 0.0;
for (j = 0; j < n; j++) {
u = i * n + j;
if (j > i) r2 += U * xk[j];
else if (j < i) r2 += L * x[j];
}
x = r2 + f;
}
for (i = 0; i < n; i++) xk = x - xk;
r = r1 * fanshu(xk, n, 1);
for (i = 0; i < n; i++) xk = x;
}
return 1;
}
void main() {
static float A[N][N],
A1[N][N],
L[N][N],
U[N][N],
f[N],
f1[N],
x[N],
xk[N];
float eps = 1e-6;
int i;
scanfm(A[0], N, N, 'A');
scanfm(f, N, 1, 'b');
printm(A[0], N, N, 'A');
printm(f, N, 1, 'b');
copym(A[0], A1[0], N, N);
copym(f, f1, N, 1);
i = Jacobi(A[0], N, f, x, eps, xk);
if (i == 0) printf("Jacobi Failed!\n");
else {
printf("Jacobi:\n");
printm(x, N, 1, 'x');
}
i = JGS(A1[0], N, f1, x, eps, xk, L[0], U[0]);
if (i == 0) printf("JGS Failed!\n");
else {
printf("JGS:\n");
printm(x, N, 1, 'x');
}
}
标签:xk,Jacbi,r1,r2,++,float,int From: https://blog.51cto.com/u_16076050/6195958