home > さまざまな言語で数値計算 > 連立一次方程式 >

さまざまな言語で数値計算

Only Do What Only You Can Do
home | index | 前へ | 次へ

ガウスの消去法

ガウスの消去法を使って, 連立一次方程式の解を求める .

例として, 連立一次方程式

を考える.

1. 前進消去

まず(1)式を 9 で割って(1)'とする.

(2)式から(1)'式に 2 を掛けたものを引くと $x$ の項が消える.

(3)式に(1)'式を足すと $x$ の項が消え, (4)式から(1)'式を引くと $x$ の項が消える.
同様にして, 2番目の式を利用して, 3番目以降の式の $y$ を消す.
3番目の式を利用して, 4番目の式の $z$ を消す.
このようにして、上三角型方程式を得る.

これを、前進消去という.

2. 後退代入

(4)式から $u$ の値が求まるので, (3)式に代入すると $z$ の値が求まる.
$u, z$ の値を (2)式に代入すると $y$ の値が求まる.
$u, z, y$ の値を (1)式に代入すると $x$ の値が求まる.
これを、後退代入という.

VBScript

Option Explicit

Private Const N = 4

Private a: a = Array(Array(-1,-2,7,-2),Array(1,-1,-2,6),Array(9,2,1,1),Array(2,8,-2,1))
Private b: b = Array(8,17,20,16)

'ピボット選択
pivoting a,b

WScript.StdOut.WriteLine "pivoting"
WScript.StdOut.WriteLine "A"
disp_matrix a
WScript.StdOut.WriteLine "B"
disp_vector b
WScript.StdOut.WriteLine

'前進消去
forward_elimination a,b

WScript.StdOut.WriteLine "forward elimination"
WScript.StdOut.WriteLine "A"
disp_matrix a
WScript.StdOut.WriteLine "B"
disp_vector b
WScript.StdOut.WriteLine

'後退代入
backward_substitution a,b

WScript.StdOut.WriteLine "X"
disp_vector b

'前進消去
Private Sub forward_elimination(ByRef a, ByRef b)
    Dim pivot, row, col
    For pivot = 0 To (N - 2)
        For row = (pivot + 1) To (N - 1)
            Dim s: s = a(row)(pivot) / a(pivot)(pivot)
            For col = pivot To (N - 1)
                a(row)(col) = a(row)(col) - a(pivot)(col) * s
            Next
            b(row) = b(row) - b(pivot) * s
        Next
    Next
End Sub
'後退代入
Private Sub backward_substitution(ByVal a, ByRef b)
    Dim row, col
    For row = (N - 1) To 0 Step -1
        For col = (N - 1) To (row + 1) Step - 1
            b(row) = b(row) - a(row)(col) * b(col)
        Next
        b(row) = b(row) / a(row)(row)
    Next
End Sub
'ピボット選択
Private Sub pivoting(ByRef a, ByRef b)
    Dim pivot, row, col
    For pivot = 0 To (N - 1)
        '各列で 一番値が大きい行を 探す
        Dim max_row: max_row =   pivot
        Dim max_val: max_val =   0
        For row = pivot To (N - 1)
            If Abs(a(row)(pivot)) > max_val Then
                '一番値が大きい行
                max_val =   Abs(a(row)(pivot))
                max_row =   row
            End If
        Next

        '一番値が大きい行と入れ替え
        If max_row <> pivot Then
            Dim tmp
            For col = 0 To (N - 1)
                tmp             =   a(max_row)(col)
                a(max_row)(col) =   a(pivot)(col)
                a(pivot)(col)   =   tmp
            Next
            tmp         =   b(max_row)
            b(max_row)  =   b(pivot)
            b(pivot)    =   tmp
        End If
    Next
End Sub
'１次元配列を表示
Private Sub disp_vector(ByVal row())
    Dim col
    For Each col In row
        WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab
    Next
    WScript.StdOut.WriteLine
End Sub
'２次元配列を表示
Private Sub disp_matrix(ByVal matrix)
    Dim row, col
    For Each row In matrix
        For Each col In row
            WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab
        Next
        WScript.StdOut.WriteLine
    Next
End Sub

Z:\>cscript //nologo Z:\1003.vbs
pivoting
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  2.0000000000    8.0000000000   -2.0000000000    1.0000000000
 -1.0000000000   -2.0000000000    7.0000000000   -2.0000000000
  1.0000000000   -1.0000000000   -2.0000000000    6.0000000000
B
 20.0000000000   16.0000000000    8.0000000000   17.0000000000

forward elimination
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.0000000000    7.5555555556   -2.2222222222    0.7777777778
  0.0000000000    0.0000000000    6.5882352941   -1.7058823529
  0.0000000000    0.0000000000    0.0000000000    5.3750000000
B
 20.0000000000   11.5555555556   12.9411764706   21.5000000000

X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

JScript

var N = 4

var a = [[-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1]]
var b = [8,17,20,16]

// ピボット選択
pivoting(a,b)

WScript.StdOut.WriteLine("pivoting")
WScript.StdOut.WriteLine("A")
disp_matrix(a)
WScript.StdOut.WriteLine("B")
disp_vector(b)
WScript.StdOut.WriteLine()

// 前進消去
forward_elimination(a,b)

WScript.StdOut.WriteLine("forward elimination")
WScript.StdOut.WriteLine("A")
disp_matrix(a)
WScript.StdOut.WriteLine("B")
disp_vector(b)
WScript.StdOut.WriteLine()

// 後退代入
backward_substitution(a,b)

WScript.StdOut.WriteLine("X")
disp_vector(b)

// 前進消去
function forward_elimination(a, b)
{
    for (pivot = 0; pivot < N - 1; pivot++)
    {
        for (row = pivot + 1; row < N; row++)
        {
            var s = a[row][pivot] / a[pivot][pivot]
            for (col = pivot; col < N; col++)
                a[row][col] -= a[pivot][col]    * s
            b[row]          -= b[pivot]         * s
        }
    }
}
// 後退代入
function backward_substitution(a, b)
{
    for (row = N - 1; row >= 0; row--)
    {
        for (col = N - 1; col > row; col--)
            b[row] -= a[row][col] * b[col]
        b[row] /= a[row][row]
    }
}
// ピボット選択
function pivoting(a, b)
{
    for (pivot = 0; pivot < N; pivot++)
    {
        // 各列で 一番値が大きい行を 探す
        var max_row =   pivot
        var max_val =   0
        for (row = pivot; row < N; row++)
        {
            if (Math.abs(a[row][pivot]) > max_val)
            {
                // 一番値が大きい行
                max_val =   Math.abs(a[row][pivot])
                max_row =   row
            }
        }

        // 一番値が大きい行と入れ替え
        if (max_row != pivot)
        {
            var tmp
            for (col = 0; col < N; col++)
            {
                tmp             =   a[max_row][col]
                a[max_row][col] =   a[pivot][col]
                a[pivot][col]   =   tmp
            }
            tmp         =   b[max_row]
            b[max_row]  =   b[pivot]
            b[pivot]    =   tmp
        }
    }
}
// １次元配列を表示
function disp_vector(row)
{
    for (var i = 0; i < N; i++)
        WScript.StdOut.Write(("              "    + row[i].toFixed(10)).slice(-14) + "\t")
    WScript.StdOut.WriteLine()
}
// ２次元配列を表示
function disp_matrix(matrix)
{
    for (i = 0; i < N; i++)
    {
        for (j = 0; j < N; j++)
            WScript.StdOut.Write(("              "    + matrix[i][j].toFixed(10)).slice(-14) + "\t")
        WScript.StdOut.WriteLine()
    }
}

Z:\>cscript //nologo Z:\1003.js
pivoting
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  2.0000000000    8.0000000000   -2.0000000000    1.0000000000
 -1.0000000000   -2.0000000000    7.0000000000   -2.0000000000
  1.0000000000   -1.0000000000   -2.0000000000    6.0000000000
B
 20.0000000000   16.0000000000    8.0000000000   17.0000000000

forward elimination
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.0000000000    7.5555555556   -2.2222222222    0.7777777778
  0.0000000000    0.0000000000    6.5882352941   -1.7058823529
  0.0000000000    0.0000000000    0.0000000000    5.3750000000
B
 20.0000000000   11.5555555556   12.9411764706   21.5000000000

X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

PowerShell

set-variable -name N -value 3 -option constant

# 前進消去
function forward_elimination($a, $b)
{
    foreach ($pivot in 0..($N - 1))
    {
        foreach ($row in ($pivot + 1)..$N)
        {
            $s = $a[$row][$pivot] / $a[$pivot][$pivot]
            foreach ($col in $pivot..$N)
            {
                $a[$row][$col] -= $a[$pivot][$col]    * $s
            }
            $b[$row]           -= $b[$pivot]          * $s
        }
    }
}
# 後退代入
function backward_substitution($a, $b)
{
    foreach ($row in $N..0)
    {
        if (($row + 1) -le $N)
        {
            foreach ($col in $N..($row + 1))
            {
                $b[$row] -= $a[$row][$col] * $b[$col]
            }
        }
        $b[$row] /= $a[$row][$row]
    }
}
# ピボット選択
function pivoting($a, $b)
{
    foreach ($pivot in 0..$N)
    {
        # 各列で 一番値が大きい行を 探す
        $max_row =   $pivot
        $max_val =   0
        foreach ($row in $pivot..$N)
        {
            if ([Math]::Abs($a[$row][$pivot]) -gt $max_val)
            {
                # 一番値が大きい行
                $max_val =   [Math]::Abs($a[$row][$pivot])
                $max_row =   $row
            }
        }

        # 一番値が大きい行と入れ替え
        if ($max_row -ne $pivot)
        {
            foreach ($col in 0..$N)
            {
                $tmp               =   $a[$max_row][$col]
                $a[$max_row][$col] =   $a[$pivot][$col]
                $a[$pivot][$col]   =   $tmp
            }
            $tmp         =   $b[$max_row]
            $b[$max_row] =   $b[$pivot]
            $b[$pivot]   =   $tmp
        }
    }
}
# １次元配列を表示
function disp_vector($row)
{
    foreach ($col in $row)
    {
        Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline
    }
    Write-Host
}
# ２次元配列を表示
function disp_matrix($matrix)
{
    foreach ($row in $matrix)
    {
        foreach ($col in $row)
        {
            Write-Host ([String]::Format("{0,14:F10}", $col)) -nonewline
        }
        Write-Host
    }
}

