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

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

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

LU分解法

LU分解法を使って, 連立一次方程式の解を求める .

連立一次方程式

を行列で,

と表す. このとき

とすると, 最初の連立方程式は $ Ax = b $ と表すことができる.
$ LU = A $ となる上三角行列 $U$, 下三角行列 $L$ を考えると,
$ LUx = b $ であり, $ Ux = y $ とおくと $ Ly = b $ となる.

LU分解法では, 係数行列 $A$ を上三角行列 $U$, 下三角行列 $L$ に分解し,
$ Ly = b $ から $y$ を求め, $ Ux = y $ から $x$ を求める.

上三角行列 $U$ は, ガウスの消去法で使った上三角行列そのもので,当然求め方も同じ.
下三角行列 $L$ は, ガウスの消去法の計算途中で使った値を保存しておくだけ.
手順はガウスの消去法とほとんど同じだが, 行列 $b$ の値だけを変えて繰り返し実行する場合,
ガウスの消去法より高速に計算を行える.

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 "LU"
disp_matrix a

'Ly=b から y を求める (前進代入)
Private y: y = Array(0.0,0.0,0.0,0.0)
forward_substitution a,b,y

WScript.StdOut.WriteLine "Y"
disp_vector y

'Ux=y から x を求める (後退代入)
Private x: x = Array(0.0,0.0,0.0,0.0)
backward_substitution a,y,x

WScript.StdOut.WriteLine "X"
disp_vector x

