施密特正交化c语言,C语言实现矩阵的LU分解、施密特正交化、Givens分解、Householder分解...

施密特正交化算法（Schmidt orthogonalization algorithm）可以用 C 语言实现如下：

#include <stdio.h>
#include <math.h>

#define M 3
#define N 3

void schmidt(double A[M][N], double Q[M][N], double R[N][N]) {
    int i, j, k;
    double temp, norm;

    for (i = 0; i < N; i++) {
        for (j = 0; j < M; j++) {
            Q[j][i] = A[j][i];
        }
        for (j = 0; j < i; j++) {
            temp = 0.0;
            for (k = 0; k < M; k++) {
                temp += Q[k][i] * Q[k][j];
            }
            for (k = 0; k < M; k++) {
                Q[k][i] -= temp * Q[k][j];
            }
        }
        norm = 0.0;
        for (j = 0; j < M; j++) {
            norm += Q[j][i] * Q[j][i];
        }
        norm = sqrt(norm);
        for (j = 0; j < M; j++) {
            Q[j][i] /= norm;
        }
    }

    for (i = 0; i < N; i++) {
        for (j = i; j < N; j++) {
            R[i][j] = 0.0;
            for (k = 0; k < M; k++) {
                R[i][j] += Q[k][i] * A[k][j];
            }
        }
    }
}

LU 分解算法可以用 C 语言实现如下：

#include <stdio.h>

#define N 3

void lu(double A[N][N], double L[N][N], double U[N][N]) {
    int i, j, k;
    double sum;

    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++) {
            if (j < i) {
                sum = 0.0;
                for (k = 0; k < j; k++) {
                    sum += L[j][k] * U[k][i];
                }
                L[j][i] = (A[j][i] - sum) / U[j][j];
            } else if (j == i) {
                L[j][i] = 1.0;
                sum = 0.0;
                for (k = 0; k < j; k++) {
                    sum += L[j][k] * U[k][i];
                }
                U[j][i] = A[j][i] - sum;
            } else {
                U[j][i] = 0.0;
                sum = 0.0;
                for (k = 0; k < i; k++) {
                    sum += L[j][k] * U[k][i];
                }
                L[j][i] = (A[j][i] - sum) / U[i][i];
            }
        }
    }
}

Givens 分解算法可以用 C 语言实现如下：

#include <stdio.h>
#include <math.h>

#define N 3

void givens(double A[N][N], double Q[N][N], double R[N][N]) {
    int i, j, k;
    double c, s, temp;

    for (i = 0; i < N; i++) {
        for (j = i + 1; j < N; j++) {
            if (A[i][j] != 0.0) {
                if (fabs(A[i][i]) >= fabs(A[j][j])) {
                    temp = A[j][i] / A[i][i];
                    c = 1.0 / sqrt(1.0 + temp * temp);
                    s = temp * c;
                } else {
                    temp = A[i][j] / A[j][j];
                    s = 1.0 / sqrt(1.0 + temp * temp);
                    c = temp * s;
                }
                for (k = 0; k < N; k++) {
                    temp = A[i][k];
                    A[i][k] = c * temp + s * A[j][k];
                    A[j][k] = -s * temp + c * A[j][k];
                    temp = Q[i][k];
                    Q[i][k] = c * temp + s * Q[j][k];
                    Q[j][k] = -s * temp + c * Q[j][k];
                }
            }
        }
    }

    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++) {
            if (i <= j) {
                R[i][j] = A[i][j];
            } else {
                R[i][j] = 0.0;
            }
        }
    }
}

Householder 分解算法可以用 C 语言实现如下：

#include <stdio.h>
#include <math.h>

#define N 3

void householder(double A[N][N], double Q[N][N], double R[N][N]) {
    int i, j, k;
    double alpha, norm;
    double u[N], v[N];

    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++) {
            Q[i][j] = A[i][j];
        }
    }

    for (k = 0; k < N - 1; k++) {
        norm = 0.0;
        for (i = k; i < N; i++) {
            norm += Q[i][k] * Q[i][k];
        }
        alpha = sqrt(norm);
        if (Q[k][k] > 0.0) {
            alpha = -alpha;
        }
        norm = 0.0;
        for (i = k; i < N; i++) {
            u[i] = Q[i][k];
            norm += u[i] * u[i];
        }
        norm = sqrt(norm);
        for (i = k; i < N; i++) {
            u[i] /= norm;
        }
        u[k] += 1.0;
        norm = 0.0;
        for (i = k; i < N; i++) {
            v[i] = 0.0;
            for (j = k; j < N; j++) {
                v[i] += Q[i][j] * u[j];
            }
            norm += v[i] * v[i];
        }
        norm = sqrt(norm);
        for (i = k; i < N; i++) {
            v[i] /= norm;
        }
        for (i = k; i < N; i++) {
            for (j = k; j < N; j++) {
                Q[i][j] -= 2.0 * u[i] * v[j];
            }
        }
    }

    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++) {
            if (i <= j) {
                R[i][j] = 0.0;
                for (k = 0; k < N; k++) {
                    R[i][j] += Q[k][i] * A[k][j];
                }
            } else {
                R[i][j] = Q[i][j];
            }
        }
    }
}