Jacbi_J
2015-09-17 11:46
405 查看
#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 , A1 , L , U , f , f1 , x , xk ; 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'); } }
相关文章推荐
- EventBus小实例-传值、控制其他页控件显示
- VB.NET合并图片
- HBase单机模式安装
- Android之socket客户端
- pdf怎么转换成jpg格式欣赏
- java基本数据类型默认值
- 我的第一篇博客
- IE7中图片链接有边框问题解决方法
- 正则表达式
- hibernate4操作数据库blob案例
- Mac下通过Homebrew安装配置MySQL
- c#中动态输入和计算 代数表达式
- cocos2d-x android游戏启动后需要多次退出
- An SSL error has occurred and a secure connection to the server cannot be made
- 2015 Selenium大会
- Dedecms中列表页缩略图怎么大小控制
- android下访问已存在的sqlite数据库文件的办法
- 运维进行时-修改数据库密码
- HDU 5441 Travel(离线 + 带权并查集)——2015 ACM/ICPC Asia Regional Changchun Online
- IEulerRK