$a = (-1.0,-2.0,7.0,-2.0), (1.0,-1.0,-2.0,6.0), (9.0,2.0,1.0,1.0), (2.0,8.0,-2.0,1.0)
$b = 8.0, 17.0, 20.0, 16.0

# ピボット選択
pivoting $a $b

Write-Host "pivoting"
Write-Host "A"
disp_matrix $a
Write-Host "B"
disp_vector $b
Write-Host ""

# 前進消去
forward_elimination $a $b

Write-Host "forward elimination"
Write-Host "A"
disp_matrix $a
Write-Host "B"
disp_vector $b
Write-Host ""

# 後退代入
backward_substitution $a $b

Write-Host "X"
disp_vector $b

Z:\>powershell -file Z:\1003.ps1
pivoting
A
  9.0000000000  2.0000000000  1.0000000000  1.0000000000
  2.0000000000  8.0000000000 -2.0000000000  1.0000000000
 -1.0000000000 -2.0000000000  7.0000000000 -2.0000000000
  1.0000000000 -1.0000000000 -2.0000000000  6.0000000000
B
 20.0000000000 16.0000000000  8.0000000000 17.0000000000

forward elimination
A
  9.0000000000  2.0000000000  1.0000000000  1.0000000000
  0.0000000000  7.5555555556 -2.2222222222  0.7777777778
  0.0000000000  0.0000000000  6.5882352941 -1.7058823529
  0.0000000000  0.0000000000  0.0000000000  5.3750000000
B
 20.0000000000 11.5555555556 12.9411764706 21.5000000000

X
  1.0000000000  2.0000000000  3.0000000000  4.0000000000
TOP

Perl

use constant N => 3;

my @a = ([-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1]);
my @b = (8,17,20,16);

# ピボット選択
pivoting(\@a,\@b);
print "pivoting\n";
print "A\n";
disp_matrix(\@a);
print "B\n";
disp_vector(\@b);
print "\n";

# 前進消去
forward_elimination(\@a,\@b);
print "forward elimination\n";
print "A\n";
disp_matrix(\@a);
print "B\n";
disp_vector(\@b);
print "\n";

# 後退代入
backward_substitution(\@a,\@b);
print "X\n";
disp_vector(\@b);

# 前進消去
sub forward_elimination
{
    my ($a, $b) = @_;

    for $pivot (0..(N - 1))
    {
        for $row (($pivot + 1)..N)
        {
            my $s = $$a[$row][$pivot] / $$a[$pivot][$pivot];
            for $col ($pivot..N)
            {
                $$a[$row][$col] -= $$a[$pivot][$col]    * $s;
            }
            $$b[$row]           -= $$b[$pivot]          * $s;
        }
    }
}
# 後退代入
sub backward_substitution
{
    my ($a, $b) = @_;

    for ($row = N; $row >= 0; $row--)
    {
        for ($col = N; $col > $row; $col--)
        {
            $$b[$row] -= $$a[$row][$col] * $$b[$col];
        }
        $$b[$row] /= $$a[$row][$row];
    }
}
# ピボット選択
sub pivoting
{
    my ($a, $b) = @_;

    for $pivot (0..N)
    {
        # 各列で 一番値が大きい行を 探す
        my $max_row =   $pivot;
        my $max_val =   0;
        for $row ($pivot..N)
        {
            if (abs($$a[$row][$pivot]) > $max_val)
            {
                # 一番値が大きい行
                $max_val =   abs($$a[$row][$pivot]);
                $max_row =   $row;
            }
        }

        # 一番値が大きい行と入れ替え
        if ($max_row != $pivot)
        {
            my $tmp;
            for $col (0..N)
            {
                $tmp                =   $$a[$max_row][$col];
                $$a[$max_row][$col] =   $$a[$pivot][$col];
                $$a[$pivot][$col]   =   $tmp;
            }
            $tmp          =   $$b[$max_row];
            $$b[$max_row] =   $$b[$pivot];
            $$b[$pivot]   =   $tmp;
        }
    }
}
# ２次元配列を表示
sub disp_matrix
{
    my ($matrix) = @_;
    foreach $row (@$matrix)
    {
        foreach $col (@$row)
        {
            printf("%14.10f\t", $col);
        }
        print "\n";
    }
}
# １次元配列を表示
sub disp_vector
{
    my ($row) = @_;
    foreach $col (@$row)
    {
        printf("%14.10f\t", $col);
    }
    print "\n";
}

Z:\>perl Z:\1003.pl
pivoting
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  2.0000000000    8.0000000000   -2.0000000000    1.0000000000
 -1.0000000000   -2.0000000000    7.0000000000   -2.0000000000
  1.0000000000   -1.0000000000   -2.0000000000    6.0000000000
B
 20.0000000000   16.0000000000    8.0000000000   17.0000000000

forward elimination
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.0000000000    7.5555555556   -2.2222222222    0.7777777778
  0.0000000000    0.0000000000    6.5882352941   -1.7058823529
  0.0000000000    0.0000000000    0.0000000000    5.3750000000
B
 20.0000000000   11.5555555556   12.9411764706   21.5000000000

X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

PHP

<?php
define("N", 3);

$a = [[-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1]];
$b = [8,17,20,16];

# ピボット選択
pivoting($a,$b);
print "pivoting\n";
print "A\n";
disp_matrix($a);
print "B\n";
disp_vector($b);
print "\n";

# 前進消去
forward_elimination($a,$b);
print "forward elimination\n";
print "A\n";
disp_matrix($a);
print "B\n";
disp_vector($b);
print "\n";

# 後退代入
backward_substitution($a,$b);
print "X\n";
disp_vector($b);

# 前進消去
function forward_elimination(&$a, &$b)
{
    foreach (range(0, N - 1) as $pivot)
    {
        foreach (range($pivot + 1, N) as $row)
        {
            $s = $a[$row][$pivot] / $a[$pivot][$pivot];
            foreach (range($pivot, N) as $col)
                $a[$row][$col] -= $a[$pivot][$col]    * $s;
            $b[$row]           -= $b[$pivot]          * $s;
        }
    }
}
# 後退代入
function backward_substitution(&$a, &$b)
{
    foreach (range(N, 0) as $row)
    {
        if (($row + 1) <= N)
            foreach (range(N, $row + 1) as $col)
                $b[$row] -= $a[$row][$col] * $b[$col];
        $b[$row] /= $a[$row][$row];
    }
}
# ピボット選択
function pivoting(&$a, &$b)
{
    foreach (range(0, N) as $pivot)
    {
        # 各列で 一番値が大きい行を 探す
        $max_row =   $pivot;
        $max_val =   0;
        foreach (range($pivot, N) as $row)
        {
            if (abs($a[$row][$pivot]) > $max_val)
            {
                # 一番値が大きい行
                $max_val =   abs($a[$row][$pivot]);
                $max_row =   $row;
            }
        }

        # 一番値が大きい行と入れ替え
        if ($max_row != $pivot)
        {
            foreach (range(0, N) as $col)
            {
                $tmp               =   $a[$max_row][$col];
                $a[$max_row][$col] =   $a[$pivot][$col];
                $a[$pivot][$col]   =   $tmp;
            }
            $tmp         =   $b[$max_row];
            $b[$max_row] =   $b[$pivot];
            $b[$pivot]   =   $tmp;
        }
    }
}
# ２次元配列を表示
function disp_matrix($matrix)
{
    foreach ($matrix as $row)
    {
        foreach ($row as $col)
            printf("%14.10f\t", $col);
        print "\n";
    }
}
# １次元配列を表示
function disp_vector($row)
{
    foreach ($row as $col)
        printf("%14.10f\t", $col);
    print "\n";
}?>

Z:\>php Z:\1003.php
pivoting
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  2.0000000000    8.0000000000   -2.0000000000    1.0000000000
 -1.0000000000   -2.0000000000    7.0000000000   -2.0000000000
  1.0000000000   -1.0000000000   -2.0000000000    6.0000000000
B
 20.0000000000   16.0000000000    8.0000000000   17.0000000000

forward elimination
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.0000000000    7.5555555556   -2.2222222222    0.7777777778
  0.0000000000    0.0000000000    6.5882352941   -1.7058823529
  0.0000000000    0.0000000000    0.0000000000    5.3750000000
B
 20.0000000000   11.5555555556   12.9411764706   21.5000000000

X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

Python

# coding: Shift_JIS

N = 4

# 前進消去
def forward_elimination(a, b):
    for pivot in range(0, N - 1, 1):
        for row in range(pivot + 1, N, 1):
            s = a[row][pivot] / a[pivot][pivot]
            for col in range(pivot, N, 1):
                a[row][col] -= a[pivot][col]    * s
            b[row]          -= b[pivot]         * s

# 後退代入
def backward_substitution(a, b):
    for row in range(N - 1, -1, -1):
        for col in range(N - 1, row, -1):
            b[row] -= a[row][col] * b[col]
        b[row] /= a[row][row]

# ピボット選択
def pivoting(a, b):
    for pivot in range(0, N, 1):
        # 各列で 一番値が大きい行を 探す
        max_row =   pivot
        max_val =   0.0
        for row in range(pivot, N, 1):
            if (abs(a[row][pivot]) > max_val):
                # 一番値が大きい行
                max_val =   abs(a[row][pivot])
                max_row =   row

        # 一番値が大きい行と入れ替え
        if (max_row != pivot):
            tmp = 0.0
            for col in range(0, N, 1):
                tmp             =   a[max_row][col]
                a[max_row][col] =   a[pivot][col]
                a[pivot][col]   =   tmp
            tmp        =   b[max_row]
            b[max_row] =   b[pivot]
            b[pivot]   =   tmp

