Householder

来源：互联网发布：数控编程工资月薪过万编辑：程序博客网时间：2024/06/06 02:55
#include < stdio.h > #include < math.h > #define N 3 float Q[N][N],R[N][N],H[N][N];void printm(float a[N][N]) {    int i,    j;    printf("\n");    for (i = 0; i < N; i++) {        for (j = 0; j < N; j++) {            if (fabs(a[j]) <= 1e-5) a[j] = 0.0;            printf("%-12.5g", a[j]);        }        printf("\n");    }}void Tm(float a[N][N]) {    int i,    j;    float t;    for (i = 0; i < N; i++) for (j = i + 1; j < N; j++) {        t = a[j];        a[j] = a[j];        a[j] = t;    }}void copym(float a[N][N], float b[N][N]) {    int i,    j;    for (i = 0; i < N; i++) for (j = 0; j < N; j++) a[j] = b[j];}void mult2(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];    }    copym(b, t);}void I(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;    }}double arraym(float x[N]) {    int i;    double r;    r = 0.0;    for (i = 0; i < N; i++) r += (double)(x * x);    r = sqrt(r);    return r;}int Householder(int j) {    int i,    k;    double r;    float x[N];    for (i = 0; i < N; i++) x = 0.0;    for (i = j; i < N; i++) x = R[j];    r = arraym(x);    x[j] -= r;    r = arraym(x);    if (fabs(r) <= 1e-6) return 0;    for (i = j; i < N; i++) x /= r;    I(H);    for (i = j; i < N; i++) for (k = j; k < N; k++) H[k] += x * x[k] * ( - 2);    return 1;}void main() {    int i,    j,    t;    for (i = 0; i < N; i++) for (j = 0; j < N; j++) {        printf("A[%d][%d]=", i + 1, j + 1);        scanf("%f", &R[j]);    }    printf("\nA=\n");    printm(R);    I(Q);    for (j = 0; j < N - 1; j++) {        t = Householder(j);        if (t == 0) continue;        mult2(H, R);        mult2(H, Q);    }    Tm(Q);    printf("Q=\n");    printm(Q);    printf("R=\n");    printm(R);}
0 0