File Coverage

lib/Math/MatrixDecomposition/LU.pm

Criterion	Covered	Total	%
statement	106	125	84.8
branch	20	38	52.6
condition	10	28	35.7
subroutine	11	13	84.6
pod	5	5	100.0
total	152	209	72.7

## Math/MatrixDecomposition/LU.pm --- LU decomposition.

# Copyright (C) 2010 Ralph Schleicher.  All rights reserved.

# This program is free software; you can redistribute it and/or
# modify it under the same terms as Perl itself.

## Code:

package Math::MatrixDecomposition::LU;

use strict;
use warnings;
use Carp;
use Exporter qw(import);
use Scalar::Util qw(looks_like_number);

use Math::BLAS;
use Math::MatrixDecomposition::Util qw(mod min);

BEGIN
{
  our $VERSION = '1.05';
  our @EXPORT_OK = qw(lu);
}

# Create a LU decomposition (convenience function).
sub lu
{
  __PACKAGE__->new (@_);
}

# Standard constructor.
sub new
{
  my $class = shift;
  my $self =
    {
     # LU matrix in row major layout.
     LU => [],

     # Number of matrix rows and columns.
     rows => 0,
     columns => 0,

     # Pivot row indices (row permutations).
     pivot => [],

     # Sign of determinant.  A value of zero means that
     # the matrix is singular.
     sign => 0,
    };

  bless ($self, ref ($class) || $class);

  # Process arguments.
  $self->decompose (@_)
    if @_ > 0;

  # Return object.
  $self;
}

# LU decomposition.
sub decompose
{
  my $self = shift;

  # Check arguments.
  my $a = shift;

  my $m = @_ > 0 && looks_like_number ($_[0]) ? shift : sqrt (@$a);
  my $n = @_ > 0 && looks_like_number ($_[0]) ? shift : $m ? @$a / $m : 0;

  croak ('Invalid argument')
    if (@$a != ($m * $n)
	|| mod ($m, 1) != 0 || $m < 1
	|| mod ($n, 1) != 0 || $n < 1);

  # Get properties.
  my %prop = @_;

  $prop{capture} //= 0;

  # Index of last row/column.
  my $last_row = $m - 1;
  my $last_col = $n - 1;

  # Matrix elements.
  if ($prop{capture})
    {
      # Work in-place.
      $$self{LU} = $a;
    }
  else
    {
      # Copy matrix elements.
      $$self{LU} //= [];
      @{$$self{LU}} = @$a;
      $a = $$self{LU};
    }

  # Matrix size.
  $$self{rows} = $m;
  $$self{columns} = $n;

  # Pivot row indices (row permutations).
  my $pivot = $$self{pivot};

  @$pivot = ();

  # Sign of determinant.
  my $sign = 1;

  # Work variables.
  my ($i, $j, $k, $p, $row, $tem);

  for $k (0 .. min ($last_row, $last_col))
    {
      # Search pivot element.
      $p = $k + (blas_amax_val ($m - $k, $a,
				x_ind => $k * $n + $k,
				x_incr => $n))[0];

      push (@$pivot, $p);

      # Swap rows.
      if ($p != $k)
	{
	  blas_swap ($n, $a, $a,
		     x_ind => $k * $n,
		     y_ind => $p * $n);

	  $sign = - $sign;
	}

      # Divide by the pivot element.
      if ($$a[$k * $n + $k] == 0)
	{
	  $sign = 0;
	  next;
	}

      for $i ($k + 1 .. $last_row)
	{
	  $$a[$i * $n + $k] /= $$a[$k * $n + $k];

	  blas_axpby ($last_col - $k, $a, $a,
		      alpha => - $$a[$i * $n + $k],
		      x_ind => $k * $n + $k + 1,
		      y_ind => $i * $n + $k + 1);
	}
    }

  # Save sign of determinant.
  $$self{sign} = $sign;

  # Return object.
  $self;
}

# Solve a system of linear equations.
sub solve
{
  my $self = shift;

  my $b = shift;
  my $x = shift;

  # LU matrix elements.
  my $a = $$self{LU};

  # Number of matrix rows and columns.
  my $m = $$self{rows};
  my $n = $$self{columns};

  croak ('Invalid argument')
    if @$b == 0 || @$b % $m != 0;

  my $l = @$b / $m;

  # Copy right-hand side.
  if (defined ($x))
    {
      @$x = @$b;
    }
  else
    {
      # Work in-place.
      $x = $b;
    }

  # Index of last row/column.
  my $last_row = $m - 1;
  my $last_col = $n - 1;
  my $last_vec = $l - 1;

  # Pivot row indices.
  my $pivot = $$self{pivot};

  # Work variables.
  my ($i, $j, $k, $p, $row, $tem);

  if ($l == 1)
    {
      # Apply row permutations.
      for $i (0 .. $#$pivot)
	{
	  $p = $$pivot[$i];
	  next if $p == $i;

	  @$x[$i, $p] = @$x[$p, $i];
	}

      # Forward substitution, that is solve 'Ly = b' for 'y'.
      # Elements of 'y' are saved in 'x'.
      for $i (1 .. $last_row)
	{
	  $row = $i * $n;

	  for $j (0 .. $i - 1)
	    {
	      $$x[$i] -= $$x[$j] * $$a[$row + $j];
	    }
	}

      # Backward substitution, that is solve 'Ux = y' for 'x'.
      for $i (reverse (0 .. $last_row))
	{
	  $row = $i * $n;

	  if ($$a[$row + $i] == 0)
	    {
	      # Matrix A is singular.
	      return unless $$x[$i] == 0;
	    }
	  else
	    {
	      for $j ($i + 1 .. $last_col)
		{
		  $$x[$i] -= $$x[$j] * $$a[$row + $j];
		}

	      $$x[$i] /= $$a[$row + $i];
	    }
	}
    }
  else
    {
      # Matrix A is singular.
      return if $$self{sign} == 0;

      # Apply row permutations.
      for $i (0 .. $#$pivot)
	{
	  $p = $$pivot[$i];
	  next if $p == $i;

	  blas_swap ($l, $x, $x,
		     x_ind => $i * $l,
		     y_ind => $p * $l);
	}

      # Forward substitution.
      for $k (0 .. $last_col)
	{
	  for $i ($k + 1 .. $last_col)
	    {
	      blas_axpby ($l, $x, $x,
			  alpha => - $$a[$i * $n + $k],
			  x_ind => $k * $l,
			  y_ind => $i * $l);
	    }
	}

      # Backward substitution.
      for $k (reverse (0 .. $last_col))
	{
	  blas_rscale ($l, $x,
		       alpha => $$a[$k * $n + $k],
		       x_ind => $k * $l);

	  for $i (0 .. $k - 1)
	    {
	      blas_axpby ($l, $x, $x,
			  alpha => - $$a[$i * $n + $k],
			  x_ind => $k * $l,
			  y_ind => $i * $l);
	    }
	}
    }

  # Resize result.
  splice (@$x, $n * $l)
    if $m > $n;

  # Fix rounding errors.
  for (@$x)
    {
      $_ = 0 + "$_";
    }

  # Return value.
  $x;
}

# Determinant.
sub det
{
  my $self = shift;

  my $m = $$self{rows};
  my $n = $$self{columns};

  croak ('Matrix has to be square')
    if $m != $n;

  # Calculate determinant.
  my $det = $$self{sign};

  if ($det)
    {
      my $u = $$self{LU};

      my $ind = 0;
      my $incr = $n + 1;

      for (1 .. $n)
	{
	  $det *= $$u[$ind];
	  $ind += $incr;
	}
    }

  $det;
}

1;

__END__