# １次元配列を表示
def disp_matrix(matrix):
    for row in matrix:
        for col in row:
            print "%14.10f\t" % col,
        print ""

# ２次元配列を表示
def disp_vector(row):
    for col in row:
        print "%14.10f\t" % col,
    print ""

a = [[-1.0, -2.0,  7.0, -2.0], [1.0, -1.0, -2.0, 6.0], [ 9.0,  2.0, 1.0,  1.0], [2.0, 8.0, -2.0, 1.0], ]
b =  [ 8.0, 17.0, 20.0, 16.0]
c =  [ 0.0,  0.0,  0.0,  0.0]

# ピボット選択
pivoting(a,b)

print "pivoting"
print "A"
disp_matrix(a)
print "B"
disp_vector(b)
print ""

# 前進消去
forward_elimination(a,b)

print "forward elimination"
print "A"
disp_matrix(a)
print "B"
disp_vector(b)
print ""

# 後退代入
backward_substitution(a,b)

print "X"
disp_vector(b)

Z:\>python Z:\1003.py
pivoting
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  2.0000000000    8.0000000000   -2.0000000000    1.0000000000
 -1.0000000000   -2.0000000000    7.0000000000   -2.0000000000
  1.0000000000   -1.0000000000   -2.0000000000    6.0000000000
B
 20.0000000000   16.0000000000    8.0000000000   17.0000000000

forward elimination
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.0000000000    7.5555555556   -2.2222222222    0.7777777778
  0.0000000000    0.0000000000    6.5882352941   -1.7058823529
  0.0000000000    0.0000000000    0.0000000000    5.3750000000
B
 20.0000000000   11.5555555556   12.9411764706   21.5000000000

X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

Ruby

# coding: Shift_JIS

N = 3

# １次元配列を表示
def disp_matrix(matrix)
    matrix.each do |row|
        row.each do |col|
            printf("%14.10f\t", col)
        end
        puts ""
    end
end
# ２次元配列を表示
def disp_vector(row)
    row.each do |col|
        printf("%14.10f\t", col)
    end
    puts ""
end
# 前進消去
def forward_elimination(a, b)
    0.upto(N - 1) do |pivot|
        (pivot + 1).upto(N) do |row|
            s = a[row][pivot] / a[pivot][pivot]
            pivot.upto(N) do |col|
                a[row][col] -= a[pivot][col] * s
            end
            b[row]          -= b[pivot]      * s
        end
    end
end
# 後退代入
def backward_substitution(a, b)
    N.downto(0) do |row|
        N.downto(row + 1) do |col|
            b[row] -= a[row][col] * b[col]
        end
        b[row] /= a[row][row]
    end
end
# ピボット選択
def pivoting(a, b)
    (0..N).each do |pivot|
        # 各列で 一番値が大きい行を 探す
        max_row =   pivot
        max_val =   0
        (pivot..N).each do |row|
            if ((a[row][pivot]).abs > max_val)
                # 一番値が大きい行
                max_val =   (a[row][pivot]).abs
                max_row =   row
            end
        end

        # 一番値が大きい行と入れ替え
        if (max_row != pivot)
            tmp = 0
            (0..N).each do |col|
                tmp             =   a[max_row][col]
                a[max_row][col] =   a[pivot][col]
                a[pivot][col]   =   tmp
            end
            tmp        =   b[max_row]
            b[max_row] =   b[pivot]
            b[pivot]   =   tmp
        end
    end
end


a = [[-1.0,-2.0,7.0,-2.0],[1.0,-1.0,-2.0,6.0],[9.0,2.0,1.0,1.0],[2.0,8.0,-2.0,1.0]]
b =  [8.0,17.0,20.0,16.0]

# ピボット選択
pivoting(a,b)

puts "pivoting"
puts "A"
disp_matrix(a)
puts "B"
disp_vector(b)
puts ""

# 前進消去
forward_elimination(a,b)

puts "forward elimination"
puts "A"
disp_matrix(a)
puts "B"
disp_vector(b)
puts ""

# 後退代入
backward_substitution(a,b)

puts "X"
disp_vector(b)

Z:\>ruby Z:\1003.rb
pivoting
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  2.0000000000    8.0000000000   -2.0000000000    1.0000000000
 -1.0000000000   -2.0000000000    7.0000000000   -2.0000000000
  1.0000000000   -1.0000000000   -2.0000000000    6.0000000000
B
 20.0000000000   16.0000000000    8.0000000000   17.0000000000

forward elimination
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.0000000000    7.5555555556   -2.2222222222    0.7777777778
  0.0000000000    0.0000000000    6.5882352941   -1.7058823529
  0.0000000000    0.0000000000    0.0000000000    5.3750000000
B
 20.0000000000   11.5555555556   12.9411764706   21.5000000000

X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

Groovy

N = 3

def a = [[ 9.0,  2.0, 1.0,  1.0], [2.0, 8.0, -2.0, 1.0], [-1.0, -2.0, 7.0, -2.0], [1.0, -1.0, -2.0, 6.0]]  as double[][]
def b =  [20.0, 16.0, 8.0, 17.0] as double[]
def c =  [ 0.0,  0.0, 0.0,  0.0] as double[]

// ピボット選択
pivoting(a,b)

println("pivoting")
println("A")
disp_matrix(a)
println("B")
disp_vector(b)
println()

// 前進消去
forward_elimination(a,b)

println("forward elimination")
println("A")
disp_matrix(a)
println("B")
disp_vector(b)
println()

// 後退代入
backward_substitution(a,b)

println("X")
disp_vector(b)

// 前進消去
def forward_elimination(a, b) {
    for (pivot in 0..(N - 1)) {
        for (row in (pivot + 1)..N) {
            s = a[row][pivot] / a[pivot][pivot]
            for (col in pivot..N)
                a[row][col] -= a[pivot][col] * s
            b[row]          -= b[pivot]      * s
        }
    }
}
// 後退代入
def backward_substitution(a, b) {
    for (row = N; row >= 0; row--) {
        for (col = N; col > row; col--)
            b[row] -= a[row][col] * b[col]
        b[row] /= a[row][row]
    }
}
// ピボット選択
def pivoting(a, b) {
    for (pivot in 0..N) {
        // 各列で 一番値が大きい行を 探す
        max_row =   pivot
        max_val =   0
        for (row in pivot..N) {
            if (Math.abs(a[row][pivot]) > max_val) {
                // 一番値が大きい行
                max_val =   Math.abs(a[row][pivot])
                max_row =   row
            }
        }

        // 一番値が大きい行と入れ替え
        if (max_row != pivot) {
            tmp = 0
            for (col in 0..N) {
                tmp             =   a[max_row][col]
                a[max_row][col] =   a[pivot][col]
                a[pivot][col]   =   tmp
            }
            tmp        =   b[max_row]
            b[max_row] =   b[pivot]
            b[pivot]   =   tmp
        }
    }
}
// １次元配列を表示
def disp_matrix(matrix) {
    for (row in matrix) {
        for (col in row) {
            printf ("%14.10f\t" , col)
        }
        println()
    }
}
// ２次元配列を表示
def disp_vector(row) {
    for (col in row) {
        printf ("%14.10f\t" , col)
    }
    println()
}

Z:\>groovy Groovy1003.groovy
pivoting
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  2.0000000000    8.0000000000   -2.0000000000    1.0000000000
 -1.0000000000   -2.0000000000    7.0000000000   -2.0000000000
  1.0000000000   -1.0000000000   -2.0000000000    6.0000000000
B
 20.0000000000   16.0000000000    8.0000000000   17.0000000000

forward elimination
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.0000000000    7.5555555556   -2.2222222222    0.7777777778
  0.0000000000    0.0000000000    6.5882352941   -1.7058823529
  0.0000000000    0.0000000000    0.0000000000    5.3750000000
B
 20.0000000000   11.5555555556   12.9411764706   21.5000000000

X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

Pascal

program Pas1003(arg);
{$MODE delphi}

uses
    SysUtils, Math;

const
    N = 3;

type
    TwoDimArray = array [0..N,0..N] of Double;

// １次元配列を表示
procedure disp_vector(row:array of Double);
var
    i:Integer;
