qyz
来源：互联网 发布：java 图片合成pdf 编辑：程序博客网 时间：2024/05/22 15:43
#include < stdio.h > #include < math.h > #define N1 2#define M1 3#define F(1e-6)#if M1 < N1#define T M#define M N1#define N M1#else#define M M1#define N N1#endif#define STR { if (fabs(a[j]) <= F * 10) a[j] = 0.0; printf("%-10.5g", a[j]);}float A[m][N],U[m][m],V[N][N];int p,q;double arraym(float * x, int k) { int i; double r = 0.0; for (i = 0; i < k; i++) r += x * x; r = sqrt(r); return r;}void mmmn2(float a[m][m], float b[m][N]) { int i, j, k; float t[m][N]; for (i = 0; i < M; i++) for (j = 0; j < N; j++) { t[j] = 0.0; for (k = 0; k < M; k++) t[j] += a[k] * b[k][j]; } for (i = 0; i < M; i++) for (j = 0; j < N; j++) b[j] = t[j];}void mmmm1(float a[m][m], float b[m][m]) { int i, j, k; float t[m][m]; for (i = 0; i < M; i++) for (j = 0; j < M; j++) { t[j] = 0.0; for (k = 0; k < M; k++) t[j] += a[k] * b[k][j]; } for (i = 0; i < M; i++) for (j = 0; j < M; j++) a[j] = t[j];}void mnnn1(float a[m][N], float b[N][N]) { int i, j, k; float t[m][N]; for (i = 0; i < M; i++) for (j = 0; j < N; j++) { t[j] = 0.0; for (k = 0; k < N; k++) t[j] += a[k] * b[k][j]; } for (i = 0; i < M; i++) for (j = 0; j < N; j++) a[j] = t[j];}void nnnn1(float a[N][N], float b[N][N]) { int i, j, k; float t[N][N]; for (i = 0; i < N; i++) for (j = 0; j < N; j++) { t[j] = 0.0; for (k = 0; k < N; k++) t[j] += a[k] * b[k][j]; } for (i = 0; i < N; i++) for (j = 0; j < N; j++) a[j] = t[j];}void Imm(float a[m][m]) { int i, j; for (i = 0; i < M; i++) for (j = 0; j < M; j++) { if (i == j) a[j] = 1.0; else a[j] = 0.0; }}void Inn(float a[N][N]) { int i, j; for (i = 0; i < N; i++) for (j = 0; j < N; j++) { if (i == j) a[j] = 1.0; else a[j] = 0.0; }}float fanshu(void) { int i, j; float r = 0.0; for (i = p; i <= q; i++) for (j = p; j <= q; j++) r += A[j] * A[j]; r = sqrt(r); return r;}float miu(void) { float x[4], g1, g2, det, q2; if (q - p == 1) q2 = 0; else q2 = A[q - 2][q - 1]; x[0] = A[q - 1][q - 1] * A[q - 1][q - 1] + q2 * q2; x[3] = A[q][q] * A[q][q] + A[q - 1][q] * A[q - 1][q]; x[1] = A[q - 1][q - 1] * A[q - 1][q]; x[2] = x[1]; det = (x[0] - x[3]) * (x[0] - x[3]) + 4 * x[1] * x[1]; det = sqrt(det); g1 = (x[0] + x[3] + det) / 2; g2 = (x[0] + x[3] - det) / 2; if (fabs(x[3] - g1) > fabs(x[3] - g2)) g1 = g2; return g1;}void printmm(float a[m][m]) { int i, j; for (i = 0; i < M; i++) { for (j = 0; j < M; j++) STR printf("\n"); }}void printnn(float a[N][N]) { int i, j; for (i = 0; i < N; i++) { for (j = 0; j < N; j++) STR printf("\n"); }}void printmn(float a[m][N]) { int i, j; for (i = 0; i < M; i++) { for (j = 0; j < N; j++) STR printf("\n"); }}void printnm(float a[m][N]) { int i, j; for (j = 0; j < N; j++) { for (i = 0; i < M; i++) STR printf("\n"); }}int HouseM(int jj) { int i, j; double r; float x[m], u1[m][m]; for (i = 0; i < M; i++) x = 0.0; for (i = jj; i < M; i++) x = A[jj]; r = arraym(x, M); x[jj] -= r; r = arraym(x, M); if (r <= F) return 0; for (i = jj; i < M; i++) x /= r; Imm(u1); for (i = jj; i < M; i++) for (j = jj; j < M; j++) u1[j] += x * x[j] * ( - 2); mmmn2(u1, A); mmmm1(U, u1); return 1;}int HouseN(int jj) { int i, j, k; double r; float x[N], v1[N][N]; k = jj + 1; for (i = 0; i < N; i++) x = 0.0; for (i = k; i < N; i++) x = A[jj]; r = arraym(x, N); x[k] -= r; r = arraym(x, N); if (r <= F) return 0; for (i = k; i < N; i++) x /= r; Inn(v1); for (i = k; i < N; i++) for (j = k; j < N; j++) v1[j] += x * x[j] * ( - 2); mnnn1(A, v1); nnnn1(V, v1); return 1;}void Givens(void) { int k; float a, b, r, g1[m][m], g2[N][N]; a = A[p][p] * A[p][p] - miu(); b = A[p][p] * A[p][p + 1]; r = a * a + b * b; r = sqrt(r); Inn(g2); g2[p][p] = a / r; g2[p + 1][p + 1] = g2[p][p]; g2[p + 1][p] = b / r; g2[p][p + 1] = -g2[p + 1][p]; mnnn1(A, g2); nnnn1(V, g2); for (k = p; k < q; k++) { a = A[k][k]; b = A[k + 1][k]; r = a * a + b * b; r = sqrt(r); Imm(g1); g1[k][k] = a / r; g1[k + 1][k + 1] = g1[k][k]; g1[k][k + 1] = b / r; g1[k + 1][k] = -g1[k][k + 1]; mmmn2(g1, A); g1[k][k + 1] = g1[k + 1][k]; g1[k + 1][k] = -g1[k][k + 1]; mmmm1(U, g1); if (k < q - 1) { a = A[k][k + 1]; b = A[k][k + 2]; r = a * a + b * b; r = sqrt(r); Inn(g2); g2[k + 1][k + 1] = a / r; g2[k + 2][k + 2] = g2[k + 1][k + 1]; g2[k + 2][k + 1] = b / r; g2[k + 1][k + 2] = -g2[k + 2][k + 1]; mnnn1(A, g2); nnnn1(V, g2); } }}void paixu(void) { int i, j, k, h; float t; for (i = 0; i < M; i++) for (j = 0; j < N; j++) if (fabs(A[j]) <= F * 10) A[j] = 0.0; for (i = 0; i < M; i++) for (j = 0; j < M; j++) if (fabs(U[j]) <= F * 10) U[j] = 0.0; for (i = 0; i < N; i++) for (j = 0; j < N; j++) if (fabs(V[j]) <= F * 10) V[j] = 0.0; for (i = 0; i < N - 1; i++) { k = i; for (j = i + 1; j < N; j++) if (A[j][j] > A[k][k]) k = j; if (k > i) { t = A; A = A[k][k]; A[k][k] = t; for (h = 0; h < M; h++) { t = U[h]; U[h] = U[h][k]; U[h][k] = t; } for (h = 0; h < N; h++) { t = V[h]; V[h] = V[h][k]; V[h][k] = t; } } }}void main(void) { int i, j;#ifdef T for (i = 0; i < N; i++) for (j = 0; j < M; j++) { printf("A[%d][%d]=", i + 1, j + 1); scanf("%f", &A[j]); } printf("A=\n"); printnm(A);# else for (i = 0; i < M; i++) for (j = 0; j < N; j++) { printf("A[%d][%d]=", i + 1, j + 1); scanf("%f", &A[j]); } printf("A=\n"); printmn(A);#endif Imm(U); Inn(V); for (j = 0; j < N; j++) { HouseM(j); if (j < N - 2) HouseN(j); } p = 0; q = N - 1; while (q > 0) { line1: for (i = p; i < q; i++) if (fabs(A) <= F * (fabs(A) + fabs(A))) A = 0.0; for (i = q; i > 0; i--) if (fabs(A) > F) break; q = i; for (i = q; i > 0; i--) if (fabs(A) <= F) break; p = i; if (q == 0) break; for (i = p; i < q; i++) if (fabs(A) <= F * fanshu()) { A = 0.0; goto line1; } Givens(); } paixu();#ifdef T printf("U'AV=\n"); printnm(A); printf("U=\n"); printnn(V); printf("V=\n"); printmm(U);# else printf("U'AV=\n"); printmn(A); printf("U=\n"); printmm(U); printf("V=\n"); printnn(V);#endif} 0 0
