C++ source code to solve simultaneous equations with Gauss-Jorden(pivot) method.

gauss_jordern_pivot.cpp

/*********************************************
 * 連立方程式の解法 ( ガウス・ジョルダン法（ピボット選択法） )
 *********************************************/
#include <iostream>  // for cout
#include <math.h>    // for fabs()
#include <stdio.h>   // for printf()
#include <stdlib.h>  // for exit()

// 元の数定義
#define N 3  // 4

using namespace std;

/*
 * 計算クラス
 */
class Calc
{
    double a[N][N + 1];  // 係数用配列

    // 各種変数
    double p, d;         // ピボット係数、ピボット行ｘ係数
    double max, dummy;   // 最大絶対値、入れ替え時ダミー
    int i, j, k;         // LOOP インデックス
    int s;               // 入替行

    public:
        // 連立方程式を解く（ガウス・ジョルダン法（ピボット選択法））
        void calcGaussJordenPivot();
};

/*
 * 連立方程式を解く（ガウス・ジョルダン法（ピボット選択法））
 */
void Calc::calcGaussJordenPivot()
{
    // 係数
    static double a[N][N + 1] = {
        { 2.0, -3.0,  1.0,  5.0},
        { 1.0,  1.0, -1.0,  2.0},
        { 3.0,  5.0, -7.0,  0.0}
        //{ 1.0, -2.0,  3.0, -4.0,  5.0},
        //{-2.0,  5.0,  8.0, -3.0,  9.0},
        //{ 5.0,  4.0,  7.0,  1.0, -1.0},
        //{ 9.0,  7.0,  3.0,  5.0,  4.0}
    };

    // 元の連立方程式をコンソール出力
    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++)
            printf("%+fx%d ", a[i][j], j + 1);
        printf("= %+f\n", a[i][N]);
    }

    for (k = 0; k < N; k++) {
        // 行入れ替え
        max = 0; s = k;
        for (j = k; j < N; j++) {
            if (fabs(a[j][k]) > max) {
                max = fabs(a[j][k]);
                s = j;
            }
        }
        if (max == 0) {
            printf("解けない！");
            exit(1);
        }
        for (j = 0; j <= N; j++) {
            dummy   = a[k][j];
            a[k][j] = a[s][j];
            a[s][j] = dummy;
        }

        // ピボット係数
        p = a[k][k];

        // ピボット行を p で除算
        for (j = k; j < N + 1; j++)
            a[k][j] /= p;

        // ピボット列の掃き出し
        for (i = 0; i < N; i++) {
            if (i != k) {
                d = a[i][k];
                for (j = k; j < N + 1; j++)
                    a[i][j] -= d * a[k][j];
            }
        }
    }

    // 結果出力
    for (k = 0; k < N; k++)
        printf("x%d = %f\n", k + 1, a[k][N]);
}

/*
 * メイン処理
 */
int main()
{
    try
    {
        // 計算クラスインスタンス化
        Calc objCalc;

        // 連立方程式を解く（ガウス・ジョルダン法（ピボット選択法））
        objCalc.calcGaussJordenPivot();
    }
    catch (...) {
        cout << "例外発生！" << endl;
        return -1;
    }

    // 正常終了
    return 0;
}