求解线性方程组--全主元消去法

求解线性方程组--全主元消去法。

源代码用C++实现~~这里采用的是数值分析里的算法。

代码在VC6.0下编译通过，经测试没大问题。

/* 求解线性方程组--全主元消去法 */

#include <iostream>
using namespace std;

/* 二维数组动态分配模板 */
template <class T>
T** Allocation2D(int m, int n)
{
    T **a;
    a = new T*[m];
    for (int i = 0; i < m; i++) {
        a[i] = new T[n];
    }
    return a;
}

/* 一维数组动态分配模板 */
template <class T>
T* Allocation1D(int n)
{
    T *a;
    a = new T[n];
    return a;
}

int main()
{
    /* 循环变量 */
    int i, j, k;

    /* 系数矩阵的行数 */
    int n;

    /* 增广矩阵 */
    float** a;

    /* 动态生成增广矩阵 */
    cout << "输入系数矩阵的N值，N：";
    cin >> n;
    a = Allocation2D<float>(n, n + 1);

    /* 输入增广矩阵的各值 */
    cout << endl << "输入增广矩阵的各值：/n";
    for (i = 0; i < n; i++) {
        for (j = 0; j < n + 1; j++) {
            cin >> a[i][j];
        }
    }

    /* e为记录列交换过程的数组，用在最后理顺解向量 */
    int* e;
    e = Allocation1D<int>(n);
    for (i = 0; i < n; i++) {
        e[i] = i;
    }

    float temp;
    int row, col;
    for (k = 0; k < n - 1; k++) {
        temp = 0;
        row = 0;
        col = 0;

        /* 选主元 */
        for (i = k; i < n; i++) {
            for (j = k; j < n; j++) {
                if (a[i][j] > temp) {
                    temp = a[i][j];
                    row = i;
                    col = j;
                }
            }
        }

        if (0 == temp) {
            cout << "系数矩阵为奇异阵！/n";
            return 0;
        }

        /* 行交换 */
        if (row != k) {
            for (j = k; j < n + 1; j++) {
                temp = a[k][j];
                a[k][j] = a[row][j];
                a[row][j] = temp;
            }
        }
        /* 列交换 */
        e[k] = k;
        if (col != k) {
            e[k] = col;
            for (i = 0; i < n; i++) {
                temp = a[i][k];
                a[i][k] = a[i][col];
                a[i][col] = temp;
            }
        }
        for (i = k + 1; i < n; i++) {
            for (j = k + 1; j < n + 1; j++) {
                a[i][j] = a[i][j] - a[i][k] * a[k][j] / a[k][k];
            }
        }
    }

    /* 回代过程 */
    a[n - 1][n] = a[n - 1][n] / a[n - 1][n - 1];
    for (k = n - 2; k >= 0; k--) {
        temp = 0;
        for (j = k + 1; j < n; j++) {
            temp = temp + a[k][j] * a[j][n];
        }
        a[k][n] = (a[k][n] - temp) / a[k][k];
    }

    /* 理顺解向量 */
    for (k = n - 2; k >= 0; k--) {
        if (e[k] != k) {
            temp = a[e[k]][n];
            a[e[k]][n] = a[k][n];
            a[k][n] = temp;
        }
    }

    /* 输出过程 */
    cout << "解向量为：/n";
    for (i = 0; i < n; i++) {
        cout << "x" << i << ": " << a[i][n] << endl;
    }
    return 0;
}