begin
    for i := Low(row) to High(row) do
        write(format('%14.10f'#9, [row[i]]));
    writeln();
end;
// ２次元配列を表示
procedure disp_matrix(matrix:TwoDimArray);
var
    i,j:Integer;
begin
    for i := Low(matrix) to High(matrix) do
    begin
        for j := Low(matrix) to High(matrix) do
            write(format('%14.10f'#9, [matrix[i,j]]));
        writeln();
    end;
end;
// 前進消去
procedure forward_elimination(var a:TwoDimArray; var b:array of Double);
var
    s:Double;
    pivot, row, col:Integer;
begin
    for pivot := Low(a) to High(a) - 1 do
    begin
        for row := (pivot + 1) to High(a) do
        begin
            s := a[row,pivot] / a[pivot,pivot];
            for col := pivot to High(a) do
                a[row,col] := a[row,col] - a[pivot,col] * s;
            b[row]         := b[row]     - b[pivot]     * s;
        end;
    end;
end;
// 後退代入
procedure backward_substitution(a:TwoDimArray; var b:array of Double);
var
    row, col:Integer;
begin
    for row := High(a) downto Low(a) do
    begin
        for col := High(a) downto row + 1 do
            b[row] := b[row] - a[row,col] * b[col];
        b[row] := b[row] / a[row,row];
    end;
end;
// ピボット選択
procedure pivoting(var a:TwoDimArray; var b:array of Double);
var
    max_val, tmp:Double;
    pivot, row, col, max_row:Integer;
begin
    for pivot := Low(a) to High(a) do
    begin
        // 各列で 一番値が大きい行を 探す
        max_row :=   pivot;
        max_val :=   0;
        for row := pivot to High(a) do
        begin
            if Abs(a[row,pivot]) > max_val then
            begin
                // 一番値が大きい行
                max_val :=   Abs(a[row,pivot]);
                max_row :=   row;
            end;
        end;

        // 一番値が大きい行と入れ替え
        if max_row <> pivot then
        begin
            tmp := 0;
            for col := Low(a) to High(a) do
            begin
                tmp            :=   a[max_row,col];
                a[max_row,col] :=   a[pivot,col];
                a[pivot,col]   :=   tmp;
            end;
            tmp        :=   b[max_row];
            b[max_row] :=   b[pivot];
            b[pivot]   :=   tmp;
        end;
    end;
end;

var
    a :TwoDimArray = ((-1.0, -2.0, 7.0, -2.0), (1.0, -1.0, -2.0, 6.0), ( 9.0,  2.0, 1.0,  1.0), (2.0, 8.0, -2.0, 1.0));
    b :array [0..N] of Double = (8.0, 17.0, 20.0, 16.0);
begin
    // ピボット選択
    pivoting(a,b);

    writeln('pivoting');
    writeln('A');
    disp_matrix(a);
    writeln('B');
    disp_vector(b);
    writeln();

    // 前進消去
    forward_elimination(a,b);

    writeln('forward elimination');
    writeln('A');
    disp_matrix(a);
    writeln('B');
    disp_vector(b);
    writeln();

    // 後退代入
    backward_substitution(a,b);

    writeln('X');
    disp_vector(b);
end.

Z:\>fpc -v0 -l- Pas1003.pp

Z:\>Pas1003
pivoting
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  2.0000000000    8.0000000000   -2.0000000000    1.0000000000
 -1.0000000000   -2.0000000000    7.0000000000   -2.0000000000
  1.0000000000   -1.0000000000   -2.0000000000    6.0000000000
B
 20.0000000000   16.0000000000    8.0000000000   17.0000000000

forward elimination
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.0000000000    7.5555555556   -2.2222222222    0.7777777778
  0.0000000000    0.0000000000    6.5882352941   -1.7058823529
  0.0000000000    0.0000000000    0.0000000000    5.3750000000
B
 20.0000000000   11.5555555556   12.9411764706   21.5000000000

X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

Ada

with TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions;
use  TEXT_IO, Ada.Long_Float_Text_IO, Ada.Numerics.Long_Elementary_Functions;

procedure Ada1003 is

    N : Constant Integer := 3;

    type Long_Float_Array       is array (0..N)       of Long_Float;
    type Long_Float_TwoDimArray is array (0..N, 0..N) of Long_Float;

    a:Long_Float_TwoDimArray := ((-1.0, -2.0, 7.0, -2.0), (1.0, -1.0, -2.0, 6.0), ( 9.0,  2.0, 1.0,  1.0), (2.0, 8.0, -2.0, 1.0));
    b:Long_Float_Array := (8.0, 17.0, 20.0, 16.0);

    -- １次元配列を表示
    procedure disp_vector(row:Long_Float_Array) is
    begin
        for i in row'Range loop
            Put(row(i), Fore=>3, Aft=>10, Exp=>0);
            Put(Ascii.HT);
        end loop;
        New_Line;
    end disp_vector;

    -- ２次元配列を表示
    procedure disp_matrix(matrix:Long_Float_TwoDimArray) is
    begin
        for i in matrix'Range loop
            for j in matrix'Range loop
                Put(matrix(i,j), Fore=>3, Aft=>10, Exp=>0);
                Put(Ascii.HT);
            end loop;
            New_Line;
        end loop;
    end disp_matrix;

    -- 前進消去
    procedure forward_elimination(a:in out Long_Float_TwoDimArray; b:in out Long_Float_Array) is
        s:Long_Float;
    begin
        for pivot in a'First..(a'Last - 1) loop
            for row in (pivot + 1)..a'Last loop
                s := a(row,pivot) / a(pivot,pivot);
                for col in pivot..a'Last loop
                    a(row,col) := a(row,col) - a(pivot,col) * s;
                end loop;
                b(row)         := b(row)     - b(pivot)     * s;
            end loop;
        end loop;
    end forward_elimination;

    -- 後退代入
    procedure backward_substitution(a:Long_Float_TwoDimArray; b:in out Long_Float_Array) is
    begin
        for row in reverse a'Range loop
            for col in reverse (row + 1)..a'Last loop
                b(row) := b(row) - a(row,col) * b(col);
            end loop;
            b(row) := b(row) / a(row,row);
        end loop;
    end backward_substitution;

    -- ピボット選択
    procedure pivoting(a:in out Long_Float_TwoDimArray; b:in out Long_Float_Array) is
        max_val, tmp:Long_Float;
        max_row:Integer;
    begin
        for pivot in a'Range loop
            -- 各列で 一番値が大きい行を 探す
            max_row :=   pivot;
            max_val :=   0.0;
            for row in  pivot..a'Last loop
                if Abs(a(row,pivot)) > max_val then
                    -- 一番値が大きい行
                    max_val :=   Abs(a(row,pivot));
                    max_row :=   row;
                end if;
            end loop;

            -- 一番値が大きい行と入れ替え
            if max_row /= pivot then
                tmp := 0.0;
                for col in a'Range loop
                    tmp            :=   a(max_row,col);
                    a(max_row,col) :=   a(pivot,col);
                    a(pivot,col)   :=   tmp;
                end loop;
                tmp        :=   b(max_row);
                b(max_row) :=   b(pivot);
                b(pivot)   :=   tmp;
            end if;
        end loop;
    end pivoting;

begin
    -- ピボット選択
    pivoting(a,b);

    Put_Line("pivoting");
    Put_Line("A");
    disp_matrix(a);
    Put_Line("B");
    disp_vector(b);
    Put_Line("");

    -- 前進消去
    forward_elimination(a,b);

    Put_Line("forward elimination");
    Put_Line("A");
    disp_matrix(a);
    Put_Line("B");
    disp_vector(b);
    Put_Line("");

    -- 後退代入
    backward_substitution(a,b);

    Put_Line("X");
    disp_vector(b);
end Ada1003;

xxxxxx@yyyyyy /Z
$ gnatmake Ada1003.adb

xxxxxx@yyyyyy /Z
$ Ada1003
pivoting
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  2.0000000000    8.0000000000   -2.0000000000    1.0000000000
 -1.0000000000   -2.0000000000    7.0000000000   -2.0000000000
  1.0000000000   -1.0000000000   -2.0000000000    6.0000000000
B
 20.0000000000   16.0000000000    8.0000000000   17.0000000000

forward elimination
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.0000000000    7.5555555556   -2.2222222222    0.7777777778
 -0.0000000000    0.0000000000    6.5882352941   -1.7058823529
  0.0000000000   -0.0000000000    0.0000000000    5.3750000000
B
 20.0000000000   11.5555555556   12.9411764706   21.5000000000

X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

VB.NET

Option Explicit

Module VB1003
    Private Const N As Integer = 3

    Public Sub Main()
        Dim a(,) As Double = {{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}}
        Dim b()  As Double = {8,17,20,16}

        'ピボット選択
        pivoting(a,b)

        Console.WriteLine("pivoting")
        Console.WriteLine("A")
        disp_matrix(a)
        Console.WriteLine("B")
        disp_vector(b)
        Console.WriteLine

        '前進消去
        forward_elimination(a,b)

        Console.WriteLine("forward elimination")
        Console.WriteLine("A")
        disp_matrix(a)
        Console.WriteLine("B")
        disp_vector(b)
        Console.WriteLine

        '後退代入
        backward_substitution(a,b)

        Console.WriteLine("X")
        disp_vector(b)
    End Sub

    '１次元配列の表示
    Private Sub disp_vector(ByVal row() As Double)
        For Each col As Double In row
            Console.Write(String.Format("{0,14:F10}{1}", col, vbTab))
        Next
        Console.WriteLine()
    End Sub

    '２次元配列の表示
    Private Sub disp_matrix(ByVal matrix(,) As Double)
        Dim i As Integer = 0
        For Each col As Double In matrix
            Console.Write(String.Format("{0,14:F10}{1}", col, vbTab))
            i += 1
            If i > N Then
                i = 0
                Console.WriteLine()
            End If
        Next
    End Sub

    '前進消去
    Private Sub forward_elimination(ByRef a(,) As Double, ByRef b() As Double)
        For pivot As Integer = 0 To (N - 1)
            For row As Integer = (pivot + 1) To N
                Dim s As Double = a(row,pivot) / a(pivot,pivot)
                For col As Integer = pivot To N
                    a(row,col) -= a(pivot,col) * s
                Next
                b(row) -= b(pivot) * s
            Next
        Next
    End Sub

    '後退代入
    Private Sub backward_substitution(ByVal a(,) As Double, ByRef b() As Double)
        For row As Integer = N To 0 Step -1
            For col As Integer = N To (row + 1) Step - 1
                b(row) -= a(row,col) * b(col)
            Next
            b(row) /= a(row,row)
        Next
    End Sub

    'ピボット選択
    Private Sub pivoting(ByRef a(,) As Double, ByRef b() As Double)
        For pivot As Integer = 0 To N
            '各列で 一番値が大きい行を 探す
            Dim max_row As Integer = pivot
            Dim max_val As Double = 0
            For row As Integer = pivot To N
                If Math.Abs(a(row,pivot)) > max_val Then
                    '一番値が大きい行
                    max_val =   Math.Abs(a(row,pivot))
                    max_row =   row
                End If
            Next

            '一番値が大きい行と入れ替え
            If max_row <> pivot Then
                Dim tmp As Double
                For col As Integer = 0 To N
                    tmp            =   a(max_row,col)
                    a(max_row,col) =   a(pivot,col)
                    a(pivot,col)   =   tmp
                Next
                tmp         =   b(max_row)
                b(max_row)  =   b(pivot)
                b(pivot)    =   tmp
            End If
        Next
    End Sub
End Module

Z:\>vbc -nologo VB1003.vb

Z:\>VB1003
pivoting
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  2.0000000000    8.0000000000   -2.0000000000    1.0000000000
 -1.0000000000   -2.0000000000    7.0000000000   -2.0000000000
  1.0000000000   -1.0000000000   -2.0000000000    6.0000000000
B
 20.0000000000   16.0000000000    8.0000000000   17.0000000000

forward elimination
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.0000000000    7.5555555556   -2.2222222222    0.7777777778
  0.0000000000    0.0000000000    6.5882352941   -1.7058823529
  0.0000000000    0.0000000000    0.0000000000    5.3750000000
B
 20.0000000000   11.5555555556   12.9411764706   21.5000000000

X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

C#

using System;

public class CS1003
{
    private const int N = 4;

    // ガウスの消去法
    public static void Main()
    {
        double[,] a = {{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}};
        double[]  b = {8,17,20,16};

        // ピボット選択
        pivoting(a,b);

        Console.WriteLine("pivoting");
        Console.WriteLine("A");
        disp_matrix(a);
        Console.WriteLine("B");
        disp_vector(b);
        Console.WriteLine();

        // 前進消去
        forward_elimination(a,b);

        Console.WriteLine("forward elimination");
        Console.WriteLine("A");
        disp_matrix(a);
        Console.WriteLine("B");
        disp_vector(b);
        Console.WriteLine();

        // 後退代入
        backward_substitution(a,b);

        Console.WriteLine("X");
        disp_vector(b);
    }
    // 前進消去
    private static void forward_elimination(double[,] a, double[] b)
    {
        for (int pivot = 0; pivot < N - 1; pivot++)
        {
            for (int row = pivot + 1; row < N; row++)
            {
                double  s   = a[row,pivot] / a[pivot,pivot];
                for (int col = pivot; col < N; col++)
                    a[row,col] -= a[pivot,col]    * s;
                b[row]         -= b[pivot]        * s;
            }
        }
    }
    // 後退代入
    private static void backward_substitution(double[,] a, double[] b)
    {
        for (int row = N - 1; row >= 0; row--)
        {
            for (int col = N - 1; col > row; col--)
                b[row] -= a[row,col] * b[col];
            b[row] /= a[row,row];
        }
    }
    // ピボット選択
    private static void pivoting(double[,] a, double[] b)
    {
        for(int pivot = 0; pivot < N; pivot++)
        {
            // 各列で 一番値が大きい行を 探す
            int     max_row =   pivot;
            double  max_val =   0;
            for (int row = pivot; row < N; row++)
            {
                if (Math.Abs(a[row,pivot]) > max_val)
                {
                    // 一番値が大きい行
                    max_val =   Math.Abs(a[row,pivot]);
                    max_row =   row;
                }
            }

            // 一番値が大きい行と入れ替え
            if (max_row != pivot) {
                double tmp;
                for (int col = 0; col < N; col++)
                {
                    tmp             =   a[max_row,col];
                    a[max_row,col] =   a[pivot,col];
                    a[pivot,col]   =   tmp;
                }
                tmp         =   b[max_row];
                b[max_row]  =   b[pivot];
                b[pivot]    =   tmp;
            }
        }
    }
    // １次元配列を表示
    private static void disp_vector(double[] row)
    {
        foreach (double col in row)
            Console.Write(string.Format("{0,14:F10}\t", col));
        Console.WriteLine();
    }
    // ２次元配列を表示
    private static void disp_matrix(double[,] matrix)
    {
        int i = 0;
        foreach (double col in matrix)
        {
            Console.Write(string.Format("{0,14:F10}\t", col));
            if (++i >= N)
            {
                i = 0;
                Console.WriteLine();
            }
        }
    }
}

Z:\>csc -nologo CS1003.cs

Z:\>CS1003
pivoting
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  2.0000000000    8.0000000000   -2.0000000000    1.0000000000
 -1.0000000000   -2.0000000000    7.0000000000   -2.0000000000
  1.0000000000   -1.0000000000   -2.0000000000    6.0000000000
B
 20.0000000000   16.0000000000    8.0000000000   17.0000000000

forward elimination
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.0000000000    7.5555555556   -2.2222222222    0.7777777778
  0.0000000000    0.0000000000    6.5882352941   -1.7058823529
  0.0000000000    0.0000000000    0.0000000000    5.3750000000
B
 20.0000000000   11.5555555556   12.9411764706   21.5000000000

X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

Java

import java.lang.*;

public class Java1003 {

    private static final int N = 4;

    // ガウスの消去法
    public static void main(String []args) {
        double[][] a = {{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}};
        double[]   b = {8,17,20,16};

        // ピボット選択
        pivoting(a,b);

        System.out.println("pivoting");
        System.out.println("A");
        disp_matrix(a);
        System.out.println("B");
        disp_vector(b);
        System.out.println();

        // 前進消去
        forward_elimination(a,b);

        System.out.println("forward elimination");
        System.out.println("A");
        disp_matrix(a);
        System.out.println("B");
        disp_vector(b);
        System.out.println();

        // 後退代入
        backward_substitution(a,b);

        System.out.println("X");
        disp_vector(b);
    }
    // 前進消去
    private static void forward_elimination(double[][] a, double[] b) {
        for (int pivot = 0; pivot < N - 1; pivot++) {
            for (int row = pivot + 1; row < N; row++) {
                double  s   = a[row][pivot] / a[pivot][pivot];
                for (int col = pivot; col < N; col++)
                    a[row][col] -= a[pivot][col]    * s;
                b[row]          -= b[pivot]         * s;
            }
        }
    }
    // 後退代入
    private static void backward_substitution(double[][] a, double[] b) {
        for (int row = N - 1; row >= 0; row--) {
            for (int col = N - 1; col > row; col--)
                b[row] -= a[row][col] * b[col];
            b[row] /= a[row][row];
        }
    }
    // ピボット選択
    private static void pivoting(double[][] a, double[] b) {
        for(int pivot = 0; pivot < N; pivot++) {
            // 各列で 一番値が大きい行を 探す
            int     max_row =   pivot;
            double  max_val =   0;
            for (int row = pivot; row < N; row++) {
                if (Math.abs(a[row][pivot]) > max_val) {
                    // 一番値が大きい行
                    max_val =   Math.abs(a[row][pivot]);
                    max_row =   row;
                }
            }

            // 一番値が大きい行と入れ替え
            if (max_row != pivot) {
                double tmp;
                for (int col = 0; col < N; col++) {
                    tmp             =   a[max_row][col];
                    a[max_row][col] =   a[pivot][col];
                    a[pivot][col]   =   tmp;
                }
                tmp         =   b[max_row];
                b[max_row]  =   b[pivot];
                b[pivot]    =   tmp;
            }
        }
    }
    // １次元配列を表示
    private static void disp_vector(double[] row) {
        for (double col: row)
            System.out.print(String.format("%14.10f\t", col));
        System.out.println();
    }
    // ２次元配列を表示
    private static void disp_matrix(double[][] matrix) {
        for (double[] row: matrix) {
            for (double col: row)
                System.out.print(String.format("%14.10f\t", col));
            System.out.println();
        }
    }
}

Z:\>javac Java1003.java

Z:\>java Java1003
pivoting
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  2.0000000000    8.0000000000   -2.0000000000    1.0000000000
 -1.0000000000   -2.0000000000    7.0000000000   -2.0000000000
  1.0000000000   -1.0000000000   -2.0000000000    6.0000000000
B
 20.0000000000   16.0000000000    8.0000000000   17.0000000000

forward elimination
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.0000000000    7.5555555556   -2.2222222222    0.7777777778
  0.0000000000    0.0000000000    6.5882352941   -1.7058823529
  0.0000000000    0.0000000000    0.0000000000    5.3750000000
B
 20.0000000000   11.5555555556   12.9411764706   21.5000000000

X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

C++

#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;

const int N = 4;

// ピボット選択
void pivoting(double a[N][N], double b[N]);
// 前進消去
void forward_elimination(double a[N][N], double b[N]);
// 後退代入
void backward_substitution(double a[N][N], double b[N]);
// １次元配列を表示
void disp_vector(double row[N]);
// ２次元配列を表示
void disp_matrix(double matrix[N][N]);

// ガウスの消去法
int main()
{
    double a[N][N] = {{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}};
    double b[N]    = {8,17,20,16};

    // ピボット選択
    pivoting(a,b);

    cout << "pivoting" << endl;
    cout << "A" << endl;
    disp_matrix(a);
    cout << "B" << endl;
    disp_vector(b);
    cout << endl;

    // 前進消去
    forward_elimination(a,b);

    cout << "forward elimination" << endl;
    cout << "A" << endl;
    disp_matrix(a);
    cout << "B" << endl;
    disp_vector(b);
    cout << endl;

    // 後退代入
    backward_substitution(a,b);

    cout << "X" << endl;
    disp_vector(b);

    return 0;
}
// 前進消去
void forward_elimination(double a[N][N], double b[N])
{
    for (int pivot = 0; pivot < N - 1; pivot++)
    {
        for (int row = pivot + 1; row < N; row++)
        {
            double  s   = a[row][pivot] / a[pivot][pivot];
            for (int col = pivot; col < N; col++)
                a[row][col] -= a[pivot][col]    * s;
            b[row]          -= b[pivot]         * s;
        }
    }
}
// 後退代入
void backward_substitution(double a[N][N], double b[N])
{
    for (int row = N - 1; row >= 0; row--)
    {
        for (int col = N - 1; col > row; col--)
            b[row] -= a[row][col] * b[col];
        b[row] /= a[row][row];
    }
}
// ピボット選択
void pivoting(double a[N][N], double b[N])
{
    for(int pivot = 0; pivot < N; pivot++)
    {
        // 各列で 一番値が大きい行を 探す
        int     max_row =   pivot;
        double  max_val =   0;
        for (int row = pivot; row < N; row++)
        {
            if (fabs(a[row][pivot]) > max_val)
            {
                // 一番値が大きい行
                max_val =   fabs(a[row][pivot]);
                max_row =   row;
            }
        }

        // 一番値が大きい行と入れ替え
        if (max_row != pivot)
        {
            double tmp;
            for (int col = 0; col < N; col++)
            {
                tmp             =   a[max_row][col];
                a[max_row][col] =   a[pivot][col];
                a[pivot][col]   =   tmp;
            }
            tmp         =   b[max_row];
            b[max_row]  =   b[pivot];
            b[pivot]    =   tmp;
        }
    }
}
// １次元配列を表示
void disp_vector(double row[N])
{
    for (int i = 0; i < N; i++)
        cout << setw(14) << fixed << setprecision(10) << row[i] << "\t";
    cout << endl;
}
// ２次元配列を表示
void disp_matrix(double matrix[N][N])
{
    for (int i = 0; i < N; i++)
    {
        for (int j = 0; j < N; j++)
            cout << setw(14) << fixed << setprecision(10) << matrix[i][j] << "\t";
        cout << endl;
    }
}

Z:\>bcc32 -q CP1003.cpp
cp1003.cpp:

Z:\>CP1003
pivoting
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  2.0000000000    8.0000000000   -2.0000000000    1.0000000000
 -1.0000000000   -2.0000000000    7.0000000000   -2.0000000000
  1.0000000000   -1.0000000000   -2.0000000000    6.0000000000
B
 20.0000000000   16.0000000000    8.0000000000   17.0000000000

forward elimination
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.0000000000    7.5555555556   -2.2222222222    0.7777777778
 -0.0000000000    0.0000000000    6.5882352941   -1.7058823529
  0.0000000000   -0.0000000000    0.0000000000    5.3750000000
B
 20.0000000000   11.5555555556   12.9411764706   21.5000000000

X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

Objective-C

#import <Foundation/Foundation.h>
#import <math.h>

const int N = 4;

// ピボット選択
void pivoting(double a[N][N], double b[N]);
// 前進消去
void forward_elimination(double a[N][N], double b[N]);
// 後退代入
void backward_substitution(double a[N][N], double b[N]);
// １次元配列を表示
void disp_vector(double row[N]);
// ２次元配列を表示
void disp_matrix(double matrix[N][N]);

// ガウスの消去法
int main()
{
    double a[4][4] = {{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}};
    double b[4]    = {8,17,20,16};

    // ピボット選択
    pivoting(a,b);

    printf("pivoting\n");
    printf("A\n");
    disp_matrix(a);
    printf("B\n");
    disp_vector(b);
    printf("\n");

    // 前進消去
    forward_elimination(a,b);

    printf("forward elimination\n");
    printf("A\n");
    disp_matrix(a);
    printf("B\n");
    disp_vector(b);
    printf("\n");

    // 後退代入
    backward_substitution(a,b);

    printf("X\n");
    disp_vector(b);

    return 0;
}
// 前進消去
void forward_elimination(double a[N][N], double b[N])
{
    int pivot, row, col;
    for (pivot = 0; pivot < N - 1; pivot++)
    {
        for (row = pivot + 1; row < N; row++)
        {
            double  s   = a[row][pivot] / a[pivot][pivot];
            for (col = pivot; col < N; col++)
                a[row][col] -= a[pivot][col]    * s;
            b[row]          -= b[pivot]         * s;
        }
    }
}
// 後退代入
void backward_substitution(double a[N][N], double b[N])
{
    int row, col;
    for (row = N - 1; row >= 0; row--)
    {
        for (col = N - 1; col > row; col--)
            b[row] -= a[row][col] * b[col];
        b[row] /= a[row][row];
    }
}
// ピボット選択
void pivoting(double a[N][N], double b[N])
{
    int pivot, row, col;
    for(pivot = 0; pivot < N; pivot++)
    {
        // 各列で 一番値が大きい行を 探す
        int     max_row =   pivot;
        double  max_val =   0;
        for (row = pivot; row < N; row++)
        {
            if (fabs(a[row][pivot]) > max_val)
            {
                // 一番値が大きい行
                max_val =   fabs(a[row][pivot]);
                max_row =   row;
            }
        }

        // 一番値が大きい行と入れ替え
        if (max_row != pivot)
        {
            double tmp;
            for (col = 0; col < N; col++)
            {
                tmp             =   a[max_row][col];
                a[max_row][col] =   a[pivot][col];
                a[pivot][col]   =   tmp;
            }
            tmp         =   b[max_row];
            b[max_row]  =   b[pivot];
            b[pivot]    =   tmp;
        }
    }
}
// １次元配列を表示
void disp_vector(double row[N])
{
    int i;
    for (i = 0; i < N; i++)
    {
        printf("%14.10f\t", row[i]);
    }
    printf("\n");
}
// ２次元配列を表示
void disp_matrix(double matrix[N][N])
{
    int i, j;
    for (i = 0; i < N; i++)
    {
        for (j = 0; j < N; j++)
            printf("%14.10f\t", matrix[i][j]);
        printf("\n");
    }
}

xxxxxx@yyyyyy /Z
$ gcc -o OC1003 OC1003.m -lobjc -lgnustep-base -I $INCLUDE -L $LIB $CFLAGS

xxxxxx@yyyyyy /Z
$ OC1003
pivoting
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  2.0000000000    8.0000000000   -2.0000000000    1.0000000000
 -1.0000000000   -2.0000000000    7.0000000000   -2.0000000000
  1.0000000000   -1.0000000000   -2.0000000000    6.0000000000
B
 20.0000000000   16.0000000000    8.0000000000   17.0000000000

forward elimination
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.0000000000    7.5555555556   -2.2222222222    0.7777777778
 -0.0000000000    0.0000000000    6.5882352941   -1.7058823529
  0.0000000000   -0.0000000000    0.0000000000    5.3750000000
B
 20.0000000000   11.5555555556   12.9411764706   21.5000000000

X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

D

import std.stdio;
import std.math;

const int N = 4;

// ガウスの消去法
void main(string[] args)
{
    double[N][N] a = [[-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1]];
    double[N]    b = [8,17,20,16];

    // ピボット選択
    pivoting(a,b);

    writefln("pivoting");
    writefln("A");
    disp_matrix(a);
    writefln("B");
    disp_vector(b);
    writefln("");

    // 前進消去
    forward_elimination(a,b);

    writefln("forward elimination");
    writefln("A");
    disp_matrix(a);
    writefln("B");
    disp_vector(b);
    writefln("");

    // 後退代入
    backward_substitution(a,b);

    writefln("X");
    disp_vector(b);
}
// 前進消去
void forward_elimination(ref double[N][N] a, ref double[N] b)
{
    foreach (pivot; 0..N)
    {
        foreach (row; (pivot + 1)..N)
        {
            double s = a[row][pivot] / a[pivot][pivot];
            foreach (col; pivot..N)
                a[row][col] -= a[pivot][col]    * s;
            b[row]          -= b[pivot]         * s;
        }
    }
}
// 後退代入
void backward_substitution(ref double[N][N] a, ref double[N] b)
{
    foreach_reverse (row; 0..N)
    {
        foreach_reverse (col; (row+1)..N)
            b[row] -= a[row][col] * b[col];
        b[row] /= a[row][row];
    }
}
// ピボット選択
void pivoting(ref double[N][N] a, ref double[N] b)
{
    foreach (pivot; 0..N)
    {
        // 各列で 一番値が大きい行を 探す
        int     max_row =   pivot;
        double  max_val =   0;
        foreach (row; pivot..N)
        {
            if (fabs(a[row][pivot]) > max_val)
            {
                // 一番値が大きい行
                max_val =   fabs(a[row][pivot]);
                max_row =   row;
            }
        }

        // 一番値が大きい行と入れ替え
        if (max_row != pivot)
        {
            double tmp;
            foreach (col; 0..N)
            {
                tmp             =   a[max_row][col];
                a[max_row][col] =   a[pivot][col];
                a[pivot][col]   =   tmp;
            }
            tmp         =   b[max_row];
            b[max_row]  =   b[pivot];
            b[pivot]    =   tmp;
        }
    }
}
// １次元配列を表示
void disp_vector(double[] row)
{
    foreach(col; row)
        writef("%14.10f\t", col);
    writefln("");
}
// ２次元配列を表示
void disp_matrix(double[N][N] matrix)
{
    foreach(row; matrix)
    {
        foreach(col; row)
            writef("%14.10f\t", col);
        writefln("");
    }
}

Z:\>dmd D1003.d

Z:\>D1003
pivoting
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  2.0000000000    8.0000000000   -2.0000000000    1.0000000000
 -1.0000000000   -2.0000000000    7.0000000000   -2.0000000000
  1.0000000000   -1.0000000000   -2.0000000000    6.0000000000
B
 20.0000000000   16.0000000000    8.0000000000   17.0000000000

forward elimination
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.0000000000    7.5555555556   -2.2222222222    0.7777777778
  0.0000000000    0.0000000000    6.5882352941   -1.7058823529
  0.0000000000    0.0000000000    0.0000000000    5.3750000000
B
 20.0000000000   11.5555555556   12.9411764706   21.5000000000

X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

Go

package main

import "fmt"
import "math"

const N = 4

// ガウスの消去法
func main() {
    var a [N][N]float64 = [N][N]float64{{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}}
    var b []float64     = []float64{8,17,20,16}

    // ピボット選択
    pivoting(&a,b)

    fmt.Println("pivoting")
    fmt.Println("A")
    disp_matrix(a);
    fmt.Println("B")
    disp_vector(b);
    fmt.Println()

    // 前進消去
    forward_elimination(&a,b)

    fmt.Println("forward elimination")
    fmt.Println("A")
    disp_matrix(a);
    fmt.Println("B")
    disp_vector(b);
    fmt.Println()

    // 後退代入
    backward_substitution(&a,b)

    fmt.Println("X")
    disp_vector(b);
}
// 前進消去
func forward_elimination(a *[N][N]float64, b []float64) {
    for pivot := 0; pivot < N - 1; pivot++ {
        for row := pivot + 1; row < N; row++ {
            var s = a[row][pivot] / a[pivot][pivot]
            for col := pivot; col < N; col++ {
                a[row][col] -= a[pivot][col]    * s
            }
            b[row]          -= b[pivot]         * s
        }
    }
}
// 後退代入
func backward_substitution(a *[N][N]float64, b []float64) {
    for row := N - 1; row >= 0; row-- {
        for col := N - 1; col > row; col-- {
            b[row] -= a[row][col] * b[col]
        }
        b[row] /= a[row][row]
    }
}
// ピボット選択
func pivoting(a *[N][N]float64, b []float64) {
    for pivot := 0; pivot < N; pivot++ {
        // 各列で 一番値が大きい行を 探す
        var max_row         = pivot
        var max_val float64 = 0.0
        for row := pivot; row < N; row++ {
            if math.Fabs(a[row][pivot]) > max_val {
                // 一番値が大きい行
                max_val =   math.Fabs(a[row][pivot])
                max_row =   row
            }
        }

        // 一番値が大きい行と入れ替え
        if (max_row != pivot) {
            var tmp float64 = 0.0
            for col := 0; col < N; col++ {
                tmp             =   a[max_row][col]
                a[max_row][col] =   a[pivot][col]
                a[pivot][col]   =   tmp
            }
            tmp         =   b[max_row]
            b[max_row]  =   b[pivot]
            b[pivot]    =   tmp
        }
    }
}
// １次元配列を表示
func disp_vector(row[]float64) {
    for _, col := range row {
        fmt.Printf("%14.10f\t", col)
    }
    fmt.Println("")
}
// ２次元配列を表示
func disp_matrix(matrix[N][N]float64) {
    for _, row := range matrix {
        for _, col := range row {
            fmt.Printf("%14.10f\t", col)
        }
        fmt.Println("")
    }
}

Z:\>8g GO1003.go

Z:\>8l -o GO1003.exe GO1003.8

Z:\>GO1003
pivoting
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  2.0000000000    8.0000000000   -2.0000000000    1.0000000000
 -1.0000000000   -2.0000000000    7.0000000000   -2.0000000000
  1.0000000000   -1.0000000000   -2.0000000000    6.0000000000
B
 20.0000000000   16.0000000000    8.0000000000   17.0000000000

forward elimination
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.0000000000    7.5555555556   -2.2222222222    0.7777777778
  0.0000000000    0.0000000000    6.5882352941   -1.7058823529
  0.0000000000    0.0000000000    0.0000000000    5.3750000000
B
 20.0000000000   11.5555555556   12.9411764706   21.5000000000

X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

Scala

object Scala1003 {
    val N = 3

    def main(args: Array[String]) {
        var a:Array[Array[Double]] = Array(Array(-1,-2,7,-2),Array(1,-1,-2,6),Array(9,2,1,1),Array(2,8,-2,1))
        var b:Array[Double] = Array(8,17,20,16)

        // ピボット選択
        pivoting(a,b)

        println("pivoting")
        println("A")
        disp_matrix(a)
        println("B")
        disp_vector(b)
        println()

        // 前進消去
        forward_elimination(a,b)

        println("forward elimination")
        println("A")
        disp_matrix(a)
        println("B")
        disp_vector(b)
        println()

        // 後退代入
        backward_substitution(a,b)

        println("X")
        disp_vector(b)
    }
    // 前進消去
    def forward_elimination(a:Array[Array[Double]], b:Array[Double]) {
        for (pivot <- 0 to N - 1) {
            for (row <- pivot + 1 to N) {
                var s:Double = a(row)(pivot) / a(pivot)(pivot)
                for (col <- pivot to N)
                    a(row)(col) -= a(pivot)(col)    * s
                b(row)          -= b(pivot)         * s
            }
        }
    }
    // 後退代入
    def backward_substitution(a:Array[Array[Double]], b:Array[Double]) {
        for (row <- (0 to N).reverse) {
            for (col <- (row + 1 to N).reverse)
                b(row) -= a(row)(col) * b(col)
            b(row) /= a(row)(row)
        }
    }
    // ピボット選択
    def pivoting(a:Array[Array[Double]], b:Array[Double]) {
        for (pivot <- 0 to N) {
            // 各列で 一番値が大きい行を 探す
            var max_row:Integer =   pivot
            var max_val:Double  =   0
            for (row <- pivot to N) {
                if (Math.abs(a(row)(pivot)) > max_val) {
                    // 一番値が大きい行
                    max_val =   Math.abs(a(row)(pivot))
                    max_row =   row
                }
            }

            // 一番値が大きい行と入れ替え
            if (max_row != pivot) {
                var tmp:Double = 0
                for (col <- 0 to N) {
                    tmp             =   a(max_row)(col)
                    a(max_row)(col) =   a(pivot)(col)
                    a(pivot)(col)   =   tmp
                }
                tmp         =   b(max_row)
                b(max_row)  =   b(pivot)
                b(pivot)    =   tmp
            }
        }
    }

    // １次元配列を表示
    def disp_vector(row:Array[Double]) = {
        row.foreach { col =>
            print("%14.10f\t".format(col))
        }
        println()
    }
    // ２次元配列を表示
    def disp_matrix(matrix:Array[Array[Double]]) = {
        matrix.foreach { row =>
            row.foreach { col =>
                print("%14.10f\t".format(col))
            }
            println()
        }
    }
}

Z:\>scala Scala1003.scala
pivoting
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  2.0000000000    8.0000000000   -2.0000000000    1.0000000000
 -1.0000000000   -2.0000000000    7.0000000000   -2.0000000000
  1.0000000000   -1.0000000000   -2.0000000000    6.0000000000
B
 20.0000000000   16.0000000000    8.0000000000   17.0000000000

forward elimination
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.0000000000    7.5555555556   -2.2222222222    0.7777777778
  0.0000000000    0.0000000000    6.5882352941   -1.7058823529
  0.0000000000    0.0000000000    0.0000000000    5.3750000000
B
 20.0000000000   11.5555555556   12.9411764706   21.5000000000

X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

F#

module Fs1003

open System

let N = 3

// １次元配列を表示
let disp_vector (row:float[]) =
    row
    |> Array.iter (fun x -> printf "%14.10f" x)
    printfn ""

// ２次元配列を表示
let disp_matrix (matrix:float[][]) =
    matrix
    |> Array.iter
        (fun row ->
            row
            |> Array.iter (fun col -> printf "%14.10f" col)
            printfn ""
        )

// 前進消去
let forward_elimination (a:float[][]) (b:float[]) =
    for pivot in [0..N - 1] do
        for row in [pivot + 1 .. N] do
            let s:float = a.[row].[pivot] / a.[pivot].[pivot]
            for col in [pivot..N] do
                a.[row].[col] <- a.[row].[col] - a.[pivot].[col]    * s
            b.[row]           <- b.[row]       - b.[pivot]          * s

// 後退代入
let backward_substitution (a:float[][]) (b:float[]) =
    for row in (List.rev [0..N]) do
        for col in (List.rev [row + 1 .. N]) do
            b.[row] <- b.[row] - a.[row].[col] * b.[col]
        b.[row] <- b.[row] / a.[row].[row]

// ピボット選択
let pivoting (a:float[][]) (b:float[]) =
    for pivot in [0..N] do
        // 各列で 一番値が大きい行を 探す
        let mutable max_row:int = pivot
        let mutable max_val:float = 0.0
        for row in [pivot..N] do
            if (Math.Abs(a.[row].[pivot]) > max_val) then
                // 一番値が大きい行
                max_val <- Math.Abs(a.[row].[pivot])
                max_row <- row

        // 一番値が大きい行と入れ替え
        if (max_row <> pivot) then
            let mutable tmp:float = 0.0
            for col in [0..N] do
                tmp               <-   a.[max_row].[col]
                a.[max_row].[col] <-   a.[pivot].[col]
                a.[pivot].[col]   <-   tmp
            tmp         <-   b.[max_row]
            b.[max_row] <-   b.[pivot]
            b.[pivot]   <-   tmp

let mutable a:float[][] = [| [|-1.0;-2.0;7.0;-2.0 |]; [|1.0;-1.0;-2.0;6.0 |]; [| 9.0;2.0;1.0;1.0 |]; [|2.0;8.0;-2.0;1.0 |] |]
let mutable b:float[] = [| 8.0;17.0;20.0;16.0 |]

// ピボット選択
pivoting a b

printfn "pivoting"
printfn "A"
disp_matrix a
printfn "B"
disp_vector b
printfn ""

// 前進消去
forward_elimination a b

printfn "forward elimination"
printfn "A"
disp_matrix a
printfn "B"
disp_vector b
printfn ""

// 後退代入
backward_substitution a b

printfn "X"
disp_vector b

exit 0

Z:\>fsi  --nologo --quiet Fs1003.fs
pivoting
A
  9.0000000000  2.0000000000  1.0000000000  1.0000000000
  2.0000000000  8.0000000000 -2.0000000000  1.0000000000
 -1.0000000000 -2.0000000000  7.0000000000 -2.0000000000
  1.0000000000 -1.0000000000 -2.0000000000  6.0000000000
B
 20.0000000000 16.0000000000  8.0000000000 17.0000000000

forward elimination
A
  9.0000000000  2.0000000000  1.0000000000  1.0000000000
  0.0000000000  7.5555555556 -2.2222222222  0.7777777778
  0.0000000000  0.0000000000  6.5882352941 -1.7058823529
  0.0000000000  0.0000000000  0.0000000000  5.3750000000
B
 20.0000000000 11.5555555556 12.9411764706 21.5000000000

X
  1.0000000000  2.0000000000  3.0000000000  4.0000000000
TOP

Clojure

(def N 4)

(def a [[-1.0 -2.0 7.0 -2.0] [1.0 -1.0 -2.0 6.0] [9.0 2.0 1.0 1.0] [2.0 8.0 -2.0 1.0]])
(def b  [8.0 17.0 20.0 16.0])

;１次元配列を表示
(defn disp_vector [row]
    (doseq [col row]
        (print (format "%14.10f\t" col)))
    (println))

;２次元配列を表示
(defn disp_matrix [matrix]
    (doseq [row matrix]
        (doseq [col row]
            (print (format "%14.10f\t" col)))
        (println)))

;各列で 一番値が大きい行を 探す
(defn get_max_row [row col a max_row max_val]
    ;一番値が大きい行
    (def ws1
        (if (< max_val (. Math abs (nth (nth a row) col)))
            (vector row     (. Math abs (nth (nth a row) col)))
            (vector max_row max_val)))

    (def max_row2 (first  ws1))
    (def max_val2 (second ws1))

    (if (>= row (dec (count a)))
        (vector max_row2 max_val2)
        (get_max_row (inc row) col a max_row2 max_val2)))

;ピボット選択
(defn pivoting [pivot a b a2 b2]
    (def ws1 (get_max_row 0 pivot a 0 0.0))
    (def max_row (first  ws1))
    (def max_val (second ws1))

    (def a3 (cons (nth a max_row) a2))
    (def b3 (cons (nth b max_row) b2))
    (def a4 (concat (take max_row a) (drop (inc max_row) a)))
    (def b4 (concat (take max_row b) (drop (inc max_row) b)))

    (if (>= pivot (dec N))
        (vector (reverse a3) (reverse b3))
        (pivoting (inc pivot) a4 b4 a3 b3)))

;前進消去
(defn forward_elim_loop [pivot row a b]
    (def s (/ (nth (nth a row) pivot) (nth (nth a pivot) pivot)))

    (def a_map (map
        (fn [a_piv a_row] (- a_row (* a_piv s)))
            (drop pivot (nth a pivot)) (drop pivot (nth a row))))
    (def a1 (concat (take pivot (nth a row)) a_map))
    (def a2 (concat (take row   a)           (cons (vec a1) (drop (inc row) a))))

    (def x (- (nth b row) (* (nth b pivot) s)))
    (def b2 (concat (take row b) (cons x (drop (inc row) b))))
    (if (< row (dec N))
        (forward_elim_loop pivot (inc row) a2 b2)
        (vector a2 b2)))

(defn forward_elimination [pivot a b]
    (def ws1
        (if (< pivot (dec N))
            (forward_elim_loop pivot (inc pivot) a b)
            (vector a b)))
    (def a2 (first ws1))
    (def b2 (second ws1))

    (if (< pivot (dec N))
        (forward_elimination (inc pivot) a2 b2)
        (vector a2 b2)))

;後退代入
(defn backward_substitution [row a b]
    (def a_map (map
        (fn [a_col b_col] (* a_col b_col))
            (drop (inc row) (nth a row)) (drop (inc row) b)))
    (def s (apply + a_map))
    (def x (/ (- (nth b row) s) (nth (nth a row) row)))
    (def b2 (concat (take row b) (cons x (drop (inc row) b))))
    (if (> row 0)
        (cons x (backward_substitution (dec row) a b2))
        (cons x [])))

;ピボット選択
(def ws1 (pivoting 0 a b [] []))
(def a1 (first ws1))
(def b1 (vec (second ws1)))

(println "pivoting")
(println "A")
(disp_matrix a1)
(println "B")
(disp_vector b1)
(println "")

;前進消去
(def ws2 (forward_elimination 0 a1 b1))
(def a2 (first ws2))
(def b2 (vec (second ws2)))

(println "forward elimination")
(println "A")
(disp_matrix a2)
(println "B")
(disp_vector b2)
(println "")

;後退代入
(def x (backward_substitution (dec N) a2 b2))

(println "X")
(disp_vector (reverse x))

Z:\>java -cp C:\ProgramFiles\clojure-1.5.1\clojure-1.5.1.jar clojure.main Clj1003.clj
pivoting
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  2.0000000000    8.0000000000   -2.0000000000    1.0000000000
 -1.0000000000   -2.0000000000    7.0000000000   -2.0000000000
  1.0000000000   -1.0000000000   -2.0000000000    6.0000000000
B
 20.0000000000   16.0000000000    8.0000000000   17.0000000000

forward elimination
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.0000000000    7.5555555556   -2.2222222222    0.7777777778
  0.0000000000    0.0000000000    6.5882352941   -1.7058823529
  0.0000000000    0.0000000000    0.0000000000    5.3750000000
B
 20.0000000000   11.5555555556   12.9411764706   21.5000000000

X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

Haskell

import Text.Printf
import Debug.Trace
import Control.Monad

n = 4::Int

-- ２次元配列を表示
disp_matrix::[[Double]]->IO()
disp_matrix matrix = do
    forM_ matrix $ \row -> do
        forM_ row $ \elem -> do
            printf "%14.10f\t" elem
        putStrLn ""

-- １次元配列を表示
disp_vector::[Double]->IO()
disp_vector vector = do
    forM_ vector $ \elem -> do
        printf "%14.10f\t" elem
    putStrLn ""

-- 各列で 一番値が大きい行を 探す
get_max_row::Int->Int->[[Double]]->Int->Double->(Int, Double)
get_max_row row col a max_row max_val =
    let
        -- 一番値が大きい行
        (max_row2, max_val2) =
            if max_val < abs(a!!row!!col)
                then (row,     abs(a!!row!!col))
                else (max_row, max_val)
    in
        if row >= length(a) - 1
            then
                (max_row2, max_val2)
            else
                get_max_row (row+1) col a max_row2 max_val2

-- ピボット選択
pivoting::Int->[[Double]]->[Double]->[[Double]]->[Double]->([[Double]],[Double])
pivoting pivot a b a2 b2 =
    let
        (max_row, max_val) = get_max_row 0 pivot a 0 0.0
        a3 = (a!!max_row):a2
        b3 = (b!!max_row):b2
        a4 = (take max_row a) ++ (drop (max_row+1) a)
        b4 = (take max_row b) ++ (drop (max_row+1) b)
    in
        if pivot >= (n - 1)
            then
                (reverse a3, reverse b3)
            else
                pivoting (pivot+1) a4 b4 a3 b3

-- 前進消去
forward_elim_loop::Int->Int->[[Double]]->[Double]->([[Double]],[Double])
forward_elim_loop pivot row a b =
    let
        s = a!!row!!pivot / a!!pivot!!pivot

        a_zip   = zip (drop pivot (a!!pivot)) (drop pivot (a!!row))
        a_map   = map (\(a_pivot, a_row) -> a_row - a_pivot * s) a_zip
        new_row = (take pivot (a!!row)) ++ a_map
        a2      = (take row   a       ) ++ (new_row:(drop (row+1) a))

        x  = b!!row - b!!pivot * s
        b2 = (take row b) ++ (x:(drop (row+1) b))
    in
        if row < (n - 1)
            then
                forward_elim_loop pivot (row+1) a2 b2
            else
                (a2, b2)

forward_elimination::Int->[[Double]]->[Double]->([[Double]],[Double])
forward_elimination pivot a b =
    let
        (a2, b2) =
            if pivot < (n - 1)
                then forward_elim_loop pivot (pivot+1) a b
                else (a, b)
    in
        if pivot < (n - 1)
            then
                forward_elimination (pivot+1) a2 b2
            else
                (a2, b2)

-- 後退代入
backward_substitution::Int->[[Double]]->[Double]->[Double]
backward_substitution row a b =
    let
        a_zip = zip (drop (row+1) (a!!row)) (drop (row+1) b)
        s  = sum $ map(\(a_col, b_col) -> a_col * b_col) a_zip
        x  = (b!!row - s) / (a!!row!!row)
        b2 = (take row b) ++ (x:(drop (row+1) b))
    in
        if row > 0
            then x:(backward_substitution (row-1) a b2)
            else x:[]

main = do
    let a  = [[-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1::Double]]
    let b  = [8,17,20,16::Double]

    -- ピボット選択
    let (a1, b1) = pivoting 0 a b [] []

    putStrLn "pivoting"
    putStrLn "A"
    disp_matrix a1
    putStrLn "B"
    disp_vector b1
    putStrLn ""

    -- 前進消去
    let (a2, b2) = forward_elimination 0 a1 b1

    putStrLn "forward elimination"
    putStrLn "A"
    disp_matrix a2
    putStrLn "B"
    disp_vector b2
    putStrLn ""

    -- 後退代入
    let x = backward_substitution (n-1) a2 b2

    putStrLn "X"
    disp_vector (reverse x)

Z:\>runghc Hs1003.hs
pivoting
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  2.0000000000    8.0000000000   -2.0000000000    1.0000000000
 -1.0000000000   -2.0000000000    7.0000000000   -2.0000000000
  1.0000000000   -1.0000000000   -2.0000000000    6.0000000000
B
 20.0000000000   16.0000000000    8.0000000000   17.0000000000

forward elimination
A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.0000000000    7.5555555556   -2.2222222222    0.7777777778
  0.0000000000    0.0000000000    6.5882352941   -1.7058823529
  0.0000000000    0.0000000000    0.0000000000    5.3750000000
B
 20.0000000000   11.5555555556   12.9411764706   21.5000000000

X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP
inserted by FC2 system