Source

perl-PDL-Lib-Linear-Solve / MatInv.pm

Full commit
package MatInv;

use strict;

use PDL;

require Exporter;

use vars qw(@EXPORT @ISA);

@EXPORT=qw(linear_solver);

@ISA=qw(Exporter);

my $eps = 0.00001;

sub linear_solver
{
    my $coeffs = shift;
    my $values = shift;

    my @dims = dims($coeffs);
    if ((scalar(@dims) != 2) || ($dims[0] != $dims[1]))
    {
        die "First argument is not a matrix";
    }
    
    my @v_dims = dims($values);
    if (@v_dims == 1)
    {
        $values = transpose($values);
        @v_dims = dims($values);
    }
    if ($v_dims[1] != $dims[1])
    {
        die "Second argument does not match in length";
    }

    $coeffs = $coeffs->copy();

    my $mr   = $dims[0];
    my $mc   = $mr;
    my $f;
    my $try;

    ROW: 
    for(my $i = 0; $i < $mr; $i++) {
        $try=$i;
        # make diagonal element nonzero if possible
        while(abs($coeffs->at($i,$try)) < $eps)
        {
            $try++;
            if ($try >= $mr)
            {
                last ROW;
            }
        }
        my ($temp, $mat);

        foreach $mat ($coeffs, $values)
        {
            $temp = $mat->slice(":,$i")->copy();
            $mat->slice(":,$i") .= $mat->slice(":,$try");
            $mat->slice(":,$try") .= $temp;
        }

        $f = (1 / ($coeffs->at($i,$i)));
        foreach $mat ($coeffs, $values)
        {
            $mat->slice(":,$i") *= $f;
        }
        
        # subtract multiple of designated row from other rows

        my @matrixes = 
            (
                [$coeffs, $coeffs->slice(":,$i")], 
                [$values, $values->slice(":,$i")]
            );
            
        for(my $j = 0; $j < $mr; $j++) {
            next if $i == $j;
            $f = $coeffs->at($i,$j);
            foreach $mat (@matrixes)
            {
                $mat->[0]->slice(":,$j") -= $mat->[1] * $f;
            }
        }
    }

    return ($coeffs, $values);
}

1;