'前進消去
Private Sub forward_elimination(ByRef a, ByVal 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
            a(row)(pivot) = s                                 'これが 下三角行列
            'b(row) = b(row) - b(pivot) * s                   'この値は変更しない
        Next
    Next
End Sub
'前進代入
Private Sub forward_substitution(ByVal a, ByVal b, ByRef y)
    Dim row, col
    For row = 0 To (N - 1)
        For col = 0 To row
            b(row) = b(row) - a(row)(col) * y(col)
        Next
        y(row) = b(row)
    Next
End Sub
'後退代入
Private Sub backward_substitution(ByVal a, ByVal y, ByRef x)
    Dim row, col
    For row = (N - 1) To 0 Step -1
        For col = (N - 1) To (row + 1) Step - 1
            y(row) = y(row) - a(row)(col) * x(col)
        Next
        x(row) = y(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:\1005.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

LU
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.2222222222    7.5555555556   -2.2222222222    0.7777777778
 -0.1111111111   -0.2352941176    6.5882352941   -1.7058823529
  0.1111111111   -0.1617647059   -0.3750000000    5.3750000000
Y
 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()

// LU分解
// 前進消去
forward_elimination(a,b)

WScript.StdOut.WriteLine("LU")
disp_matrix(a)

// Ly=b から y を求める (前進代入)
var y = [0,0,0,0]
forward_substitution(a,b,y)

WScript.StdOut.WriteLine("Y")
disp_vector(y)

// Ux=y から x を求める (後退代入)
var x = [0,0,0,0]
backward_substitution(a,y,x)

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

// 前進消去
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 // これが 上三角行列
            a[row][pivot] = s                    // これが 下三角行列
            // b[row]    -= b[pivot] * s         // この値は変更しない
        }
    }
}
// 前進代入
function forward_substitution(a, b, y)
{
    for (row = 0; row < N; row++)
    {
        for (col = 0; col < row; col++)
            b[row] -= a[row][col] * y[col]
        y[row] = b[row]
    }
}
// 後退代入
function backward_substitution(a, y, x)
{
    for (row = N - 1; row >= 0; row--)
    {
        for (col = N - 1; col > row; col--)
            y[row] -= a[row][col] * x[col]
        x[row] = y[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:\1005.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

LU
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.2222222222    7.5555555556   -2.2222222222    0.7777777778
 -0.1111111111   -0.2352941176    6.5882352941   -1.7058823529
  0.1111111111   -0.1617647059   -0.3750000000    5.3750000000
Y
 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 # これが 上三角行列
            }
            $a[$row][$pivot]    = $s                       # これが 下三角行列
            # $b[$row]         -= $b[$pivot]          * $s # この値は変更しない
        }
    }
}
# 前進代入
function forward_substitution($a, $b, $y)
{
    foreach ($row in 0..$N)
    {
        foreach ($col in 0..$row)
        {
            $b[$row] -= $a[$row][$col] * $y[$col]
        }
        $y[$row] = $b[$row]
    }
}
# 後退代入
function backward_substitution($a, $y, $x)
{
    foreach ($row in $N..0)
    {
        if (($row+1) -le $N)
        {
            foreach ($col in $N..($row + 1))
            {
                $y[$row] -= $a[$row][$col] * $x[$col]
            }
        }
        $x[$row] = $y[$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 "LU"
disp_matrix $a

# Ly=b から y を求める (前進代入)
$y = 0.0, 0.0, 0.0, 0.0
forward_substitution $a $b $y

Write-Host "Y"
disp_vector $y

# Ux=y から x を求める (後退代入)
$x = 0.0, 0.0, 0.0, 0.0
backward_substitution $a $y $x

Write-Host "X"
disp_vector $x

Z:\>powershell -file Z:\1005.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

LU
  9.0000000000  2.0000000000  1.0000000000  1.0000000000
  0.2222222222  7.5555555556 -2.2222222222  0.7777777778
 -0.1111111111 -0.2352941176  6.5882352941 -1.7058823529
  0.1111111111 -0.1617647059 -0.3750000000  5.3750000000
Y
 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);

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 "LU\n";
disp_matrix(\@a);

# Ly=b から y を求める (前進代入)
my @y = (0,0,0,0);
forward_substitution(\@a,\@b,\@y);
print "Y\n";
disp_vector(\@y);

# Ux=y から x を求める (後退代入)
my @x = (0,0,0,0);
backward_substitution(\@a,\@y,\@x);
print "X\n";
disp_vector(\@x);

# 前進消去
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; # これが 上三角行列
            }
            $$a[$row][$pivot]  = $s;                       # これが 下三角行列
            # $$b[$row]       -= $$b[$pivot] * $s;         # この値は変更しない
        }
    }
}
# Ly=b から y を求める (前進代入)
sub forward_substitution
{
    my ($a, $b, $y) = @_;

    for $row (0..N)
    {
        for $col (0..$row)
        {
            $$b[$row] -= $$a[$row][$col] * $$y[$col];
        }
        $$y[$row] = $$b[$row];
    }
}
# Ux=y から x を求める (後退代入)
sub backward_substitution
{
    my ($a, $y, $x) = @_;

    for ($row = N; $row >= 0; $row--)
    {
        for ($col = N; $col > $row; $col--)
        {
            $$y[$row] -= $$a[$row][$col] * $$x[$col];
        }
        $$x[$row] = $$y[$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:\1005.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

LU
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.2222222222    7.5555555556   -2.2222222222    0.7777777778
 -0.1111111111   -0.2352941176    6.5882352941   -1.7058823529
  0.1111111111   -0.1617647059   -0.3750000000    5.3750000000
Y
 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 "LU\n";
disp_matrix($a);

# Ly=b から y を求める (前進代入)
$y = [0,0,0,0];
forward_substitution($a,$b,$y);
print "Y\n";
disp_vector($y);

# Ux=y から x を求める (後退代入)
$x = [0,0,0,0];
backward_substitution($a,$y,$x);
print "X\n";
disp_vector($x);

# 前進消去
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; # これが 上三角行列
            $a[$row][$pivot]    = $s;                       # これが 下三角行列
            # $b[$row]         -= $b[$pivot]          * $s; # この値は変更しない
        }
    }
}
# Ly=b から y を求める (前進代入)
function forward_substitution($a, $b, &$y)
{
    foreach (range(0, N) as $row)
    {
        foreach (range(0, $row) as $col)
            $b[$row] -= $a[$row][$col] * $y[$col];
        $y[$row] = $b[$row];
    }
}
# Ux=y から x を求める (後退代入)
function backward_substitution($a, $y, &$x)
{
    foreach (range(N, 0) as $row)
    {
        if (($row + 1) <= N)
            foreach (range(N, $row + 1) as $col)
                $y[$row] -= $a[$row][$col] * $x[$col];
        $x[$row] = $y[$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:\1005.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

LU
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.2222222222    7.5555555556   -2.2222222222    0.7777777778
 -0.1111111111   -0.2352941176    6.5882352941   -1.7058823529
  0.1111111111   -0.1617647059   -0.3750000000    5.3750000000
Y
 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 # これが 上三角行列
            a[row][pivot]    = s                    # これが 下三角行列
            # b[row]        -= b[pivot]         * s # この値は変更しない

# 前進代入
def forward_substitution(a, b, y):
    for row in range(0, N, 1):
        for col in range(0, row, 1):
            b[row] -= a[row][col] * y[col]

        y[row] = b[row]

# 後退代入
def backward_substitution(a, y, x):
    for row in range(N - 1, -1, -1):
        for col in range(N - 1, row, -1):
            y[row] -= a[row][col] * x[col]
        x[row] = y[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]

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

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

# 前進消去
forward_elimination(a,b)

print "LU"
disp_matrix(a)

# Ly=b から y を求める (前進代入)
y = [0.0, 0.0, 0.0, 0.0]
forward_substitution(a,b,y)

print "Y"
disp_vector(y)

# Ux=y から x を求める (後退代入)
x =  [ 0.0,  0.0,  0.0,  0.0]
backward_substitution(a,y,x)

print "X"
disp_vector(x)

Z:\>python Z:\1005.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

LU
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.2222222222    7.5555555556   -2.2222222222    0.7777777778
 -0.1111111111   -0.2352941176    6.5882352941   -1.7058823529
  0.1111111111   -0.1617647059   -0.3750000000    5.3750000000
Y
 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
            a[row][pivot]    = s                 # これが 下三角行列
            # b[row]        -= b[pivot]      * s # この値は変更しない
        end
    end
end
# Ly=b から y を求める (前進代入)
def forward_substitution(a, b, y)
    (0..N).each do |row|
        (0..row).each do |col|
            b[row] -= a[row][col] * y[col]
        end
        y[row] = b[row]
    end
end
# Ux=y から x を求める (後退代入)
def backward_substitution(a, y, x)
    N.downto(0) do |row|
        N.downto(row + 1) do |col|
            y[row] -= a[row][col] * x[col]
        end
        x[row] = y[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 "LU"
disp_matrix(a)

# Ly=b から y を求める (前進代入)
y = [0.0,0.0,0.0,0.0]
forward_substitution(a,b,y)

puts "Y"
disp_vector(y)

# Ux=y から x を求める (後退代入)
x =  [0.0,0.0,0.0,0.0]
backward_substitution(a,y,x)

puts "X"
disp_vector(x)

Z:\>ruby Z:\1005.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

LU
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.2222222222    7.5555555556   -2.2222222222    0.7777777778
 -0.1111111111   -0.2352941176    6.5882352941   -1.7058823529
  0.1111111111   -0.1617647059   -0.3750000000    5.3750000000
Y
 20.0000000000   11.5555555556   12.9411764706   21.5000000000
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

Groovy

N = 3

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]]  as double[][]
def b =  [ 8.0, 17.0, 20.0, 16.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("LU")
disp_matrix(a)

// Ly=b から y を求める (前進代入)
y = [0.0,0.0,0.0,0.0]
forward_substitution(a,b,y)

println "Y"
disp_vector(y)

// Ux=y から x を求める (後退代入)
x =  [0.0,0.0,0.0,0.0]
backward_substitution(a,y,x)

println("X")
disp_vector(x)

// ピボット選択
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 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 // これが 上三角行列
            a[row][pivot]    = s                 // これが 下三角行列
            // b[row]       -= b[pivot]      * s // この値は変更しない
        }
    }
}
// Ly=b から y を求める (前進代入)
def forward_substitution(a, b, y) {
    for (row in 0..N) {
        for (col in 0..row)
            b[row] -= a[row][col] * y[col]
        y[row] = b[row]
    }
}
// Ux=y から x を求める (後退代入)
def backward_substitution(a, y, x) {
    for (row = N; row >= 0; row--) {
        for (col = N; col > row; col--)
            y[row] -= a[row][col] * x[col]
        x[row] = y[row] / a[row][row]
    }
}
// １次元配列を表示
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 Groovy1005.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

