gauss-jordan法

#!/usr/bin/perl
use strict;
use warnings;
use Data::Dumper;

my $a = [
    [1, -1, 1],
    [2, 1, -3],
    [3, 2, -1],
];

my $b = [
    [-5],
    [19],
    [16],
];

gauss_jordan($a, $b);

print Dumper($a);
print Dumper($b);

sub gauss_jordan {
    my ($a, $b) = @_;

    my $n = scalar @$a;
    my $m = scalar @{$b->[0]};

    my @check;
    my @switch;
    for my $i ( 0.. $n - 1 ) {

        # まだ使用していない行の最大値を取得する
        my $max_row;
        my $max_col;

        for my $j ( 0.. $n - 1 ) { # 行

            if ( $check[$j] ) { next };
            for my $k ( 0.. $n - 1 ) { # 列

                if ( $check[$k] ) { next };

                if ( !defined $max_row  || abs($a->[$max_row][$max_col]) < abs($a->[$j][$k]) ) {
                    $max_row = $j;
                    $max_col = $k;
                }
            }
        }

        $check[$max_col] = 1;

        if ( $max_row != $max_col ) {
            # 行入れ替え
            ($a->[$max_col], $a->[$max_row]) = ($a->[$max_row], $a->[$max_col]);
            ($b->[$max_col], $b->[$max_row]) = ($b->[$max_row], $b->[$max_col]);
            push @switch, [$max_col, $max_row];
        }

        # 対象行の処理
        my $pivnv = 1 / $a->[$max_col][$max_col];

        $a->[$max_col][$max_col] = 1;

        foreach ( @{$a->[$max_col]} ) {
            $_ *= $pivnv;
        }

        foreach ( @{$b->[$max_col]} ) {
            $_ *= $pivnv;
        }

        # 対象以外の行
        foreach my $l ( 0.. $n - 1 ) {
            next if ( $l == $max_col );

            my $dum = $a->[$l][$max_col]; # $a->[$l][$max_col] / 1

            $a->[$l][$max_col] = 0;

            foreach my $p ( 0.. $n - 1 ) {
                $a->[$l][$p] -= $a->[$max_col][$p] * $dum;
            }

            foreach my $p ( 0.. $m - 1 ) {
                $b->[$l][$p] -= $b->[$max_col][$p] * $dum;
            }
        }
    }
    for my $c (@switch) {
        foreach my $i (0.. $n - 1) {
            my $tmp = $a->[$i][$c->[0]];
            $a->[$i][$c->[0]] = $a->[$i][$c->[1]];
            $a->[$i][$c->[1]] = $tmp;
        }
    }
}

NUMERICAL RECIPED in Cのをperlで書き直してみた。 中学とかで習う連立方程式を解くやり方をそのままプログラムにした方法

bには結果が返り、 計算に必要なくなった列を単位行列として変換していくためaには逆行列が返る。