LU
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.2222222222    7.5555555556   -2.2222222222    0.7777777778
 -0.1111111111   -0.2352941176    6.5882352941   -1.7058823529
  0.1111111111   -0.1617647059   -0.3750000000    5.3750000000
Y
 20.0000000000   11.5555555556   12.9411764706   21.5000000000
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

Pascal

program Pas1005(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; // これが 上三角行列
            a[row,pivot]   := s;                             // これが 下三角行列
            // b[row]      := b[row]     - b[pivot]     * s; // この値は変更しない
        end;
    end;
end;
// 前進代入
procedure forward_substitution(a:TwoDimArray; b:array of Double; var y:array of Double);
var
    row, col:Integer;
begin
    for row := Low(a) to High(a) do
    begin
        for col := Low(a) to row do
            b[row] := b[row] - a[row,col] * y[col];
        y[row] := b[row];
    end;
end;
// 後退代入
procedure backward_substitution(a:TwoDimArray; y:array of Double; var x: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
            y[row] := y[row] - a[row,col] * x[col];
        x[row] := y[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);
    x :array [0..N] of Double = (0.0, 0.0, 0.0, 0.0);
    y :array [0..N] of Double = (0.0, 0.0, 0.0, 0.0);
begin
    // ピボット選択
    pivoting(a,b);

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

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

    writeln('LU');
    disp_matrix(a);

    // Ly=b から y を求める (前進代入)
    forward_substitution(a,b,y);

    writeln('Y');
    disp_vector(y);

    // Ux=y から x を求める (後退代入)
    backward_substitution(a,y,x);

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

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

Z:\>Pas1005
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

LU
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.2222222222    7.5555555556   -2.2222222222    0.7777777778
 -0.1111111111   -0.2352941176    6.5882352941   -1.7058823529
  0.1111111111   -0.1617647059   -0.3750000000    5.3750000000
Y
 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 Ada1005 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);
    y:Long_Float_Array := (0.0, 0.0, 0.0, 0.0);
    x:Long_Float_Array := (0.0, 0.0, 0.0, 0.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;
                a(row,pivot)   := s;                             -- これが 下三角行列
                -- b(row)      := b(row)     - b(pivot)     * s; -- この値は変更しない
            end loop;
        end loop;
    end forward_elimination;

    -- 前進代入
    procedure forward_substitution(a:Long_Float_TwoDimArray; b:in out Long_Float_Array; y:in out Long_Float_Array) is
    begin
        for row in a'Range loop
            for col in a'First..row loop
                b(row) := b(row) - a(row,col) * y(col);
            end loop;
            y(row) := b(row);
        end loop;
    end forward_substitution;

    -- 後退代入
    procedure backward_substitution(a:Long_Float_TwoDimArray; y:in out Long_Float_Array; x:in out Long_Float_Array) is
    begin
        for row in reverse a'Range loop
            for col in reverse (row + 1)..a'Last loop
                y(row) := y(row) - a(row,col) * x(col);
            end loop;
            x(row) := y(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("LU");
    disp_matrix(a);

    -- Ly=b から y を求める (前進代入)
    forward_substitution(a,b,y);

    Put_Line("Y");
    disp_vector(y);

    -- Ux=y から x を求める (後退代入)
    backward_substitution(a,y,x);

    Put_Line("X");
    disp_vector(x);
end Ada1005;

xxxxxx@yyyyyy /Z
$ gnatmake Ada1005.adb

xxxxxx@yyyyyy /Z
$ Ada1005
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

LU
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.2222222222    7.5555555556   -2.2222222222    0.7777777778
 -0.1111111111   -0.2352941176    6.5882352941   -1.7058823529
  0.1111111111   -0.1617647059   -0.3750000000    5.3750000000
Y
 20.0000000000   11.5555555556   12.9411764706   21.5000000000
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

VB.NET

Option Explicit

Module VB1005
    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("LU")
        disp_matrix(a)

        'Ly=b から y を求める (前進代入)
        Dim y()  As Double = {0,0,0,0}
        forward_substitution(a,b,y)

        Console.WriteLine("Y")
        disp_vector(y)

        'Ux=y から x を求める (後退代入)
        Dim x()  As Double = {0,0,0,0}
        backward_substitution(a,y,x)

        Console.WriteLine("X")
        disp_vector(x)
    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
                a(row,pivot) = s                   'これが 下三角行列
                'b(row) -= b(pivot) * s            'この値は変更しない
            Next
        Next
    End Sub

    '前進代入
    Private Sub forward_substitution(ByVal a(,) As Double, ByVal b() As Double, ByRef y() As Double)
        For row As Integer = 0 To N
            For col As Integer = 0 To row
                b(row) -= a(row,col) * y(col)
            Next
            y(row) = b(row)
        Next
    End Sub

    '後退代入
    Private Sub backward_substitution(ByVal a(,) As Double, ByVal y() As Double, ByRef x() As Double)
        For row As Integer = N To 0 Step -1
            For col As Integer = N To (row + 1) Step - 1
                y(row) -= a(row,col) * x(col)
            Next
            x(row) = y(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 VB1005.vb

Z:\>VB1005
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

LU
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.2222222222    7.5555555556   -2.2222222222    0.7777777778
 -0.1111111111   -0.2352941176    6.5882352941   -1.7058823529
  0.1111111111   -0.1617647059   -0.3750000000    5.3750000000
Y
 20.0000000000   11.5555555556   12.9411764706   21.5000000000
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

C#

using System;

public class CS1005
{
    private const int N = 4;

    // LU分解
    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("LU");
        disp_matrix(a);

        // Ly=b から y を求める (前進代入)
        double[] y = {0,0,0,0};
        forward_substitution(a,b,y);

        Console.WriteLine("Y");
        disp_vector(y);

        // Ux=y から x を求める (後退代入)
        double[] x = {0,0,0,0};
        backward_substitution(a,y,x);

        Console.WriteLine("X");
        disp_vector(x);
    }
    // 前進消去
    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; // これが 上三角行列
                a[row,pivot] = s;                   // これが 下三角行列
                // b[row]   -= b[pivot] * s;        // この値は変更しない
            }
        }
    }
    // 前進代入
    private static void forward_substitution(double[,] a, double[] b, double[] y)
    {
        for (int row = 0; row < N; row++)
        {
            for (int col = 0; col < row; col++)
                b[row] -= a[row,col] * y[col];
            y[row] = b[row];
        }
    }
    // 後退代入
    private static void backward_substitution(double[,] a, double[] y, double[] x)
    {
        for (int row = N - 1; row >= 0; row--)
        {
            for (int col = N - 1; col > row; col--)
                y[row] -= a[row,col] * x[col];
            x[row] = y[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 CS1005.cs

Z:\>CS1005
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

LU
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.2222222222    7.5555555556   -2.2222222222    0.7777777778
 -0.1111111111   -0.2352941176    6.5882352941   -1.7058823529
  0.1111111111   -0.1617647059   -0.3750000000    5.3750000000
Y
 20.0000000000   11.5555555556   12.9411764706   21.5000000000
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

Java

import java.lang.*;

public class Java1005 {

    private static final int N = 4;

    // LU分解
    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("LU");
        disp_matrix(a);

        // Ly=b から y を求める (前進代入)
        double y[] = {0,0,0,0};
        forward_substitution(a,b,y);

        System.out.println("Y");
        disp_vector(y);

        // Ux=y から x を求める (後退代入)
        double x[] = {0,0,0,0};
        backward_substitution(a,y,x);

        System.out.println("X");
        disp_vector(x);
    }
    // 前進消去
    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; // これが 上三角行列
                a[row][pivot] = s;                    // これが 下三角行列
                // b[row]    -= b[pivot] * s;         // この値は変更しない
            }
        }
    }
    // 前進代入
    private static void forward_substitution(double[][] a, double[] b, double[] y) {
        for (int row = 0; row < N; row++) {
            for (int col = 0; col < row; col++)
                b[row] -= a[row][col] * y[col];
            y[row] = b[row];
        }
    }
    // 後退代入
    private static void backward_substitution(double[][] a, double[] y, double[] x) {
        for (int row = N - 1; row >= 0; row--) {
            for (int col = N - 1; col > row; col--)
                y[row] -= a[row][col] * x[col];
            x[row] = y[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 Java1005.java

Z:\>java Java1005
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

LU
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.2222222222    7.5555555556   -2.2222222222    0.7777777778
 -0.1111111111   -0.2352941176    6.5882352941   -1.7058823529
  0.1111111111   -0.1617647059   -0.3750000000    5.3750000000
Y
 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]);
// 前進代入
void forward_substitution(double a[N][N], double b[N], double y[N]);
// 後退代入
void backward_substitution(double a[N][N], double y[N], double x[N]);
// １次元配列を表示
void disp_vector(double row[N]);
// ２次元配列を表示
void disp_matrix(double matrix[N][N]);

// LU分解
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);

    cout << "LU" << endl;
    disp_matrix(a);

    // Ly=b から y を求める (前進代入)
    double y[N] = {0,0,0,0};
    forward_substitution(a, b, y);

    cout << "Y" << endl;
    disp_vector(y);

    // Ux=y から x を求める (後退代入)
    double x[N] = {0,0,0,0};
    backward_substitution(a, y, x);

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

    return 0;
}
// 前進消去
void forward_elimination(double a[N][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; // これが 上三角行列
            a[row][pivot] = s;                    // これが 下三角行列
            // b[row]    -= b[pivot] * s;         // この値は変更しない
        }
    }
}
// 前進代入
void forward_substitution(double a[N][N], double b[N], double y[N])
{
    for (int row = 0; row < N; row++)
    {
        for (int col = 0; col < row; col++)
            b[row] -= a[row][col] * y[col];
        y[row] = b[row];
    }
}
// 後退代入
void backward_substitution(double a[N][N], double y[N], double x[N])
{
    for (int row = N - 1; row >= 0; row--)
    {
        for (int col = N - 1; col > row; col--)
            y[row] -= a[row][col] * x[col];
        x[row] = y[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 CP1005.cpp
cp1005.cpp:

Z:\>CP1005
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

LU
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.2222222222    7.5555555556   -2.2222222222    0.7777777778
 -0.1111111111   -0.2352941176    6.5882352941   -1.7058823529
  0.1111111111   -0.1617647059   -0.3750000000    5.3750000000
Y
 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 forward_substitution(double a[N][N], double b[N], double y[N]);
// 後退代入
void backward_substitution(double a[N][N], double y[N], double x[N]);
// １次元配列を表示
void disp_vector(double row[N]);
// ２次元配列を表示
void disp_matrix(double matrix[N][N]);

// LU分解
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("LU\n");
    disp_matrix(a);

    // Ly=b から y を求める (前進代入)
    double y[4] = {0,0,0,0};
    forward_substitution(a,b,y);

    printf("Y\n");
    disp_vector(y);

    // Ux=y から x を求める (後退代入)
    double x[4] = {0,0,0,0};
    backward_substitution(a,y,x);

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

    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; // これが 上三角行列
            a[row][pivot] = s;                    // これが 下三角行列
            // b[row]    -= b[pivot] * s;         // この値は変更しない
        }
    }
}
// 前進代入
void forward_substitution(double a[N][N], double b[N], double y[N])
{
    int row, col;
    for (row = 0; row < N; row++)
    {
        for (col = 0; col < row; col++)
            b[row] -= a[row][col] * y[col];
        y[row] = b[row];
    }
}
// 後退代入
void backward_substitution(double a[N][N], double y[N], double x[N])
{
    int row, col;
    for (row = N - 1; row >= 0; row--)
    {
        for (col = N - 1; col > row; col--)
            y[row] -= a[row][col] * x[col];
        x[row] = y[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_progress(double a[N][N], double b[N])
{
    int i, j;
    for (i = 0; i < N; i++)
    {
        for (j = 0; j < N; j++)
            printf("%14.10f\t", a[i][j]);
        printf("%14.10f\n", b[i]);
    }
    printf("\n");
}
// １次元配列を表示
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 OC1005 OC1005.m -lobjc -lgnustep-base -I $INCLUDE -L $LIB $CFLAGS

xxxxxx@yyyyyy /Z
$ OC1005
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

LU
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.2222222222    7.5555555556   -2.2222222222    0.7777777778
 -0.1111111111   -0.2352941176    6.5882352941   -1.7058823529
  0.1111111111   -0.1617647059   -0.3750000000    5.3750000000
Y
 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;

// LU分解
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("LU");
    disp_matrix(a);

    // Ly=b から y を求める (前進代入)
    double y[N] = [0,0,0,0];
    forward_substitution(a,b,y);

    writefln("Y");
    disp_vector(y);

    // Ux=y から x を求める (後退代入)
    double x[N] = [0,0,0,0];
    backward_substitution(a,y,x);

    writefln("X");
    disp_vector(x);
}
// 前進消去
void forward_elimination(ref double[N][N] a, 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; // これが 上三角行列
            a[row][pivot] = s;                    // これが 下三角行列
            // b[row]    -= b[pivot] * s;         // この値は変更しない
        }
    }
}
// 前進代入
void forward_substitution(double[N][N] a, double[N] b, ref double[N] y)
{
    foreach (row; 0..N)
    {
        foreach (col; 0..row)
            b[row] -= a[row][col] * y[col];
        y[row] = b[row];
    }
}
// 後退代入
void backward_substitution(double[N][N] a, double[N] y, ref double[N] x)
{
    foreach_reverse (row; 0..N)
    {
        foreach_reverse (col; (row+1)..N)
            y[row] -= a[row][col] * x[col];
        x[row] = y[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 D1005.d

Z:\>D1005
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

LU
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.2222222222    7.5555555556   -2.2222222222    0.7777777778
 -0.1111111111   -0.2352941176    6.5882352941   -1.7058823529
  0.1111111111   -0.1617647059   -0.3750000000    5.3750000000
Y
 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

// LU分解
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("LU")
    disp_matrix(a)

    // Ly=b から y を求める (前進代入)
    var y []float64 = []float64{0,0,0,0}
    forward_substitution(a,b,y)

    fmt.Println("Y")
    disp_vector(y)

    // Ux=y から x を求める (後退代入)
    var x []float64 = []float64{0,0,0,0}
    backward_substitution(a,y,x)

    fmt.Println("X")
    disp_vector(x)
}
// 前進消去
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  // これが 上三角行列
            }
            a[row][pivot] = s                     // これが 下三角行列
            // b[row]    -= b[pivot] * s          // この値は変更しない
        }
    }
}
// 前進代入
func forward_substitution(a [N][N]float64, b []float64, y []float64) {
    for row := 0; row < N; row++ {
        for col := 0; col < row; col++ {
            b[row] -= a[row][col] * y[col]
        }
        y[row] = b[row]
    }
}
// 後退代入
func backward_substitution(a [N][N]float64, y []float64, x []float64) {
    for row := N - 1; row >= 0; row-- {
        for col := N - 1; col > row; col-- {
            y[row] -= a[row][col] * x[col]
        }
        x[row] = y[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 GO1005.go

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

Z:\>GO1005
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

LU
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.2222222222    7.5555555556   -2.2222222222    0.7777777778
 -0.1111111111   -0.2352941176    6.5882352941   -1.7058823529
  0.1111111111   -0.1617647059   -0.3750000000    5.3750000000
Y
 20.0000000000   11.5555555556   12.9411764706   21.5000000000
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

Scala

object Scala1005 {
    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("LU")
        disp_matrix(a)

        // Ly=b から y を求める (前進代入)
        var y:Array[Double] = Array(0,0,0,0)
        forward_substitution(a,b,y)

        println("Y")
        disp_vector(y)

        // Ux=y から x を求める (後退代入)
        var x:Array[Double] = Array(0,0,0,0)
        backward_substitution(a,y,x)

        println("X")
        disp_vector(x)
    }
    // 前進消去
    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 // これが 上三角行列
                a(row)(pivot) = s                       // これが 下三角行列
                // b(row)     -= b(pivot)           * s // この値は変更しない
            }
        }
    }
    // 前進代入
    def forward_substitution(a:Array[Array[Double]], b:Array[Double], y:Array[Double]) {
        for (row <- 0 to N) {
            for (col <- 0 to row - 1)
                b(row) -= a(row)(col) * y(col)
            y(row) = b(row)
        }
    }
    // 後退代入
    def backward_substitution(a:Array[Array[Double]], y:Array[Double], x:Array[Double]) {
        for (row <- (0 to N).reverse) {
            for (col <- (row + 1 to N).reverse)
                y(row) -= a(row)(col) * x(col)
            x(row) = y(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 Scala1005.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

LU
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.2222222222    7.5555555556   -2.2222222222    0.7777777778
 -0.1111111111   -0.2352941176    6.5882352941   -1.7058823529
  0.1111111111   -0.1617647059   -0.3750000000    5.3750000000
Y
 20.0000000000   11.5555555556   12.9411764706   21.5000000000
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP

F#

module Fs1005

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 // これが 上三角行列
            a.[row].[pivot]   <- s                                      // これが 下三角行列
            // b.[row]        <- b.[row]       - b.[pivot]          * s // この値は変更しない

// 前進代入
let forward_substitution (a:float[][]) (b:float[]) (y:float[]) =
    for row in [0..N] do
        for col in 0..(row - 1) do
            b.[row] <- b.[row] - a.[row].[col] * y.[col]
        y.[row] <- b.[row]

// 後退代入
let backward_substitution (a:float[][]) (y:float[]) (x:float[]) =
    for row in (List.rev [0..N]) do
        for col in (List.rev [row + 1 .. N]) do
            y.[row] <- y.[row] - a.[row].[col] * x.[col]
        x.[row] <- y.[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 "LU"
disp_matrix a

// Ly=b から y を求める (前進代入)
let mutable y:float[] = [| 0.0;0.0;0.0;0.0 |]
forward_substitution a b y

printfn "Y"
disp_vector y

// Ux=y から x を求める (後退代入)
let mutable x:float[] = [| 0.0;0.0;0.0;0.0 |]
backward_substitution a y x

printfn "X"
disp_vector x

exit 0

Z:\>fsi  --nologo --quiet Fs1005.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

LU
  9.0000000000  2.0000000000  1.0000000000  1.0000000000
  0.2222222222  7.5555555556 -2.2222222222  0.7777777778
 -0.1111111111 -0.2352941176  6.5882352941 -1.7058823529
  0.1111111111 -0.1617647059 -0.3750000000  5.3750000000
Y
 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 (inc pivot) (nth a pivot)) (drop (inc pivot) (nth a row))))
    (def a1 (concat (take pivot (nth a row)) (cons s        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 b)
        (vector a2 b)))

(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 forward_substitution [row a b]
    (def a_map (map
        (fn [a_col b_col] (* a_col b_col))
            (take row (nth a row)) b))
    (def s (apply + a_map))
    (def y (- (nth b row) s))
    (def b2 (concat (take row b) (cons y (drop (inc row) b))))
    (if (< row (dec N))
        (cons y (forward_substitution (inc row) a b2))
        (cons y [])))

;後退代入
(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 "LU")
(disp_matrix a2)

;Ly=b から y を求める (前進代入)
(def y (forward_substitution 0 a2 b2))

(println "Y")
(disp_vector (reverse y))

;Ux=y から x を求める (後退代入)
(def x (backward_substitution (dec N) a2 y))

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

Z:\>java -cp C:\ProgramFiles\clojure-1.5.1\clojure-1.5.1.jar clojure.main Clj1005.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

LU
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.2222222222    7.5555555556   -2.2222222222    0.7777777778
 -0.1111111111   -0.2352941176    6.5882352941   -1.7058823529
  0.1111111111   -0.1617647059   -0.3750000000    5.3750000000
Y
 21.5000000000   12.9411764706   11.5555555556   20.0000000000
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+1) (a!!pivot)) (drop (pivot+1) (a!!row))
        a_map   = map (\(a_pivot, a_row) -> a_row - a_pivot * s) a_zip
        new_row = (take pivot (a!!row)) ++ (s:a_map                   )
        a2      = (take row   a       ) ++ (new_row:(drop (row+1)   a))
    in
        if row < (n - 1)
            then
                forward_elim_loop pivot (row+1) a2 b
            else
                (a2, b)

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)

-- 前進代入
forward_substitution::Int->[[Double]]->[Double]->[Double]
forward_substitution row a b =
    let
        a_zip = zip (take row (a!!row)) b
        s  = sum $ map(\(a_col, b_col) -> a_col * b_col) a_zip
        y  = b!!row - s
        b2 = (take row b) ++ (y:(drop (row+1) b))
    in
        if row < (n - 1)
            then y:(forward_substitution (row+1) a b2)
            else y:[]

-- 後退代入
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 "LU"
    disp_matrix a2

    -- Ly=b から y を求める (前進代入)
    let y = forward_substitution 0 a2 b2

    putStrLn "Y"
    disp_vector y

    -- Ux=y から x を求める (後退代入)
    let x = backward_substitution (n-1) a2 y

    putStrLn "X"
    disp_vector (reverse x)

Z:\>runghc Hs1005.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

LU
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  0.2222222222    7.5555555556   -2.2222222222    0.7777777778
 -0.1111111111   -0.2352941176    6.5882352941   -1.7058823529
  0.1111111111   -0.1617647059   -0.3750000000    5.3750000000
Y
 20.0000000000   11.5555555556   12.9411764706   21.5000000000
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000
TOP
inserted by FC2 system