File Coverage

blib/lib/Math/MatrixReal.pm

Criterion	Covered	Total	%
statement	1370	1648	83.1
branch	453	726	62.4
condition	158	341	46.3
subroutine	156	167	93.4
pod	28	100	28.0
total	2165	2982	72.6

#  Copyright (c) 1996, 1997 by Steffen Beyer. All rights reserved.
#  Copyright (c) 1999 by Rodolphe Ortalo. All rights reserved.
#  Copyright (c) 2001-2015 by Jonathan Leto. All rights reserved.
#  This package is free software; you can redistribute it and/or
#  modify it under the same terms as Perl itself.

package Math::MatrixReal;

use strict;
use warnings;
use Carp;
use Data::Dumper;
use Scalar::Util qw/reftype/;
use vars qw(@ISA @EXPORT @EXPORT_OK %EXPORT_TAGS $VERSION);
require Exporter;

@ISA = qw(Exporter);
@EXPORT = qw();
@EXPORT_OK = qw(min max);
%EXPORT_TAGS = (all => [@EXPORT_OK]);
$VERSION = '2.12';

use overload
     '.'   => '_concat',
     'neg' => '_negate',
       '~' => '_transpose',
    'bool' => '_boolean',
       '!' => '_not_boolean',
     'abs' => '_norm',
       '+' => '_add',
       '-' => '_subtract',
       '*' => '_multiply',
       '/' => '_divide',
      '**' => '_exponent',
      '+=' => '_assign_add',
      '-=' => '_assign_subtract',
      '*=' => '_assign_multiply',
     '**=' => '_assign_exponent',
      '==' => '_equal',
      '!=' => '_not_equal',
       '<' => '_less_than',
      '<=' => '_less_than_or_equal',
       '>' => '_greater_than',
      '>=' => '_greater_than_or_equal',
      'eq' => '_equal',
      'ne' => '_not_equal',
      'lt' => '_less_than',
      'le' => '_less_than_or_equal',
      'gt' => '_greater_than',
      'ge' => '_greater_than_or_equal',
       '=' => '_clone',
      '""' => '_stringify',
'fallback' =>   undef;

=head1 NAME

Math::MatrixReal - Matrix of Reals

Implements the data type "matrix of real numbers" (and consequently also
"vector of real numbers").

=head1 SYNOPSIS

my $a = Math::MatrixReal->new_random(5, 5);

my $b = $a->new_random(10, 30, { symmetric=>1, bounded_by=>[-1,1] });

my $c = $b * $a ** 3;

my $d = $b->new_from_rows( [ [ 5, 3 ,4], [3, 4, 5], [ 2, 4, 1 ] ] );

print $a;

my $row        = ($a * $b)->row(3);

my $col        = (5*$c)->col(2);

my $transpose  = ~$c;

my $transpose  = $c->transpose;

my $inverse    = $a->inverse;

my $inverse    = 1/$a;

my $inverse    = $a ** -1;

my $determinant= $a->det;

=cut

sub new
{
    croak "Usage: \$new_matrix = Math::MatrixReal->new(\$rows,\$columns);" if (@_ != 3);

    my ($self,$rows,$cols) =  @_;
    my $class = ref($self) || $self || 'Math::MatrixReal';

    croak "Math::MatrixReal::new(): number of rows must be integer > 0"
      unless ($rows > 0 and $rows == int($rows) );

    croak "Math::MatrixReal::new(): number of columns must be integer > 0"
      unless ($cols > 0 and $cols == int($cols) );

    my $this = [ [ ], $rows, $cols ];

    # Create the first empty row and pre-lengthen
    my $empty = [ ];
    $#$empty = $cols - 1;

    map { $empty->[$_] = 0.0 } ( 0 .. $cols-1 );

    # Create a row at a time
    map { $this->[0][$_] = [ @$empty ] } ( 0 .. $rows-1);

    bless $this, $class;
}

sub new_diag {
    croak "Usage: \$new_matrix = Math::MatrixReal->new_diag( [ 1, 2, 3] );" unless (@_ == 2 );
    my ($self,$diag) = @_;
    my $n = scalar @$diag;

    croak "Math::MatrixReal::new_diag(): Third argument must be an arrayref" unless (ref($diag) eq "ARRAY");

    my $matrix = Math::MatrixReal->new($n,$n);

    map { $matrix->[0][$_][$_] = shift @$diag } ( 0 .. $n-1);
    return $matrix;
}
sub new_tridiag {
    croak "Usage: \$new_matrix = Math::MatrixReal->new_tridiag( [ 1, 2, 3], [ 4, 5, 6, 7], [-1,-2,-3] );" unless (@_ == 4 );
    my ($self,$lower,$diag,$upper) = @_;
    my $matrix;
    my ($l,$n,$m) =   (scalar(@$lower),scalar(@$diag),scalar(@$upper));
    my ($k,$p)=(-1,-1);

    croak "Math::MatrixReal::new_tridiag(): Arguments must be arrayrefs" unless
	ref $diag eq 'ARRAY' && ref $lower eq 'ARRAY' && ref $upper eq 'ARRAY';
    croak "Math::MatrixReal::new_tridiag(): new_tridiag(\$lower,\$diag,\$upper) diagonal dimensions incompatible" unless
	($l == $m && $n == ($l+1));

    $matrix = Math::MatrixReal->new_diag($diag);
    $matrix = $matrix->each(
		sub {
		    my ($e,$i,$j) = @_;
		    if    (($i-$j) == -1) { $k++; return $upper->[$k];}
		    elsif (    $i  == $j) {       return $e;          }
		    elsif (($i-$j) ==  1) { $p++; return $lower->[$p];}
		}
		);
    return $matrix;
}

sub new_random {
    croak "Usage: \$new_matrix = Math::MatrixReal->new_random(\$n,\$m, { symmetric => 1, bounded_by => [-5,5], integer => 1 } );"
        if (@_ < 2);
    my ($self, $rows, $cols, $options ) = @_;
    (($options = $cols) and ($cols = $rows)) if ref $cols eq 'HASH';
    my ($min,$max) = defined $options->{bounded_by} ?  @{ $options->{bounded_by} } : ( 0, 10);
    my $integer = $options->{integer};
    $self = ref($self) || $self || 'Math::MatrixReal';

    $cols ||= $rows;
    croak "Math::MatrixReal::new_random(): number of rows must = number of cols for symmetric option"
        if ($rows != $cols and $options->{symmetric} );

    croak "Math::MatrixReal::new_random(): number of rows must be integer > 0"
        unless ($rows > 0 and  $rows == int($rows) ) && ($cols > 0 and $cols == int($cols) ) ;

    croak "Math::MatrixReal::new_random(): bounded_by interval length must be > 0"
        unless (defined $min && defined $max && $min < $max );

    croak "Math::MatrixReal::new_random(): tridiag option only for square matrices"
        if (($options->{tridiag} || $options->{tridiagonal}) && $rows != $cols);

    croak "Math::MatrixReal::new_random(): diagonal option only for square matrices "
        if (($options->{diag} || $options->{diagonal}) && ($rows != $cols));

    my $matrix = Math::MatrixReal->new($rows,$cols);
    my $random_code = sub { $integer ? int($min + rand($max-$min)) : $min + rand($max-$min) } ;

    $matrix = $options->{diag} || $options->{diagonal} ? $matrix->each_diag($random_code) :  $matrix->each($random_code);
    $matrix = $matrix->each( sub {my($e,$i,$j)=@_; ( abs($i-$j)>1 ) ?  0 : $e } ) if ($options->{tridiag} || $options->{tridiagonal} );
    $options->{symmetric} ? 0.5*($matrix + ~$matrix) : $matrix;
}

sub new_from_string#{{{
{#{{{
    croak "Usage: \$new_matrix = Math::MatrixReal->new_from_string(\$string);"
      if (@_ != 2);

    my ($self,$string)  = @_;
    my $class  = ref($self) || $self || 'Math::MatrixReal';
    my ($line,$values);
    my ($rows,$cols);
    my ($row,$col);
    my ($warn,$this);

    $warn = $rows = $cols = 0;

    $values = [ ];
	while ($string =~ m!^\s* \[ \s+ ( (?: [+-]? \d+ (?: \. \d*)? (?: E [+-]? \d+ )? \s+ )+ ) \] \s*? \n !ix) {
			$line = $1; $string = $';
			$values->[$rows] = [ ]; @{$values->[$rows]} = split(' ', $line);
			$col = @{$values->[$rows]};
	 		if ($col != $cols) {
				unless ($cols == 0) { $warn = 1; }
				if ($col > $cols) { $cols = $col; }
			}
			$rows++;
	}
	if ($string !~ m/^\s*$/) {
        chomp $string;
        my $error_msg = "Math::MatrixReal::new_from_string(): syntax error in input string: $string";
        croak $error_msg;
    }
		if ($rows == 0) { croak "Math::MatrixReal::new_from_string(): empty input string"; }
		if ($warn) { warn "Math::MatrixReal::new_from_string(): missing elements will be set to zero!\n"; }
		$this = Math::MatrixReal::new($class,$rows,$cols);
		for ( $row = 0; $row < $rows; $row++ ) {
			for ( $col = 0; $col < @{$values->[$row]}; $col++ ) {
			$this->[0][$row][$col] = $values->[$row][$col];
			}
		}
		return $this;
}#}}}#}}}

# from Math::MatrixReal::Ext1 ([e-mail])
sub new_from_cols {
    my $self = shift;
    my $extra_args = ( @_ > 1 && ref($_[-1]) eq 'HASH' ) ? pop : {};
    $extra_args->{_type} = 'column';
    $self->_new_from_rows_or_cols(@_, $extra_args );
}
# from Math::MatrixReal::Ext1 ([e-mail])
sub new_from_columns {
    my $self = shift;
    $self->new_from_cols(@_);
}
# from Math::MatrixReal::Ext1 ([e-mail])
sub new_from_rows {
    my $self = shift;
    my $extra_args = ( @_ > 1 && ref($_[-1]) eq 'HASH' ) ? pop : {};
    $extra_args->{_type} = 'row';

    $self->_new_from_rows_or_cols(@_, $extra_args );
}

sub reshape {
    my ($self, $rows, $cols, $values) = @_;

    my @cols = ();
    my $p = 0;
    for my $c (1..$cols) {
        push @cols, [@{$values}[$p .. $p + $rows - 1]];
        $p += $rows;
    }

    return $self->new_from_cols( \@cols );
}

# from Math::MatrixReal::Ext1 ([e-mail])
sub _new_from_rows_or_cols {
    my $proto = shift;
    my $class = ref($proto) || $proto;
    my $ref_to_vectors = shift;

    # these additional args are internal at the moment,  but in the future the user could pass e.g. {pad=>1} to
    # request padding
    my $args = pop;
    my $vector_type = $args->{_type};
    die "Internal ".__PACKAGE__." error" unless $vector_type =~ /^(row|column)$/;

    # step back one frame because this private method is  not how the user called it
    my $caller_subname = (caller(1))[3];

    croak "$caller_subname: need a reference to an array of ${vector_type}s" unless reftype($ref_to_vectors) eq 'ARRAY';

    my @vectors = @{$ref_to_vectors};
    my $matrix;
    my $other_type = {row=>'column', column=>'row'}->{$vector_type};
    my %matrix_dim = (
        $vector_type => scalar( @vectors ),
        $other_type  => 0,  # we will correct this in a bit
    );

    # row and column indices are one based
    my $current_vector_count = 1;
    foreach my $current_vector (@vectors) {
        # dimension is one-based, so we're
        # starting with one here and incrementing
        # as we go.  The other dimension is fixed (for now, until
        # we add the 'pad' option), and gets set later
        my $ref = ref( $current_vector ) ;

        if ( $ref eq '' ) {
            # we hope this is a properly formatted Math::MatrixReal string,
            # but if not we just let the Math::MatrixReal die() do it's
            # thing
            $current_vector = $class->new_from_string( $current_vector );
        } elsif ( $ref eq 'ARRAY' ) {
            my @array = @$current_vector;
            croak "$caller_subname: one $vector_type you gave me was a ref to an array with no elements" unless @array;
            # we need to create the right kind of string based on whether
            # they said they were sending us rows or columns:
            if ($vector_type eq 'row') {
                $current_vector = $class->new_from_string( '[ '. join( " ", @array) ." ]\n" );
            } else {
                $current_vector = $class->new_from_string( '[ '. join( " ]\n[ ", @array) ." ]\n" );
            }
        } elsif ( $ref ne 'HASH' and
                ( $current_vector->isa('Math::MatrixReal') ||
                  $current_vector->isa('Math::MatrixComplex')
                ) ) {
            # it's already a Math::MatrixReal something.
            # we don't need to do anything, it will all
            # work out
        } else {
            # we have no idea, error time!
            croak "$caller_subname: I only know how to deal with array refs, strings, and things that inherit from Math::MatrixReal\n";
        }

        # starting now we know $current_vector isa Math::MatrixReal thingy
        my @vector_dims = $current_vector->dim;

        #die unless the appropriate dimension is 1
        croak "$caller_subname: I don't accept $other_type vectors" unless ($vector_dims[ $vector_type eq 'row' ? 0 : 1 ] == 1) ;

        # the other dimension is the length of our vector
        my $length =  $vector_dims[ $vector_type eq 'row' ? 1 : 0 ];

        # set the "other" dimension to the length of this
        # vector the first time through
        $matrix_dim{$other_type} ||= $length;

        # die unless length of this vector matches the first length
        croak "$caller_subname: one $vector_type has [$length] elements and another one had [$matrix_dim{$other_type}]--all of the ${vector_type}s passed in must have the same dimension"
              unless ($length == $matrix_dim{$other_type}) ;

        # create the matrix the first time through
        $matrix ||= $class->new($matrix_dim{row}, $matrix_dim{column});

        # step along the vector assigning the value of each element
        # to the correct place in the matrix we're building
        foreach my $element_index ( 1..$length ){
            # args for vector assignment:
            # initialize both to one and reset the correct
            # one below
            my ($v_r, $v_c) = (1,1);

            # args for matrix assignment
            my ($row_index, $column_index, $value);

            if ($vector_type eq 'row') {
                $row_index           = $current_vector_count;
                $v_c = $column_index = $element_index;
            } else {
                $v_r = $row_index    = $element_index;
                $column_index        = $current_vector_count;
            }
            $value = $current_vector->element($v_r, $v_c);
            $matrix->assign($row_index, $column_index, $value);
        }
        $current_vector_count ++ ;
    }
    return $matrix;
}

sub shadow
{
    croak "Usage: \$new_matrix = \$some_matrix->shadow();" if (@_ != 1);

    my ($matrix) = @_;

    return $matrix->new($matrix->[1],$matrix->[2]);
}

=over 4

=item * $matrix->display_precision($integer)

Sets the default precision when matrices are printed or stringified.
$matrix->display_precision(0) will only show the integer part of all the
entries of $matrix and $matrix->display_precision() will return to the default
scientific display notation. This method does not effect the precision of the
calculations.

=back

=cut

sub display_precision
{
    my ($self,$n) = @_;
    if (defined $n) {
        croak "Usage: \$matrix->display_precision(\$nonnegative_integer);" if ($n < 0);
        $self->[4] = int $n;
    } else {
        $self->[4] = undef;
    }
}

sub copy
{
    croak "Usage: \$matrix1->copy(\$matrix2);"
      if (@_ != 2);

    my ($matrix1,$matrix2) = @_;
    my ($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
    my ($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
    my ($i,$j);

    croak "Math::MatrixReal::copy(): matrix size mismatch" unless $rows1 == $rows2 && $cols1 == $cols2;

    for ( $i = 0; $i < $rows1; $i++ )
    {
        my $r1            = [];
        my $r2            = $matrix2->[0][$i];
        @$r1              = @$r2;              # Copy whole array directly
        $matrix1->[0][$i] = $r1;
    }
    if (defined $matrix2->[3]) # is an LR decomposition matrix!
    {
        $matrix1->[3] = $matrix2->[3]; # $sign
        $matrix1->[4] = $matrix2->[4]; # $perm_row
        $matrix1->[5] = $matrix2->[5]; # $perm_col
    }
}

sub clone
{
    croak "Usage: \$twin_matrix = \$some_matrix->clone();" if (@_ != 1);

    my($matrix) = @_;
    my($temp);

    $temp = $matrix->new($matrix->[1],$matrix->[2]);
    $temp->copy($matrix);
    return $temp;
}

## trace() : return the sum of the diagonal elements
sub trace {
    croak "Usage: \$trace = \$matrix->trace();" if (@_ != 1);

    my $matrix = shift;
    my($rows,$cols) = ($matrix->[1],$matrix->[2]);
    my $trace = 0;

    croak "Math::MatrixReal::trace(): matrix is not quadratic" unless ($rows == $cols);

    map { $trace += $matrix->[0][$_][$_] } (0 .. $cols-1);

    return $trace;
}
sub submatrix {
    my $self = shift;
    my ($x1, $y1, $x2, $y2) = @_;
    croak "Math::MatrixReal::submatrix(): indices must be positive integers"
        unless ($x1 >= 1 && $x2 >= 1 && $y1 >=1 && $y2 >=1 );
    my($rows,$cols) = ($self->[1],$self->[2]);
    my($sr,$sc)     = ( 1+abs($x1-$x2), 1+abs($y1-$y2) );
    my $submatrix = $self->new( $sr, $sc );

    for (my $i = $x1-1; $i < $x2; $i++ ) {
        for (my $j = $y1-1; $j < $y2; $j++ ) {
            $submatrix->[0][$i-($x1-1)][$j-($y1-1)] = $self->[0][$i][$j];
        }
    }
    return $submatrix;
}
## return the minor corresponding to $r and $c
## eliminate row $r and col $c , and return the $r-1 by $c-1 matrix
sub minor {
    croak "Usage: \$minor = \$matrix->minor(\$r,\$c);" unless (@_ == 3);
    my ($matrix,$r,$c) = @_;
    my ($rows,$cols) = $matrix->dim();

    croak "Math::MatrixReal::minor(): \$matrix must be at least 2x2"
        unless ($rows > 1 and $cols > 1);
    croak "Math::MatrixReal::minor(): $r and $c must be positive"
        unless ($r > 0 and $c > 0 );
    croak "Math::MatrixReal::minor(): matrix has no $r,$c element"
        unless ($r <= $rows and $c <= $cols );

    my $minor = new Math::MatrixReal($rows-1,$cols-1);
    my ($i,$j) = (0,0);

    ## assign() might have been easier, but this should be faster
    ## the element can be in any of 4 regions compared to the eliminated
    ## row and col:
    ## above and to the left, above and to the right
    ## below and to the left, below and to the right

    for(; $i < $rows; $i++){
        for(;$j < $rows; $j++ ){
            if( $i >= $r && $j >= $c ){
                $minor->[0][$i-1][$j-1] = $matrix->[0][$i][$j];
            } elsif ( $i >= $r && $j < $c ){
                $minor->[0][$i-1][$j] = $matrix->[0][$i][$j];
            } elsif ( $i < $r && $j < $c ){
                $minor->[0][$i][$j] = $matrix->[0][$i][$j];
            } elsif ( $i < $r && $j >= $c ){
                $minor->[0][$i][$j-1] = $matrix->[0][$i][$j];
            } else {
                croak "Very bad things";
            }
        }
        $j = 0;
    }
    return ($minor);
}
sub swap_col {
    croak "Usage: \$matrix->swap_col(\$col1,\$col2); " unless (@_ == 3);
    my ($matrix,$col1,$col2) = @_;
    my ($rows,$cols) = $matrix->dim();
    my (@temp);

    croak "Math::MatrixReal::swap_col(): col index is not valid"
        unless ( $col1 <= $cols && $col2 <= $cols &&
             $col1 == int($col1) &&
             $col2 == int($col2) );
    $col1--;$col2--;
    for(my $i=0;$i < $rows;$i++){
        $temp[$i] = $matrix->[0][$i][$col1];
        $matrix->[0][$i][$col1] = $matrix->[0][$i][$col2];
        $matrix->[0][$i][$col2] =  $temp[$i];
    }
}
sub swap_row {
    croak "Usage: \$matrix->swap_row(\$row1,\$row2); " unless (@_ == 3);
    my ($matrix,$row1,$row2) = @_;
    my ($rows,$cols) = $matrix->dim();
    my (@temp);

    croak "Math::MatrixReal::swap_row(): row index is not valid"
        unless ( $row1 <= $rows && $row2 <= $rows &&
             $row1 == int($row1) &&
             $row2 == int($row2) );
    $row1--;$row2--;
    for(my $j=0;$j < $cols;$j++){
        $temp[$j] = $matrix->[0][$row1][$j];
        $matrix->[0][$row1][$j] = $matrix->[0][$row2][$j];
        $matrix->[0][$row2][$j] =  $temp[$j];
    }
}

sub assign_row {
    croak "Usage: \$matrix->assign_row(\$row,\$row_vec);"  unless (@_ == 3);
    my ($matrix,$row,$row_vec) = @_;
    my ($rows1,$cols1) = $matrix->dim();
    my ($rows2,$cols2) = $row_vec->dim();

    croak "Math::MatrixReal::assign_row(): number of columns mismatch" if ($cols1 != $cols2);
    croak "Math::MatrixReal::assign_row(): not a row vector" unless( $rows2 == 1);

    @{$matrix->[0][--$row]} = @{$row_vec->[0][0]};
    return $matrix;
}
# returns the number of zeroes in a row
sub _count_zeroes_row {
    my ($matrix) = @_;
    my ($rows,$cols) = $matrix->dim();
    my $count = 0;
    croak "_count_zeroes_row(): only 1 row, buddy" unless ($rows == 1);

    map { $count++ unless $matrix->[0][0][$_] } (0 .. $cols-1);
    return $count;
}
## divide a row by it's largest abs() element
sub _normalize_row {
    my ($matrix) = @_;
    my ($rows,$cols) = $matrix->dim();
    my $new_row = Math::MatrixReal->new(1,$cols);

    my $big = abs($matrix->[0][0][0]);
    for(my $j=0;$j < $cols; $j++ ){
        $big = $big < abs($matrix->[0][0][$j])
            ? abs($matrix->[0][0][$j]) : $big;
    }
    next unless $big;
    # now $big is biggest element in row
    for(my $j = 0;$j < $cols; $j++ ){
        $new_row->[0][0][$j]  = $matrix->[0][0][$j] / $big;
    }
    return $new_row;

}

sub cofactor {
    my ($matrix) = @_;
    my ($rows,$cols) = $matrix->dim();

    croak "Math::MatrixReal::cofactor(): Matrix is not quadratic"
        unless ($rows == $cols);

    # black magic ahead
    my $cofactor = $matrix->each(
        sub {
            my($v,$i,$j) = @_;
            ($i+$j) % 2 == 0 ? $matrix->minor($i,$j)->det() : -1*$matrix->minor($i,$j)->det();
        });
    return ($cofactor);
}

sub adjoint {
    my ($matrix) = @_;
    return ~($matrix->cofactor);
}

sub row
{
    croak "Usage: \$row_vector = \$matrix->row(\$row);"
      if (@_ != 2);

    my($matrix,$row) = @_;
    my($rows,$cols) = ($matrix->[1],$matrix->[2]);
    my($temp);

    croak "Math::MatrixReal::row(): row index out of range" if ($row < 1 || $row > $rows);

    $row--;
    $temp = $matrix->new(1,$cols);
    for ( my $j = 0; $j < $cols; $j++ )
    {
        $temp->[0][0][$j] = $matrix->[0][$row][$j];
    }
    return($temp);
}
sub col{ return (shift)->column(shift) }
sub column
{
    croak "Usage: \$column_vector = \$matrix->column(\$column);" if (@_ != 2);

    my($matrix,$col) = @_;
    my($rows,$cols) = ($matrix->[1],$matrix->[2]);
    #my($temp);
    #my($i);
    my $col_vector;

    croak "Math::MatrixReal::column(): column index out of range" if ($col < 1 || $col > $cols);

    $col--;
    $col_vector = $matrix->new($rows,1);

    map { $col_vector->[0][$_][0] = $matrix->[0][$_][$col] } (0 .. $rows-1);

    return $col_vector;
}

sub as_list
{
    croak "Usage: \$matrix->as_list();" if (@_ != 1);

    my($self) = @_;
    my($rows,$cols) = ($self->[1], $self->[2]);
    my @list;
    for(my $i = 0; $i < $rows; $i++ ){
        for(my $j = 0; $j < $rows; $j++){
            push @list, $self->[0][$i][$j];
        }
    }
    return @list;
}

sub _undo_LR
{
    croak "Usage: \$matrix->_undo_LR();"
      if (@_ != 1);

    my($self) = @_;

    undef $self->[3];
    undef $self->[4];
    undef $self->[5];
}
# brrr
sub zero
{
    croak "Usage: \$matrix->zero();" if (@_ != 1);

    my ($self) = @_;
    my ($rows,$cols) = ($self->[1],$self->[2]);

    $self->_undo_LR();

    # zero out first row
    map {  $self->[0][0][$_] = 0.0              } (0 .. $cols-1);

    # copy that to the other rows
    map {  @{$self->[0][$_]} = @{$self->[0][0]} } (0 .. $rows-1);

    return $self;
}

sub one
{
    croak "Usage: \$matrix->one();" if (@_ != 1);

    my ($self) = @_;
    my ($rows,$cols) = ($self->[1],$self->[2]);

    $self->zero(); # We rely on zero() efficiency

    map { $self->[0][$_][$_] = 1.0 } (0 .. $rows-1);

    return $self;
}

sub assign
{
    croak "Usage: \$matrix->assign(\$row,\$column,\$value);" if (@_ != 4);

    my($self,$row,$col,$value) = @_;
    my($rows,$cols) = ($self->[1],$self->[2]);

    croak "Math::MatrixReal::assign(): row index out of range" if (($row < 1) || ($row > $rows));
    croak "Math::MatrixReal::assign(): column index out of range" if (($col < 1) || ($col > $cols));

    $self->_undo_LR();
    $self->[0][--$row][--$col] = $value;
}

sub element
{
    croak "Usage: \$value = \$matrix->element(\$row,\$column);" if (@_ != 3);

    my($self,$row,$col) = @_;
    my($rows,$cols) = ($self->[1],$self->[2]);

    croak "Math::MatrixReal::element(): row index out of range" if (($row < 1) || ($row > $rows));
    croak "Math::MatrixReal::element(): column index out of range" if (($col < 1) || ($col > $cols));

    return( $self->[0][--$row][--$col] );
}

sub dim  #  returns dimensions of a matrix
{
    croak "Usage: (\$rows,\$columns) = \$matrix->dim();" if (@_ != 1);

    my($matrix) = @_;

    return( $matrix->[1], $matrix->[2] );
}

sub norm_one  #  maximum of sums of each column
{
    croak "Usage: \$norm_one = \$matrix->norm_one();" if (@_ != 1);

    my($self) = @_;
    my($rows,$cols) = ($self->[1],$self->[2]);

    my $max = 0.0;
    for (my $j = 0; $j < $cols; $j++)
    {
        my $sum = 0.0;
        for (my $i = 0; $i < $rows; $i++)
        {
            $sum += abs( $self->[0][$i][$j] );
        }
        $max = $sum if ($sum > $max);
    }
    return($max);
}
## sum of absolute value of every element
sub norm_sum {
    croak "Usage: \$norm_sum = \$matrix->norm_sum();" unless (@_ == 1);
    my ($matrix) = @_;
    my $norm = 0;
    $matrix->each( sub { $norm+=abs(shift); } );
    return $norm;
}
sub norm_frobenius {
    my ($m) = @_;
    my ($r,$c) = $m->dim;
    my $s=0;

    $m->each( sub { $s+=abs(shift)**2 } );
    return sqrt($s);
}

# Vector Norm
sub norm_p {
    my ($v,$p) = @_;
    # sanity check on $p
    croak "Math::MatrixReal:norm_p: argument must be a row or column vector"
        unless ( $v->is_row_vector || $v->is_col_vector );
    croak "Math::MatrixReal::norm_p: $p must be >= 1"
        unless ($p =~ m/Inf(inity)?/i || $p >= 1);

    if( $p =~ m/^(Inf|Infinity)$/i ){
        my $max = $v->element(1,1);
        $v->each ( sub { my $x=abs(shift); $max = $x if( $x > $max ); } );
        return $max;
    }

    my $s=0;
    $v->each( sub { $s+= (abs(shift))**$p; } );
    return $s ** (1/$p);

}
sub norm_max  #  maximum of sums of each row
{
    croak "Usage: \$norm_max = \$matrix->norm_max();" if (@_ != 1);

    my($self) = @_;
    my($rows,$cols) = ($self->[1],$self->[2]);

    my $max = 0.0;
    for (my $i = 0; $i < $rows; $i++)
    {
        my $sum = 0.0;
        for (my $j = 0; $j < $cols; $j++)
        {
            $sum += abs( $self->[0][$i][$j] );
        }
        $max = $sum if ($sum > $max);
    }
    return($max);
}

sub negate
{
    croak "Usage: \$matrix1->negate(\$matrix2);"
      if (@_ != 2);

    my($matrix1,$matrix2) = @_;
    my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
    my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);

    croak "Math::MatrixReal::negate(): matrix size mismatch"
      unless (($rows1 == $rows2) && ($cols1 == $cols2));

    $matrix1->_undo_LR();

    for (my $i = 0; $i < $rows1; $i++ )
    {
        for (my $j = 0; $j < $cols1; $j++ )
        {
            $matrix1->[0][$i][$j] = -($matrix2->[0][$i][$j]);
        }
    }
}
## each(): evaluate a coderef on each element and return a new matrix
## of said
sub each {
    croak "Usage: \$new_matrix = \$matrix->each( \&sub );" unless (@_ == 2 );
    my($matrix,$function) = @_;
    my($rows,$cols) = ($matrix->[1],$matrix->[2]);
    my($new_matrix) = $matrix->clone();

    croak "Math::MatrixReal::each(): argument is not a sub reference" unless ref($function);
    $new_matrix->_undo_LR();

    for (my $i = 0; $i < $rows; $i++ ) {
        for (my $j = 0; $j < $cols; $j++ ) {
            no strict 'refs';
            # $i,$j are 1-based as of 1.7
            $new_matrix->[0][$i][$j] = &{ $function }($matrix->[0][$i][$j],$i+1,$j+1) ;
        }
    }
    return ($new_matrix);
}
## each_diag(): same as each() but only diag elements
sub each_diag {
    croak "Usage: \$new_matrix = \$matrix->each_diag( \&sub );" unless (@_ == 2 );
    my($matrix,$function) = @_;
    my($rows,$cols) = ($matrix->[1],$matrix->[2]);
    my($new_matrix) = $matrix->clone();

    croak "Math::MatrixReal::each(): argument is not a sub reference" unless ref($function);
    croak "Matrix is not quadratic" unless ($rows == $cols);

    $new_matrix->_undo_LR();

    for (my $i = 0; $i < $rows; $i++ ) {
        for (my $j = 0; $j < $cols; $j++ ) {
            next unless ($i == $j);
            no strict 'refs';
            # $i,$j are 1-based as of 1.7
            $new_matrix->[0][$i][$j] = &{ $function }($matrix->[0][$i][$j],$i+1,$j+1) ;
        }
    }
    return ($new_matrix);
}

## Make computing the inverse more user friendly
sub inverse {
    croak "Usage: \$inverse = \$matrix->inverse();" unless (@_ == 1);
    my ($matrix) = @_;
    return $matrix->decompose_LR->invert_LR;
}

sub transpose {

    croak "Usage: \$matrix1->transpose(\$matrix2);" if (@_ != 2);

    my($matrix1,$matrix2) = @_;
    my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
    my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);

    croak "Math::MatrixReal::transpose(): matrix size mismatch"
      unless (($rows1 == $cols2) && ($cols1 == $rows2));

    $matrix1->_undo_LR();

    if ($rows1 == $cols1)
    {
        # more complicated to make in-place possible!

        for (my $i = 0; $i < $rows1; $i++)
        {
            for (my $j = ($i + 1); $j < $cols1; $j++)
            {
                my $swap              = $matrix2->[0][$i][$j];
                $matrix1->[0][$i][$j] = $matrix2->[0][$j][$i];
                $matrix1->[0][$j][$i] = $swap;
            }
            $matrix1->[0][$i][$i] = $matrix2->[0][$i][$i];
        }
    } else {                                # ($rows1 != $cols1)
        for (my $i = 0; $i < $rows1; $i++)
        {
            for (my $j = 0; $j < $cols1; $j++)
            {
                $matrix1->[0][$i][$j] = $matrix2->[0][$j][$i];
            }
        }
    }
}

sub add
{
    croak "Usage: \$matrix1->add(\$matrix2,\$matrix3);" if (@_ != 3);

    my($matrix1,$matrix2,$matrix3) = @_;
    my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
    my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
    my($rows3,$cols3) = ($matrix3->[1],$matrix3->[2]);

    croak "Math::MatrixReal::add(): matrix size mismatch"
      unless (($rows1 == $rows2) && ($rows1 == $rows3) &&
              ($cols1 == $cols2) && ($cols1 == $cols3));

    $matrix1->_undo_LR();

    for ( my $i = 0; $i < $rows1; $i++ )
    {
        for ( my $j = 0; $j < $cols1; $j++ )
        {
            $matrix1->[0][$i][$j] = $matrix2->[0][$i][$j] + $matrix3->[0][$i][$j];
        }
    }
}

sub subtract
{
    croak "Usage: \$matrix1->subtract(\$matrix2,\$matrix3);" if (@_ != 3);

    my($matrix1,$matrix2,$matrix3) = @_;
    my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
    my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
    my($rows3,$cols3) = ($matrix3->[1],$matrix3->[2]);

    croak "Math::MatrixReal::subtract(): matrix size mismatch"
      unless (($rows1 == $rows2) && ($rows1 == $rows3) &&
              ($cols1 == $cols2) && ($cols1 == $cols3));

    $matrix1->_undo_LR();

    for ( my $i = 0; $i < $rows1; $i++ )
    {
        for ( my $j = 0; $j < $cols1; $j++ )
        {
            $matrix1->[0][$i][$j] = $matrix2->[0][$i][$j] - $matrix3->[0][$i][$j];
        }
    }
}

sub multiply_scalar
{
    croak "Usage: \$matrix1->multiply_scalar(\$matrix2,\$scalar);"
      if (@_ != 3);

    my($matrix1,$matrix2,$scalar) = @_;
    my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
    my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);

    croak "Math::MatrixReal::multiply_scalar(): matrix size mismatch"
      unless (($rows1 == $rows2) && ($cols1 == $cols2));

    $matrix1->_undo_LR();

    for ( my $i = 0; $i < $rows1; $i++ )
    {
        map { $matrix1->[0][$i][$_] = $matrix2->[0][$i][$_] * $scalar } (0 .. $cols1-1);
    }
}

sub multiply
{
    croak "Usage: \$product_matrix = \$matrix1->multiply(\$matrix2);"
      if (@_ != 2);

    my($matrix1,$matrix2) = @_;
    my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
    my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);

    croak "Math::MatrixReal::multiply(): matrix size mismatch" unless ($cols1 == $rows2);

    my $temp = $matrix1->new($rows1,$cols2);
    for (my $i = 0; $i < $rows1; $i++ )
    {
        for (my $j = 0; $j < $cols2; $j++ )
        {
            my $sum = 0.0;
            for (my $k = 0; $k < $cols1; $k++ )
            {
                $sum += ( $matrix1->[0][$i][$k] * $matrix2->[0][$k][$j] );
            }
            $temp->[0][$i][$j] = $sum;
        }
    }
    return($temp);
}

sub exponent {
    croak "Usage: \$matrix_exp = \$matrix1->exponent(\$integer);" if(@_ != 2 );
    my($matrix,$argument) = @_;
    my($rows,$cols) = ($matrix->[1],$matrix->[2]);
    my($name) = "'**'";
    my($temp) = $matrix->clone();

    croak "Matrix is not quadratic" unless ($rows == $cols);
    croak "Exponent must be integer" unless ($argument =~ m/^[+-]?\d+$/ );

    return($matrix) if ($argument == 1);

    $temp->_undo_LR();

    # negative exponent is (A^-1)^n
    if( $argument < 0 ){
        my $LR = $matrix->decompose_LR();
        my $inverse = $LR->invert_LR();
        unless (defined $inverse){
            carp "Matrix has no inverse";
            return undef;
        }
        $temp = $inverse->clone();
        if( $inverse ){
            return($inverse) if ($argument == -1);
                for( 2 .. abs($argument) ){
                    $temp = multiply($inverse,$temp);
                    }
                return($temp);
            } else {
               # TODO: is this the right behaviour?
               carp "Cannot compute negative exponent, inverse does not exist";
               return undef;
        }
    # matrix to zero power is identity matrix
    } elsif( $argument == 0 ){
        $temp->one();
        return ($temp);
    }

    # if it is diagonal, just raise diagonal entries to power
    if( $matrix->is_diagonal() ){
        $temp = $temp->each_diag( sub { (shift)**$argument } );
        return ($temp);

    } else {
        for( 2 .. $argument ){
            $temp = multiply($matrix,$temp);
        }
        return ($temp);
    }
}

sub min
{

    if ( @_ == 1 ) {
        my $matrix = shift;

        croak "Usage: \$minimum = Math::MatrixReal::min(\$number1,\$number2) or $matrix->min" if (@_ > 2);
        croak "invalid" unless ref $matrix eq 'Math::MatrixReal';

        my $min = $matrix->element(1,1);
        $matrix->each( sub { my ($e,$i,$j) = @_; $min = $e if $e < $min; } );
        return $min;
    }
    $_[0] < $_[1] ? $_[0] : $_[1];
}

sub max
{
    if ( @_ == 1 ) {
        my $matrix = shift;

        croak "Usage: \$maximum = Math::MatrixReal::max(\$number1,\$number2) or $matrix->max" if (@_ > 2);
        croak "Math::MatrixReal::max(\$matrix) \$matrix is not a Math::MatrixReal matrix" unless ref $matrix eq 'Math::MatrixReal';

        my $max = $matrix->element(1,1);
        $matrix->each( sub { my ($e,$i,$j) = @_; $max = $e if $e > $max; } );
        return $max;
    }

    $_[0] > $_[1] ? $_[0] : $_[1];
}

sub kleene
{
    croak "Usage: \$minimal_cost_matrix = \$cost_matrix->kleene();" if (@_ != 1);

    my($matrix) = @_;
    my($rows,$cols) = ($matrix->[1],$matrix->[2]);

    croak "Math::MatrixReal::kleene(): matrix is not quadratic" unless ($rows == $cols);

    my $temp = $matrix->new($rows,$cols);
    $temp->copy($matrix);
    $temp->_undo_LR();
    my $n = $rows;
    for ( my $i = 0; $i < $n; $i++ )
    {
        $temp->[0][$i][$i] = min( $temp->[0][$i][$i] , 0 );
    }
    for ( my $k = 0; $k < $n; $k++ )
    {
        for ( my $i = 0; $i < $n; $i++ )
        {
            for ( my $j = 0; $j < $n; $j++ )
            {
                $temp->[0][$i][$j] = min( $temp->[0][$i][$j] ,
                                        ( $temp->[0][$i][$k] +
                                          $temp->[0][$k][$j] ) );
            }
        }
    }
    return($temp);
}

sub normalize
{
    croak "Usage: (\$norm_matrix,\$norm_vector) = \$matrix->normalize(\$vector);"
      if (@_ != 2);

    my($matrix,$vector) = @_;
    my($rows,$cols) = ($matrix->[1],$matrix->[2]);
    my($norm_matrix,$norm_vector);
    my($max,$val);
    my($i,$j,$n);

    croak "Math::MatrixReal::normalize(): matrix is not quadratic"
      unless ($rows == $cols);

    $n = $rows;

    croak "Math::MatrixReal::normalize(): vector is not a column vector"
      unless ($vector->[2] == 1);

    croak "Math::MatrixReal::normalize(): matrix and vector size mismatch"
      unless ($vector->[1] == $n);

    $norm_matrix = $matrix->new($n,$n);
    $norm_vector = $vector->new($n,1);

    $norm_matrix->copy($matrix);
    $norm_vector->copy($vector);

    $norm_matrix->_undo_LR();

    for ( $i = 0; $i < $n; $i++ )
    {
        $max = abs($norm_vector->[0][$i][0]);
        for ( $j = 0; $j < $n; $j++ )
        {
            $val = abs($norm_matrix->[0][$i][$j]);
            if ($val > $max) { $max = $val; }
        }
        if ($max != 0)
        {
            $norm_vector->[0][$i][0] /= $max;
            for ( $j = 0; $j < $n; $j++ )
            {
                $norm_matrix->[0][$i][$j] /= $max;
            }
        }
    }
    return($norm_matrix,$norm_vector);
}

sub decompose_LR
{
    croak "Usage: \$LR_matrix = \$matrix->decompose_LR();"
      if (@_ != 1);

    my($matrix) = @_;
    my($rows,$cols) = ($matrix->[1],$matrix->[2]);
    my($perm_row,$perm_col);
    my($row,$col,$max);
    my($i,$j,$k,$n);
    my($sign) = 1;
    my($swap);
    my($temp);

    croak "Math::MatrixReal::decompose_LR(): matrix is not quadratic"
      unless ($rows == $cols);

    $temp = $matrix->new($rows,$cols);
    $temp->copy($matrix);
    $n = $rows;
    $perm_row = [ ];
    $perm_col = [ ];
    for ( $i = 0; $i < $n; $i++ )
    {
        $perm_row->[$i] = $i;
        $perm_col->[$i] = $i;
    }
    NONZERO:
    for ( $k = 0; $k < $n; $k++ ) # use Gauss's algorithm:
    {
        # complete pivot-search:

        $max = 0;
        for ( $i = $k; $i < $n; $i++ )
        {
            for ( $j = $k; $j < $n; $j++ )
            {
                if (($swap = abs($temp->[0][$i][$j])) > $max)
                {
                    $max = $swap;
                    $row = $i;
                    $col = $j;
                }
            }
        }
        last NONZERO if ($max == 0); # (all remaining elements are zero)
        if ($k != $row) # swap row $k and row $row:
        {
            $sign = -$sign;
            $swap             = $perm_row->[$k];
            $perm_row->[$k]   = $perm_row->[$row];
            $perm_row->[$row] = $swap;
            for ( $j = 0; $j < $n; $j++ )
            {
                # (must run from 0 since L has to be swapped too!)

                $swap                = $temp->[0][$k][$j];
                $temp->[0][$k][$j]   = $temp->[0][$row][$j];
                $temp->[0][$row][$j] = $swap;
            }
        }
        if ($k != $col) # swap column $k and column $col:
        {
            $sign = -$sign;
            $swap             = $perm_col->[$k];
            $perm_col->[$k]   = $perm_col->[$col];
            $perm_col->[$col] = $swap;
            for ( $i = 0; $i < $n; $i++ )
            {
                $swap                = $temp->[0][$i][$k];
                $temp->[0][$i][$k]   = $temp->[0][$i][$col];
                $temp->[0][$i][$col] = $swap;
            }
        }
        for ( $i = ($k + 1); $i < $n; $i++ )
        {
            # scan the remaining rows, add multiples of row $k to row $i:

            $swap = $temp->[0][$i][$k] / $temp->[0][$k][$k];
            if ($swap != 0)
            {
                # calculate a row of matrix R:

                for ( $j = ($k + 1); $j < $n; $j++ )
                {
                    $temp->[0][$i][$j] -= $temp->[0][$k][$j] * $swap;
                }

                # store matrix L in same matrix as R:

                $temp->[0][$i][$k] = $swap;
            }
        }
    }
    $temp->[3] = $sign;
    $temp->[4] = $perm_row;
    $temp->[5] = $perm_col;
    return($temp);
}

sub solve_LR
{
    croak "Usage: (\$dimension,\$x_vector,\$base_matrix) = \$LR_matrix->solve_LR(\$b_vector);"
      if (@_ != 2);

    my($LR_matrix,$b_vector) = @_;
    my($rows,$cols) = ($LR_matrix->[1],$LR_matrix->[2]);
    my($dimension,$x_vector,$base_matrix);
    my($perm_row,$perm_col);
    my($y_vector,$sum);
    my($i,$j,$k,$n);

    croak "Math::MatrixReal::solve_LR(): not an LR decomposition matrix"
      unless ((defined $LR_matrix->[3]) && ($rows == $cols));

    $n = $rows;

    croak "Math::MatrixReal::solve_LR(): vector is not a column vector"
      unless ($b_vector->[2] == 1);

    croak "Math::MatrixReal::solve_LR(): matrix and vector size mismatch"
      unless ($b_vector->[1] == $n);

    $perm_row = $LR_matrix->[4];
    $perm_col = $LR_matrix->[5];

    $x_vector    =   $b_vector->new($n,1);
    $y_vector    =   $b_vector->new($n,1);
    $base_matrix = $LR_matrix->new($n,$n);

    # calculate "x" so that LRx = b  ==>  calculate Ly = b, Rx = y:

    for ( $i = 0; $i < $n; $i++ ) # calculate $y_vector:
    {
        $sum = $b_vector->[0][($perm_row->[$i])][0];
        for ( $j = 0; $j < $i; $j++ )
        {
            $sum -= $LR_matrix->[0][$i][$j] * $y_vector->[0][$j][0];
        }
        $y_vector->[0][$i][0] = $sum;
    }

    $dimension = 0;
    for ( $i = ($n - 1); $i >= 0; $i-- ) # calculate $x_vector:
    {
        if ($LR_matrix->[0][$i][$i] == 0)
        {
            if ($y_vector->[0][$i][0] != 0)
            {
                return(); # a solution does not exist!
            }
            else
            {
                $dimension++;
                $x_vector->[0][($perm_col->[$i])][0] = 0;
            }
        } else {
            $sum = $y_vector->[0][$i][0];
            for ( $j = ($i + 1); $j < $n; $j++ )
            {
                $sum -= $LR_matrix->[0][$i][$j] *
                    $x_vector->[0][($perm_col->[$j])][0];
            }
            $x_vector->[0][($perm_col->[$i])][0] =
                $sum / $LR_matrix->[0][$i][$i];
        }
    }
    if ($dimension)
    {
        if ($dimension == $n)
        {
            $base_matrix->one();
        } else {
            for ( $k = 0; $k < $dimension; $k++ )
            {
                $base_matrix->[0][($perm_col->[($n-$k-1)])][$k] = 1;
                for ( $i = ($n-$dimension-1); $i >= 0; $i-- )
                {
                    $sum = 0;
                    for ( $j = ($i + 1); $j < $n; $j++ )
                    {
                        $sum -= $LR_matrix->[0][$i][$j] *
                            $base_matrix->[0][($perm_col->[$j])][$k];
                    }
                    $base_matrix->[0][($perm_col->[$i])][$k] =
                        $sum / $LR_matrix->[0][$i][$i];
                }
            }
        }
    }
    return( $dimension, $x_vector, $base_matrix );
}

sub invert_LR
{
    croak "Usage: \$inverse_matrix = \$LR_matrix->invert_LR();"
      if (@_ != 1);

    my($matrix) = @_;
    my($rows,$cols) = ($matrix->[1],$matrix->[2]);
    my($inv_matrix,$x_vector,$y_vector);
    my($i,$j,$n);

    croak "Math::MatrixReal::invert_LR(): not an LR decomposition matrix"
      unless ((defined $matrix->[3]) && ($rows == $cols));

    $n = $rows;
    #print Dumper [ $matrix ];
    if ($matrix->[0][$n-1][$n-1] != 0)
    {
        $inv_matrix = $matrix->new($n,$n);
        $y_vector   = $matrix->new($n,1);
        for ( $j = 0; $j < $n; $j++ )
        {
            if ($j > 0)
            {
                $y_vector->[0][$j-1][0] = 0;
            }
            $y_vector->[0][$j][0] = 1;
            if (($rows,$x_vector,$cols) = $matrix->solve_LR($y_vector))
            {
                for ( $i = 0; $i < $n; $i++ )
                {
                    $inv_matrix->[0][$i][$j] = $x_vector->[0][$i][0];
                }
            } else {
                die "Math::MatrixReal::invert_LR(): unexpected error - please inform author!\n";
            }
        }
        return($inv_matrix);
    } else {
        warn __PACKAGE__ . qq{: matrix not invertible\n};
        return;
    }
}

sub condition
{
    # 1st matrix MUST be the inverse of 2nd matrix (or vice-versa)
    # for a meaningful result!

    # make this work when given no args

    croak "Usage: \$condition = \$matrix->condition(\$inverse_matrix);" if (@_ != 2);

    my($matrix1,$matrix2) = @_;
    my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
    my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);

    croak "Math::MatrixReal::condition(): 1st matrix is not quadratic"
      unless ($rows1 == $cols1);

    croak "Math::MatrixReal::condition(): 2nd matrix is not quadratic"
      unless ($rows2 == $cols2);

    croak "Math::MatrixReal::condition(): matrix size mismatch"
      unless (($rows1 == $rows2) && ($cols1 == $cols2));

    return( $matrix1->norm_one() * $matrix2->norm_one() );
}

## easy to use determinant
## very fast if matrix is diagonal or triangular

sub det {
    croak "Usage: \$determinant = \$matrix->det_LR();" unless (@_ == 1);
    my ($matrix) = @_;
    my ($rows,$cols) = $matrix->dim();
    my $det = 1;

    croak "Math::MatrixReal::det(): Matrix is not quadratic"
        unless ($rows == $cols);

    # diagonal will match too
    if( $matrix->is_upper_triangular() ){
        $matrix->each_diag( sub { $det*=shift; } );
    } elsif ( $matrix->is_lower_triangular() ){
        $matrix->each_diag( sub { $det*=shift; } );
    } else {
        return $matrix->decompose_LR->det_LR();
    }
    return $det;
}

sub det_LR  #  determinant of LR decomposition matrix
{
    croak "Usage: \$determinant = \$LR_matrix->det_LR();"
      if (@_ != 1);

    my($matrix) = @_;
    my($rows,$cols) = ($matrix->[1],$matrix->[2]);
    my($k,$det);

    croak "Math::MatrixReal::det_LR(): not an LR decomposition matrix"
      unless ((defined $matrix->[3]) && ($rows == $cols));

    $det = $matrix->[3];
    for ( $k = 0; $k < $rows; $k++ )
    {
        $det *= $matrix->[0][$k][$k];
    }
    return($det);
}

sub rank_LR {
    return (shift)->order_LR;
}

sub order_LR  #  order of LR decomposition matrix (number of non-zero equations)
{
    croak "Usage: \$order = \$LR_matrix->order_LR();"
      if (@_ != 1);

    my($matrix) = @_;
    my($rows,$cols) = ($matrix->[1],$matrix->[2]);
    my($order);

    croak "Math::MatrixReal::order_LR(): not an LR decomposition matrix"
      unless ((defined $matrix->[3]) && ($rows == $cols));

    ZERO:
    for ( $order = ($rows - 1); $order >= 0; $order-- )
    {
        last ZERO if ($matrix->[0][$order][$order] != 0);
    }
    return(++$order);
}

sub scalar_product
{
    croak "Usage: \$scalar_product = \$vector1->scalar_product(\$vector2);"
      if (@_ != 2);

    my($vector1,$vector2) = @_;
    my($rows1,$cols1) = ($vector1->[1],$vector1->[2]);
    my($rows2,$cols2) = ($vector2->[1],$vector2->[2]);

    croak "Math::MatrixReal::scalar_product(): 1st vector is not a column vector"
      unless ($cols1 == 1);

    croak "Math::MatrixReal::scalar_product(): 2nd vector is not a column vector"
      unless ($cols2 == 1);

    croak "Math::MatrixReal::scalar_product(): vector size mismatch"
      unless ($rows1 == $rows2);

    my $sum = 0;
    map { $sum +=  $vector1->[0][$_][0] * $vector2->[0][$_][0] } ( 0 .. $rows1-1);
    return $sum;
}

sub vector_product
{
    croak "Usage: \$vector_product = \$vector1->vector_product(\$vector2);" if (@_ != 2);

    my($vector1,$vector2) = @_;
    my($rows1,$cols1) = ($vector1->[1],$vector1->[2]);
    my($rows2,$cols2) = ($vector2->[1],$vector2->[2]);
    my($temp);
    my($n);

    croak "Math::MatrixReal::vector_product(): 1st vector is not a column vector"
      unless ($cols1 == 1);

    croak "Math::MatrixReal::vector_product(): 2nd vector is not a column vector"
      unless ($cols2 == 1);

    croak "Math::MatrixReal::vector_product(): vector size mismatch"
      unless ($rows1 == $rows2);

    $n = $rows1;

    croak "Math::MatrixReal::vector_product(): only defined for 3 dimensions"
      unless ($n == 3);

    $temp = $vector1->new($n,1);
    $temp->[0][0][0] = $vector1->[0][1][0] * $vector2->[0][2][0] -
                       $vector1->[0][2][0] * $vector2->[0][1][0];
    $temp->[0][1][0] = $vector1->[0][2][0] * $vector2->[0][0][0] -
                       $vector1->[0][0][0] * $vector2->[0][2][0];
    $temp->[0][2][0] = $vector1->[0][0][0] * $vector2->[0][1][0] -
                       $vector1->[0][1][0] * $vector2->[0][0][0];
    return($temp);
}

sub length
{
    croak "Usage: \$length = \$vector->length();" if (@_ != 1);

    my($vector) = @_;
    my($rows,$cols) = ($vector->[1],$vector->[2]);
    my($k,$comp,$sum);

    croak "Math::MatrixReal::length(): vector is not a row or column vector"
      unless ($cols == 1 || $rows ==1 );

    $vector = ~$vector if ($rows == 1 );

    $sum = 0;
    for ( $k = 0; $k < $rows; $k++ )
    {
        $comp = $vector->[0][$k][0];
        $sum += $comp * $comp;
    }
    return sqrt $sum;
}

sub _init_iteration
{
    croak "Usage: \$which_norm = \$matrix->_init_iteration();"
      if (@_ != 1);

    my($matrix) = @_;
    my($rows,$cols) = ($matrix->[1],$matrix->[2]);
    my($ok,$max,$sum,$norm);
    my($i,$j,$n);

    croak "Math::MatrixReal::_init_iteration(): matrix is not quadratic"
      unless ($rows == $cols);

    $ok = 1;
    $n = $rows;
    for ( $i = 0; $i < $n; $i++ )
    {
        if ($matrix->[0][$i][$i] == 0) { $ok = 0; }
    }
    if ($ok)
    {
        $norm = 1; # norm_one
        $max = 0;
        for ( $j = 0; $j < $n; $j++ )
        {
            $sum = 0;
            for ( $i = 0; $i < $j; $i++ )
            {
                $sum += abs($matrix->[0][$i][$j]);
            }
            for ( $i = ($j + 1); $i < $n; $i++ )
            {
                $sum += abs($matrix->[0][$i][$j]);
            }
            $sum /= abs($matrix->[0][$j][$j]);
            if ($sum > $max) { $max = $sum; }
        }
        $ok = ($max < 1);
        unless ($ok)
        {
            $norm = -1; # norm_max
            $max = 0;
            for ( $i = 0; $i < $n; $i++ )
            {
                $sum = 0;
                for ( $j = 0; $j < $i; $j++ )
                {
                    $sum += abs($matrix->[0][$i][$j]);
                }
                for ( $j = ($i + 1); $j < $n; $j++ )
                {
                    $sum += abs($matrix->[0][$i][$j]);
                }
                $sum /= abs($matrix->[0][$i][$i]);
                if ($sum > $max) { $max = $sum; }
            }
            $ok = ($max < 1)
        }
    }
    if ($ok) { return($norm); }
    else     { return(0); }
}

sub solve_GSM  #  Global Step Method
{
    croak "Usage: \$xn_vector = \$matrix->solve_GSM(\$x0_vector,\$b_vector,\$epsilon);"
      if (@_ != 4);

    my($matrix,$x0_vector,$b_vector,$epsilon) = @_;
    my($rows1,$cols1) = (   $matrix->[1],   $matrix->[2]);
    my($rows2,$cols2) = ($x0_vector->[1],$x0_vector->[2]);
    my($rows3,$cols3) = ( $b_vector->[1], $b_vector->[2]);
    my($norm,$sum,$diff);
    my($xn_vector);
    my($i,$j,$n);

    croak "Math::MatrixReal::solve_GSM(): matrix is not quadratic"
      unless ($rows1 == $cols1);

    $n = $rows1;

    croak "Math::MatrixReal::solve_GSM(): 1st vector is not a column vector"
      unless ($cols2 == 1);

    croak "Math::MatrixReal::solve_GSM(): 2nd vector is not a column vector"
      unless ($cols3 == 1);

    croak "Math::MatrixReal::solve_GSM(): matrix and vector size mismatch"
      unless (($rows2 == $n) && ($rows3 == $n));

    return() unless ($norm = $matrix->_init_iteration());

    $xn_vector = $x0_vector->new($n,1);

    $diff = $epsilon + 1;
    while ($diff >= $epsilon)
    {
        for ( $i = 0; $i < $n; $i++ )
        {
            $sum = $b_vector->[0][$i][0];
            for ( $j = 0; $j < $i; $j++ )
            {
                $sum -= $matrix->[0][$i][$j] * $x0_vector->[0][$j][0];
            }
            for ( $j = ($i + 1); $j < $n; $j++ )
            {
                $sum -= $matrix->[0][$i][$j] * $x0_vector->[0][$j][0];
            }
            $xn_vector->[0][$i][0] = $sum / $matrix->[0][$i][$i];
        }
        $x0_vector->subtract($x0_vector,$xn_vector);
        if ($norm > 0) { $diff = $x0_vector->norm_one(); }
        else           { $diff = $x0_vector->norm_max(); }
        for ( $i = 0; $i < $n; $i++ )
        {
            $x0_vector->[0][$i][0] = $xn_vector->[0][$i][0];
        }
    }
    return($xn_vector);
}

sub solve_SSM  #  Single Step Method
{
    croak "Usage: \$xn_vector = \$matrix->solve_SSM(\$x0_vector,\$b_vector,\$epsilon);"
      if (@_ != 4);

    my($matrix,$x0_vector,$b_vector,$epsilon) = @_;
    my($rows1,$cols1) = (   $matrix->[1],   $matrix->[2]);
    my($rows2,$cols2) = ($x0_vector->[1],$x0_vector->[2]);
    my($rows3,$cols3) = ( $b_vector->[1], $b_vector->[2]);
    my($norm,$sum,$diff);
    my($xn_vector);
    my($i,$j,$n);

    croak "Math::MatrixReal::solve_SSM(): matrix is not quadratic"
      unless ($rows1 == $cols1);

    $n = $rows1;

    croak "Math::MatrixReal::solve_SSM(): 1st vector is not a column vector"
      unless ($cols2 == 1);

    croak "Math::MatrixReal::solve_SSM(): 2nd vector is not a column vector"
      unless ($cols3 == 1);

    croak "Math::MatrixReal::solve_SSM(): matrix and vector size mismatch"
      unless (($rows2 == $n) && ($rows3 == $n));

    return() unless ($norm = $matrix->_init_iteration());

    $xn_vector = $x0_vector->new($n,1);
    $xn_vector->copy($x0_vector);

    $diff = $epsilon + 1;
    while ($diff >= $epsilon)
    {
        for ( $i = 0; $i < $n; $i++ )
        {
            $sum = $b_vector->[0][$i][0];
            for ( $j = 0; $j < $i; $j++ )
            {
                $sum -= $matrix->[0][$i][$j] * $xn_vector->[0][$j][0];
            }
            for ( $j = ($i + 1); $j < $n; $j++ )
            {
                $sum -= $matrix->[0][$i][$j] * $xn_vector->[0][$j][0];
            }
            $xn_vector->[0][$i][0] = $sum / $matrix->[0][$i][$i];
        }
        $x0_vector->subtract($x0_vector,$xn_vector);
        if ($norm > 0) { $diff = $x0_vector->norm_one(); }
        else           { $diff = $x0_vector->norm_max(); }
        for ( $i = 0; $i < $n; $i++ )
        {
            $x0_vector->[0][$i][0] = $xn_vector->[0][$i][0];
        }
    }
    return($xn_vector);
}

sub solve_RM  #  Relaxation Method
{
    croak "Usage: \$xn_vector = \$matrix->solve_RM(\$x0_vector,\$b_vector,\$weight,\$epsilon);"
      if (@_ != 5);

    my($matrix,$x0_vector,$b_vector,$weight,$epsilon) = @_;
    my($rows1,$cols1) = (   $matrix->[1],   $matrix->[2]);
    my($rows2,$cols2) = ($x0_vector->[1],$x0_vector->[2]);
    my($rows3,$cols3) = ( $b_vector->[1], $b_vector->[2]);
    my($norm,$sum,$diff);
    my($xn_vector);
    my($i,$j,$n);

    croak "Math::MatrixReal::solve_RM(): matrix is not quadratic"
      unless ($rows1 == $cols1);

    $n = $rows1;

    croak "Math::MatrixReal::solve_RM(): 1st vector is not a column vector"
      unless ($cols2 == 1);

    croak "Math::MatrixReal::solve_RM(): 2nd vector is not a column vector"
      unless ($cols3 == 1);

    croak "Math::MatrixReal::solve_RM(): matrix and vector size mismatch"
      unless (($rows2 == $n) && ($rows3 == $n));

    return() unless ($norm = $matrix->_init_iteration());

    $xn_vector = $x0_vector->new($n,1);
    $xn_vector->copy($x0_vector);

    $diff = $epsilon + 1;
    while ($diff >= $epsilon)
    {
        for ( $i = 0; $i < $n; $i++ )
        {
            $sum = $b_vector->[0][$i][0];
            for ( $j = 0; $j < $i; $j++ )
            {
                $sum -= $matrix->[0][$i][$j] * $xn_vector->[0][$j][0];
            }
            for ( $j = ($i + 1); $j < $n; $j++ )
            {
                $sum -= $matrix->[0][$i][$j] * $xn_vector->[0][$j][0];
            }
            $xn_vector->[0][$i][0] = $weight * ( $sum / $matrix->[0][$i][$i] )
                                   + (1 - $weight) * $xn_vector->[0][$i][0];
        }
        $x0_vector->subtract($x0_vector,$xn_vector);
        if ($norm > 0) { $diff = $x0_vector->norm_one(); }
        else           { $diff = $x0_vector->norm_max(); }
        for ( $i = 0; $i < $n; $i++ )
        {
            $x0_vector->[0][$i][0] = $xn_vector->[0][$i][0];
        }
    }
    return($xn_vector);
}

# Core householder reduction routine (when eigenvector
# are wanted).
# Adapted from: Numerical Recipes, 2nd edition.
sub _householder_vectors ($)
{
    my ($Q) = @_;
    my ($rows, $cols) = ($Q->[1], $Q->[2]);

    # Creates tridiagonal
    # Set up tridiagonal needed elements
    my $d = []; # N Diagonal elements 0...n-1
    my $e = []; # N-1 Off-Diagonal elements 0...n-2

    my @p = ();
    for (my $i = ($rows-1); $i > 1; $i--)
    {
    my $scale = 0.0;
    # Computes norm of one column (below diagonal)
    for (my $k = 0; $k < $i; $k++)
    {
        $scale += abs($Q->[0][$i][$k]);
    }
    if ($scale == 0.0)
    { # skip the transformation
        $e->[$i-1] = $Q->[0][$i][$i-1];
    }
    else
    {
        my $h = 0.0;
        for (my $k = 0; $k < $i; $k++)
        { # Used scaled Q for transformation
            $Q->[0][$i][$k] /= $scale;
            # Form sigma in h
            $h += $Q->[0][$i][$k] * $Q->[0][$i][$k];
        }
        my $t1 = $Q->[0][$i][$i-1];
        my $t2 = (($t1 >= 0.0) ? -sqrt($h) : sqrt($h));
        $e->[$i-1] = $scale * $t2; # Update off-diagonals
        $h -= $t1 * $t2;
        $Q->[0][$i][$i-1] -= $t2;
        my $f = 0.0;
        for (my $j = 0; $j < $i; $j++)
        {
            $Q->[0][$j][$i] = $Q->[0][$i][$j] / $h;
            my $g = 0.0;
            for (my $k = 0; $k <= $j; $k++)
            {
                $g += $Q->[0][$j][$k] * $Q->[0][$i][$k];
            }
            for (my $k = $j+1; $k < $i; $k++)
            {
                $g += $Q->[0][$k][$j] * $Q->[0][$i][$k];
            }
            # Form elements of P
            $p[$j] = $g / $h;
            $f += $p[$j] * $Q->[0][$i][$j];
        }
        my $hh = $f / ($h + $h);
        for (my $j = 0; $j < $i; $j++)
        {
            my $t3 = $Q->[0][$i][$j];
            my $t4 = $p[$j] - $hh * $t3;
            $p[$j] = $t4;
            for (my $k = 0; $k <= $j; $k++)
            {
                $Q->[0][$j][$k] -= $t3 * $p[$k]
                + $t4 * $Q->[0][$i][$k];
            }
        }
    }
    }
    # Updates for i == 0,1
    $e->[0] = $Q->[0][1][0];
    $d->[0] = $Q->[0][0][0]; # i==0
    $Q->[0][0][0] = 1.0;
    $d->[1] = $Q->[0][1][1]; # i==1
    $Q->[0][1][1] = 1.0;
    $Q->[0][1][0] = $Q->[0][0][1] = 0.0;
    for (my $i = 2; $i < $rows; $i++)
    {
        for (my $j = 0; $j < $i; $j++)
        {
            my $g = 0.0;
            for (my $k = 0; $k < $i; $k++)
            {
                $g += $Q->[0][$i][$k] * $Q->[0][$k][$j];
            }
            for (my $k = 0; $k < $i; $k++)
            {
                $Q->[0][$k][$j] -= $g * $Q->[0][$k][$i];
            }
        }
        $d->[$i] = $Q->[0][$i][$i];
        # Reset row and column of Q for next iteration
        $Q->[0][$i][$i] = 1.0;
        for (my $j = 0; $j < $i; $j++)
        {
            $Q->[0][$i][$j] = $Q->[0][$j][$i] = 0.0;
        }
    }
    return ($d, $e);
}

# Computes sqrt(a*a + b*b), but more carefully...
sub _pythag ($$)
{
    my ($a, $b) = @_;
    my $aa = abs($a);
    my $ab = abs($b);
    if ($aa > $ab)
    {
        # NB: Not needed!: return 0.0 if ($aa == 0.0);
        my $t = $ab / $aa;
        return ($aa * sqrt(1.0 + $t*$t));
    } else {
        return 0.0 if ($ab == 0.0);
        my $t = $aa / $ab;
        return ($ab * sqrt(1.0 + $t*$t));
    }
}

# QL algorithm with implicit shifts to determine the eigenvalues
# of a tridiagonal matrix. Internal routine.
sub _tridiagonal_QLimplicit
{
    my ($EV, $d, $e) = @_;
    my ($rows, $cols) = ($EV->[1], $EV->[2]);

    $e->[$rows-1] = 0.0;
    # Start real computation
    for (my $l = 0; $l < $rows; $l++)
    {
        my $iter = 0;
        my $m;
        OUTER:
        do {
            for ($m = $l; $m < ($rows - 1); $m++)
            {
                my $dd = abs($d->[$m]) + abs($d->[$m+1]);
                last if ((abs($e->[$m]) + $dd) == $dd);
            }
            if ($m != $l)
            {
                ## why only allow 30 iterations?
                croak("Too many iterations!") if ($iter++ >= 30);
                my $g = ($d->[$l+1] - $d->[$l])
                    / (2.0 * $e->[$l]);
                my $r = _pythag($g, 1.0);
                $g = $d->[$m] - $d->[$l]
                    + $e->[$l] / ($g + (($g >= 0.0) ? abs($r) : -abs($r)));
                my ($p,$s,$c) = (0.0, 1.0,1.0);
            for (my $i = ($m-1); $i >= $l; $i--)
            {
                my $ii = $i + 1;
                my $f = $s * $e->[$i];
                my $t = _pythag($f, $g);
                $e->[$ii] = $t;
                if ($t == 0.0)
                {
                    $d->[$ii] -= $p;
                    $e->[$m] = 0.0;
                    next OUTER;
                }
                my $b = $c * $e->[$i];
                $s = $f / $t;
                $c = $g / $t;
                $g = $d->[$ii] - $p;
                my $t2 = ($d->[$i] - $g) * $s + 2.0 * $c * $b;
                $p = $s * $t2;
                $d->[$ii] = $g + $p;
                $g = $c * $t2 - $b;
                for (my $k = 0; $k < $rows; $k++)
                {
                    my $t1 = $EV->[0][$k][$ii];
                    my $t2 = $EV->[0][$k][$i];
                    $EV->[0][$k][$ii] = $s * $t2 + $c * $t1;
                    $EV->[0][$k][$i] = $c * $t2 - $s * $t1;
                }
            }
            $d->[$l] -= $p;
            $e->[$l] = $g;
            $e->[$m] = 0.0;
            }
        } while ($m != $l);
    }
    return;
}

# Core householder reduction routine (when eigenvector
# are NOT wanted).
sub _householder_values ($)
{
    my ($Q) = @_; # NB: Q is destroyed on output...
    my ($rows, $cols) = ($Q->[1], $Q->[2]);

    # Creates tridiagonal
    # Set up tridiagonal needed elements
    my $d = []; # N Diagonal elements 0...n-1
    my $e = []; # N-1 Off-Diagonal elements 0...n-2

    my @p = ();
    for (my $i = ($rows - 1); $i > 1; $i--)
    {
        my $scale = 0.0;
        for (my $k = 0; $k < $i; $k++)
        {
            $scale += abs($Q->[0][$i][$k]);
        }
        if ($scale == 0.0)
        { # skip the transformation
            $e->[$i-1] = $Q->[0][$i][$i-1];
        }
        else
        {
            my $h = 0.0;
            for (my $k = 0; $k < $i; $k++)
            { # Used scaled Q for transformation
                $Q->[0][$i][$k] /= $scale;
                # Form sigma in h
                $h += $Q->[0][$i][$k] * $Q->[0][$i][$k];
            }
            my $t = $Q->[0][$i][$i-1];
            my $t2 = (($t >= 0.0) ? -sqrt($h) : sqrt($h));
            $e->[$i-1] = $scale * $t2; # Updates off-diagonal
            $h -= $t * $t2;
            $Q->[0][$i][$i-1] -= $t2;
            my $f = 0.0;
            for (my $j = 0; $j < $i; $j++)
            {
                my $g = 0.0;
                for (my $k = 0; $k <= $j; $k++)
                {
                    $g += $Q->[0][$j][$k] * $Q->[0][$i][$k];
                }
                for (my $k = $j+1; $k < $i; $k++)
                {
                    $g += $Q->[0][$k][$j] * $Q->[0][$i][$k];
                }
                # Form elements of P
                $p[$j] = $g / $h;
                $f += $p[$j] * $Q->[0][$i][$j];
            }
            my $hh = $f / ($h + $h);
            for (my $j = 0; $j < $i; $j++)
            {
                my $t = $Q->[0][$i][$j];
                my $g = $p[$j] - $hh * $t;
                $p[$j] = $g;
                for (my $k = 0; $k <= $j; $k++)
                {
                    $Q->[0][$j][$k] -= $t * $p[$k]
                    + $g * $Q->[0][$i][$k];
                }
            }
        }
    }
    # Updates for i==1
    $e->[0] =  $Q->[0][1][0];
    # Updates diagonal elements
    for (my $i = 0; $i < $rows; $i++)
    {
        $d->[$i] =  $Q->[0][$i][$i];
    }
    return ($d, $e);
}

# QL algorithm with implicit shifts to determine the
# eigenvalues ONLY. This is O(N^2) only...
sub _tridiagonal_QLimplicit_values
{
    my ($M, $d, $e) = @_; # NB: M is not touched...
    my ($rows, $cols) = ($M->[1], $M->[2]);

    $e->[$rows-1] = 0.0;
    # Start real computation
    for (my $l = 0; $l < $rows; $l++)
    {
        my $iter = 0;
        my $m;
        OUTER:
        do {
            for ($m = $l; $m < ($rows - 1); $m++)
            {
                my $dd = abs($d->[$m]) + abs($d->[$m+1]);
                last if ((abs($e->[$m]) + $dd) == $dd);
            }
            if ($m != $l)
            {
                croak("Too many iterations!") if ($iter++ >= 30);
                my $g = ($d->[$l+1] - $d->[$l])
                    / (2.0 * $e->[$l]);
                my $r = _pythag($g, 1.0);
                $g = $d->[$m] - $d->[$l]
                    + $e->[$l] / ($g + (($g >= 0.0) ? abs($r) : -abs($r)));
                my ($p,$s,$c) = (0.0, 1.0,1.0);
                for (my $i = ($m-1); $i >= $l; $i--)
                {
                    my $ii = $i + 1;
                    my $f = $s * $e->[$i];
                    my $t = _pythag($f, $g);
                    $e->[$ii] = $t;
                    if ($t == 0.0)
                    {
                        $d->[$ii] -= $p;
                        $e->[$m] = 0.0;
                        next OUTER;
                    }
                    my $b = $c * $e->[$i];
                    $s = $f / $t;
                    $c = $g / $t;
                    $g = $d->[$ii] - $p;
                    my $t2 = ($d->[$i] - $g) * $s + 2.0 * $c * $b;
                    $p = $s * $t2;
                    $d->[$ii] = $g + $p;
                    $g = $c * $t2 - $b;
                }
                $d->[$l] -= $p;
                $e->[$l] = $g;
                $e->[$m] = 0.0;
            }
        } while ($m != $l);
    }
    return;
}

# Householder reduction of a real, symmetric matrix A.
# Returns a tridiagonal matrix T and the orthogonal matrix
# Q effecting the transformation between A and T.
sub householder ($)
{
    my ($A) = @_;
    my ($rows, $cols) = ($A->[1], $A->[2]);

    croak "Matrix is not quadratic"
        unless ($rows = $cols);
    croak "Matrix is not symmetric"
        unless ($A->is_symmetric());

    # Copy given matrix TODO: study if we should do in-place modification
    my $Q = $A->clone();

    # Do the computation of tridiagonal elements and of
    # transformation matrix
    my ($diag, $offdiag) = $Q->_householder_vectors();

    # Creates the tridiagonal matrix
    my $T = $A->shadow();
    for (my $i = 0; $i < $rows; $i++)
    { # Set diagonal
        $T->[0][$i][$i] = $diag->[$i];
    }
    for (my $i = 0; $i < ($rows-1); $i++)
    { # Set off diagonals
        $T->[0][$i+1][$i] = $offdiag->[$i];
        $T->[0][$i][$i+1] = $offdiag->[$i];
    }
    return ($T, $Q);
}

# QL algorithm with implicit shifts to determine the eigenvalues
# and eigenvectors of a real tridiagonal matrix - or of a matrix
# previously reduced to tridiagonal form.
sub tri_diagonalize ($;$)
{
    my ($T,$Q) = @_; # Q may be 0 if the original matrix is really tridiagonal

    my ($rows, $cols) = ($T->[1], $T->[2]);

    croak "Matrix is not quadratic"
        unless ($rows = $cols);
    croak "Matrix is not tridiagonal"
        unless ($T->is_tridiagonal()); # DONE

    my $EV;
    # Obtain/Creates the todo eigenvectors matrix
    if ($Q)
    {
        $EV = $Q->clone();
    }
    else
    {
        $EV = $T->shadow();
        $EV->one();
    }
    # Allocates diagonal vector
    my $diag = [ ];
    # Initializes it with T
    for (my $i = 0; $i < $rows; $i++)
    {
        $diag->[$i] = $T->[0][$i][$i];
    }
    # Allocate temporary vector for off-diagonal elements
    my $offdiag = [ ];
    for (my $i = 1; $i < $rows; $i++)
    {
        $offdiag->[$i-1] = $T->[0][$i][$i-1];
    }

    # Calls the calculus routine
    $EV->_tridiagonal_QLimplicit($diag, $offdiag);

    # Allocate eigenvalues vector
    my $v = Math::MatrixReal->new($rows,1);
    # Fills it
    for (my $i = 0; $i < $rows; $i++)
    {
        $v->[0][$i][0] = $diag->[$i];
    }
    return ($v, $EV);
}

# Main routine for diagonalization of a real symmetric
# matrix M. Operates by transforming M into a tridiagonal
# matrix and then obtaining the eigenvalues and eigenvectors
# for that matrix (taking into account the transformation to
# tridiagonal).
sub sym_diagonalize ($)
{
    my ($M) = @_;
    my ($rows, $cols) = ($M->[1], $M->[2]);

    croak "Matrix is not quadratic"
        unless ($rows = $cols);
    croak "Matrix is not symmetric"
        unless ($M->is_symmetric());

    # Copy initial matrix
    # TODO: study if we should allow in-place modification
    my $VEC = $M->clone();

    # Do the computation of tridiagonal elements and of
    # transformation matrix
    my ($diag, $offdiag) = $VEC->_householder_vectors();
    # Calls the calculus routine for diagonalization
    $VEC->_tridiagonal_QLimplicit($diag, $offdiag);

    # Allocate eigenvalues vector
    my $val = Math::MatrixReal->new($rows,1);
    # Fills it
    for (my $i = 0; $i < $rows; $i++)
    {
        $val->[0][$i][0] = $diag->[$i];
    }
    return ($val, $VEC);
}

# Householder reduction of a real, symmetric matrix A.
# Returns a tridiagonal matrix T equivalent to A.
sub householder_tridiagonal ($)
{
    my ($A) = @_;
    my ($rows, $cols) = ($A->[1], $A->[2]);

    croak "Matrix is not quadratic"
        unless ($rows = $cols);
    croak "Matrix is not symmetric"
        unless ($A->is_symmetric());

    # Copy given matrix
    my $Q = $A->clone();

    # Do the computation of tridiagonal elements and of
    # transformation matrix
    # Q is destroyed after reduction
    my ($diag, $offdiag) = $Q->_householder_values();

    # Creates the tridiagonal matrix in Q (avoid allocation)
    my $T = $Q;
    $T->zero();
    for (my $i = 0; $i < $rows; $i++)
    { # Set diagonal
        $T->[0][$i][$i] = $diag->[$i];
    }
    for (my $i = 0; $i < ($rows-1); $i++)
    { # Set off diagonals
        $T->[0][$i+1][$i] = $offdiag->[$i];
        $T->[0][$i][$i+1] = $offdiag->[$i];
    }
    return $T;
}

# QL algorithm with implicit shifts to determine ONLY
# the eigenvalues a real tridiagonal matrix - or of a
# matrix previously reduced to tridiagonal form.
sub tri_eigenvalues ($;$)
{
    my ($T) = @_;
    my ($rows, $cols) = ($T->[1], $T->[2]);

    croak "Matrix is not quadratic"
        unless ($rows = $cols);
    croak "Matrix is not tridiagonal"
        unless ($T->is_tridiagonal() ); # DONE

    # Allocates diagonal vector
    my $diag = [ ];
    # Initializes it with T
    for (my $i = 0; $i < $rows; $i++)
    {
        $diag->[$i] = $T->[0][$i][$i];
    }
    # Allocate temporary vector for off-diagonal elements
    my $offdiag = [ ];
    for (my $i = 1; $i < $rows; $i++)
    {
        $offdiag->[$i-1] = $T->[0][$i][$i-1];
    }

    # Calls the calculus routine (T is not touched)
    $T->_tridiagonal_QLimplicit_values($diag, $offdiag);

    # Allocate eigenvalues vector
    my $v = Math::MatrixReal->new($rows,1);
    # Fills it
    for (my $i = 0; $i < $rows; $i++)
    {
        $v->[0][$i][0] = $diag->[$i];
    }
    return $v;
}

## more general routine than sym_eigenvalues
sub eigenvalues ($){
    my ($matrix) = @_;
    my ($rows,$cols) = $matrix->dim();

    croak "Matrix is not quadratic" unless ($rows == $cols);

    if($matrix->is_upper_triangular() || $matrix->is_lower_triangular() ){
        my $l = Math::MatrixReal->new($rows,1);
        map { $l->[0][$_][0] = $matrix->[0][$_][$_] } (0 .. $rows-1);
        return $l;
    }

    return sym_eigenvalues($matrix) if $matrix->is_symmetric();

    carp "Math::MatrixReal::eigenvalues(): Matrix is not symmetric or triangular";
    return undef;

}
# Main routine for diagonalization of a real symmetric
# matrix M. Operates by transforming M into a tridiagonal
# matrix and then obtaining the eigenvalues and eigenvectors
# for that matrix (taking into account the transformation to
# tridiagonal).
sub sym_eigenvalues ($)
{
    my ($M) = @_;
    my ($rows, $cols) = ($M->[1], $M->[2]);

    croak "Matrix is not quadratic" unless ($rows == $cols);
    croak "Matrix is not symmetric" unless ($M->is_symmetric);

    # Copy matrix in temporary
    my $A = $M->clone();
    # Do the computation of tridiagonal elements and of
    # transformation matrix. A is destroyed
    my ($diag, $offdiag) = $A->_householder_values();
    # Calls the calculus routine for diagonalization
    # (M is not touched)
    $M->_tridiagonal_QLimplicit_values($diag, $offdiag);

    # Allocate eigenvalues vector
    my $val = Math::MatrixReal->new($rows,1);
    # Fills it
    map { $val->[0][$_][0] = $diag->[$_] } ( 0 .. $rows-1);

    return $val;
}
#TODO: docs+test
sub is_positive_definite {
    my ($matrix) = @_;
    my ($r,$c) = $matrix->dim;

    croak "Math::MatrixReal::is_positive_definite(): Matrix is not square" unless ($r == $c);
    # must have positive (i.e REAL) eigenvalues to be positive definite
    return 0 unless $matrix->is_symmetric;

    my $ev = $matrix->eigenvalues;
    my $pos = 1;
    $ev->each(sub { my $x = shift; if ($x <= 0){ $pos=0;return; } } );
    return $pos;
}
#TODO: docs+test
sub is_positive_semidefinite {
    my ($matrix) = @_;
    my ($r,$c) = $matrix->dim;

    croak "Math::MatrixReal::is_positive_semidefinite(): Matrix is not square" unless ($r == $c);
    # must have nonnegative (i.e REAL) eigenvalues to be positive semidefinite
    return 0 unless $matrix->is_symmetric;

    my $ev = $matrix->eigenvalues;
    my $pos = 1;
    $ev->each(sub { my $x = shift; if ($x < 0){ $pos=0;return; } } );
    return $pos;
}
sub is_row { return (shift)->is_row_vector }
sub is_col { return (shift)->is_col_vector }

sub is_row_vector {
    my ($m) = @_;
    my $r = $m->[1];
    $r == 1 ? 1 : 0;
}
sub is_col_vector {
    my ($m) = @_;
    my $c = $m->[2];
    $c == 1 ? 1 : 0;
}

sub is_orthogonal($) {
    my ($matrix) = @_;

    return 0 unless $matrix->is_quadratic;

    my $one = $matrix->shadow();
    $one->one;

    abs(~$matrix * $matrix - $one) < 1e-12 ? return 1 : return 0;

}

sub is_positive($) {
    my ($m) = @_;
    my $pos = 1;
    $m->each( sub { if( (shift) <= 0){ $pos = 0;return;} } );
    return $pos;
}
sub is_negative($) {
    my ($m) = @_;
    my $neg = 1;
    $m->each( sub { if( (shift) >= 0){ $neg = 0;return;} } );
    return $neg;
}


sub is_periodic($$) {
    my ($m,$k) = @_;
        return 0 unless $m->is_quadratic();
    abs($m**(int($k)+1) - $m) < 1e-12 ? return 1 : return 0;
}
sub is_idempotent($) {
    return (shift)->is_periodic(1);
}

# Boolean check routine to see if a matrix is
# symmetric
sub is_symmetric ($)
{
  my ($M) = @_;
  my ($rows, $cols) = ($M->[1], $M->[2]);
  # if it is not quadratic it cannot be symmetric...
  return 0 unless ($rows == $cols);
  # skip when $i=$j?
  for (my $i = 1; $i < $rows; $i++)
    {
      for (my $j = 0; $j < $i; $j++)
        {
          return 0 unless ($M->[0][$i][$j] == $M->[0][$j][$i]);
        }
    }
  return 1;
}
# Boolean check to see if matrix is tridiagonal
sub is_tridiagonal ($) {
    my ($M) = @_;
    my ($rows,$cols) = ($M->[1],$M->[2]);
    my ($i,$j) = (0,0);
    # if it is not quadratic it cannot be tridiag
    return 0 unless ($rows == $cols);

    for(;$i < $rows; $i++ ){
        for(;$j < $cols; $j++ ){
            #print "debug: testing $i,$j = " . $M->[0][$i][$j] . "\n";
            # skip diag and diag+-1
            next if ($i == $j);
            next if ($i+1 == $j);
            next if ($i-1 == $j);
            return 0 if $M->[0][$i][$j];
        }
        $j = 0;
    }
    return 1;
}
# Boolean check to see if matrix is upper triangular
# i.e all nonzero elements are above main diagonal
sub is_upper_triangular {
    my ($M) = @_;
    my ($rows,$cols) = $M->dim();
    my ($i,$j) = (1,0);
    return 0 unless ($rows == $cols);

    for(;$i < $rows; $i++ ){
        for(;$j < $cols;$j++ ){
            next if ($i <= $j);
            return 0 if $M->[0][$i][$j];
        }
        $j = 0;
    }
    return 1;
}
# Boolean check to see if matrix is lower triangular
# i.e all nonzero elements are lower main diagonal
sub is_lower_triangular {
    my ($M) = @_;
    my ($rows,$cols) = $M->dim();
    my ($i,$j) = (0,1);
    return 0 unless ($rows == $cols);

    for(;$i < $rows; $i++ ){
        for(;$j < $cols;$j++ ){
            next if ($i >= $j);
            return 0 if $M->[0][$i][$j];
        }
        $j = 0;
    }
    return 1;
}

# Boolean check to see if matrix is diagonal
sub is_diagonal ($) {
    my ($M) = @_;
    my ($rows,$cols) = ($M->[1],$M->[2]);
    my ($i,$j) = (0,0);
    return 0 unless ($rows == $cols );
    for(;$i < $rows; $i++ ){
        for(;$j < $cols; $j++ ){
            # skip diag elements
            next if ($i == $j);
            return 0 if $M->[0][$i][$j];
        }
        $j = 0;
    }
    return 1;
}
sub is_quadratic ($) {
    croak "Usage: \$matrix->is_quadratic()" unless (@_ == 1);
    my ($matrix) = @_;
    $matrix->[1] == $matrix->[2] ? return 1 : return 0;
}

sub is_square($) {
    croak "Usage: \$matrix->is_square()" unless (@_ == 1);
    return (shift)->is_quadratic();
}

sub is_LR($) {
    croak "Usage: \$matrix->is_LR()" unless (@_ == 1);
    return (shift)->[3] ? 1 : 0;
}

sub is_normal{
    my ($matrix,$eps) = @_;
    my ($rows,$cols) = $matrix->dim;

    $eps ||= 1e-8;

    (~$matrix * $matrix - $matrix * ~$matrix < $eps ) ? 1 : 0;

}

sub is_skew_symmetric{
    my ($m) = @_;
    my ($rows, $cols) = $m->dim;
    # if it is not quadratic it cannot be skew symmetric...
    return 0 unless ($rows == $cols);
    for (my $i = 1; $i < $rows; $i++) {
        for (my $j = 0; $j < $i; $j++) {
            return 0 unless ($m->[0][$i][$j] == -$m->[0][$j][$i]);
        }
    }
    return 1;

}
####
sub is_gramian{
    my ($m) = @_;
    my ($rows,$cols) = $m->dim;
    my $neg=0;
    # gramian matrix must be symmetric
    return 0 unless $m->is_symmetric;

    # must have all non-negative eigenvalues
    my $ev = $m->eigenvalues;
    $ev->each(sub { $neg++ if ((shift)<0) } );

    return $neg ? 0 : 1;
}
sub is_binary{
    my ($m) = @_;
    my ($rows, $cols) = $m->dim;
    for (my $i = 0; $i < $rows; $i++) {
        for (my $j = 0; $j < $cols; $j++) {
            return 0 unless ($m->[0][$i][$j] == 1 || $m->[0][$i][$j] == 0);
        }
    }
    return 1;

}
sub as_scilab {
    return (shift)->as_matlab;
}

sub as_matlab {
    my ($m) = shift;
    my %args = (
        format => "%s",
        name => "",
        semi => 0,
        @_);
    my ($row,$col) = $m->dim;
    my $s = "";

    if( $args{name} ){
        $s = "$args{name} = ";
    }
    $s .= "[";
    $m->each(
        sub { my($x,$i,$j) = @_;
            $s .= sprintf(" $args{format}",$x);
            $s .= ";\n" if( $j == $col && $i != $row);
        }
        );
    $s .= "]";
    $s .= ";" if $args{semi};
    return $s;
}
#TODO: docs+test
sub as_yacas{
    my ($m) = shift;
    my %args = (
        format => "%s",
        name => "",
        semi => 0,
        @_);
    my ($row,$col) = $m->dim;
    my $s = "";

    if( $args{name} ){
        $s = "$args{name} := ";
    }
    $s .= "{";

    $m->each(
        sub { my($x,$i,$j) = @_;
            $s .= "{" if ($j == 1);
                        $s .= sprintf("$args{format}",$x);
                        $s .= "," if( $j != $col );
            $s .= "}," if ($j == $col && $i != $row);
        }
    );
    $s .= "}}";

    return $s;
}

sub as_latex{
    my ($m) = shift;
    my %args = (
        format       => "%s",
        name         => "",
        align        => "c",
        display_math => 0,
    @_);


    my ($row,$col) = $m->dim;
    my $inside;
    my $s = <
\\left( \\begin{array}{%COLS%}
%INSIDE%\\end{array} \\right)
LATEX
    $args{align} = lc $args{align};
    if( $args{align} !~ m/^(c|l|r)$/ ){
        croak "Math::MatrixReal::as_latex(): Invalid alignment '$args{align}'";
    }

    $s =~ s/%COLS%/$args{align} x $col/em;

    if( $args{name} ){
        $s = "$args{name} = $s";
    }
    $m->each(
        sub {
            my ($x,$i,$j) = @_;
            $x = sprintf($args{format},$x);

            # last element in each row gets a \\
            if ($j == $col && $i != $row){
                $inside .= "$x \\\\"."\n";
            # the annoying last line has neither
            } elsif( $j == $col && $i == $row){
                $inside .= "$x\n";
            } else {
                $inside .= "$x&";
            }
        }
    );
    if($args{displaymath}){
            $s = "\\[$s\\]";
    } else {
            $s = "\$$s\$";
    }
    $s =~ s/%INSIDE%/$inside/gm;
    return $s;
}
####
sub spectral_radius
{
    my ($matrix) = @_;
    my ($r,$c) = $matrix->dim;
    my $ev = $matrix->eigenvalues;
    my $radius=0;
    $ev->each(sub { my $x = shift; $radius = $x if (abs($x) > $radius); } );
    return $radius;
}

sub maximum {
	my ($matrix) = @_;
	my ($rows, $columns) = $matrix->dim;

	my $max = [];
	my $max_p = [];

	if ($rows == 1) {
		($max, $max_p) = _max_column($matrix->row(1)->_transpose, $columns);
	} elsif ($columns == 1) {
		($max, $max_p) = _max_column($matrix->column(1), $rows);
	} else {
		for my $c (1..$columns) {
			my ($m, $mp) = _max_column($matrix->column($c), $rows);
			push @$max, $m;
			push @$max_p, $mp;
		}
	}
	return wantarray ? ($max, $max_p) : $max
}

sub _max_column {
	# passing $rows allows for some extra (minimal) efficiency
	my ($column, $rows) = @_;

	my ($m, $mp) = ($column->element(1, 1), 1);
	for my $l (1..$rows) {
		if ($column->element($l, 1) > $m) {
			$m = $column->element($l, 1);
			$mp = $l;
		}
	}
	return ($m, $mp);
}

sub minimum {
	my ($matrix) = @_;
	my ($rows, $columns) = $matrix->dim;

	my $min = [];
	my $min_p = [];

	if ($rows == 1) {
		($min, $min_p) = _min_column($matrix->row(1)->_transpose, $columns);
	} elsif ($columns == 1) {
		($min, $min_p) = _min_column($matrix->column(1), $rows);
	} else {
		for my $c (1..$columns) {
			my ($m, $mp) = _min_column($matrix->column($c), $rows);
			push @$min, $m;
			push @$min_p, $mp;
		}
	}
	return wantarray ? ($min, $min_p) : $min
}

sub _min_column {
	# passing $rows allows for some extra (minimal) efficiency
	my ($column, $rows) = @_;

	my ($m, $mp) = ($column->element(1, 1), 1);
	for my $l (1..$rows) {
		if ($column->element($l, 1) < $m) {
			$m = $column->element($l, 1);
			$mp = $l;
		}
	}
	return ($m, $mp);
}



                ########################################
                #                                      #
                # define overloaded operators section: #
                #                                      #
                ########################################
sub _concat
{
    my($object,$argument,$flag) = @_;
    my($orows,$ocols) 		= ($object->[1],$object->[2]);
    my($name)			= "concat";


    if ((defined $argument) && ref($argument) && (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/)) {
    	my($arows,$acols) 		= ($argument->[1],$argument->[2]);
        croak "Math::MatrixReal: Matrices must have same number of rows in concatenation" unless ($orows == $arows);
     	my $result = $object->new($orows,$ocols+$acols);
        for ( my $i = 0; $i < $arows; $i++ ) {
            for ( my $j = 0; $j < $ocols + $acols; $j++ ) {
		$result->[0][$i][$j] = ( $j <  $ocols ) ? $object->[0][$i][$j] : $argument->[0][$i][$j - $ocols] ;
            }
        }
        return $result;
    } elsif (defined $argument) {
	return "$object" . $argument;

    } else {
        croak "Math::MatrixReal $name: wrong argument type";
    }
}
sub _negate
{
    my($object) = @_;

    my $temp = $object->new($object->[1],$object->[2]);
    $temp->negate($object);
    return($temp);
}

sub _transpose
{
    my ($object) = @_;
    my $temp = $object->new($object->[2],$object->[1]);
    $temp->transpose($object);
    return $temp;
}

sub _boolean
{
    my($object) = @_;
    my($rows,$cols) = ($object->[1],$object->[2]);

    my $result = 0;

    BOOL:
    for ( my $i = 0; $i < $rows; $i++ )
    {
        for ( my $j = 0; $j < $cols; $j++ )
        {
            if ($object->[0][$i][$j] != 0)
            {
                $result = 1;
                last BOOL;
            }
        }
    }
    return($result);
}
#TODO: ugly copy+paste
sub _not_boolean
{
    my ($object) = @_;
    my ($rows,$cols) = ($object->[1],$object->[2]);

    my $result = 1;
    NOTBOOL:
    for ( my $i = 0; $i < $rows; $i++ )
    {
        for ( my $j = 0; $j < $cols; $j++ )
        {
            if ($object->[0][$i][$j] != 0)
            {
                $result = 0;
                last NOTBOOL;
            }
        }
    }
    return($result);
}

sub _stringify
{
    my ($self) = @_;
    my ($rows,$cols) = ($self->[1],$self->[2]);

    my $precision = $self->[4];

    my $format = !defined $precision ? '% #-19.12E ' : '% #-19.'.$precision.'f ';
    $format = '% #-12d' if defined $precision && $precision == 0;

    my $s = '';
    for ( my $i = 0; $i < $rows; $i++ )
    {
        $s .= "[ ";
        for ( my $j = 0; $j < $cols; $j++ )
        {
            $s .= sprintf $format , $self->[0][$i][$j];
        }
        $s .= "]\n";
    }
    return $s;
}

sub _norm
{
    my ($self) = @_;

    return $self->norm_one() ;
}

sub _add
{
    my($object,$argument,$flag) = @_;
    my($name) = "'+'";

    if ((defined $argument) && ref($argument) &&
        (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
    {
        if (defined $flag)
        {
            my $temp = $object->new($object->[1],$object->[2]);
            $temp->add($object,$argument);
            return($temp);
        }
        else
        {
            $object->add($object,$argument);
            return($object);
        }
    }
    else
    {
        croak "Math::MatrixReal $name: wrong argument type";
    }
}

sub _subtract
{
    my($object,$argument,$flag) = @_;
    my($name) = "'-'";

    if ((defined $argument) && ref($argument) &&
        (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
    {
        if (defined $flag)
        {
            my $temp = $object->new($object->[1],$object->[2]);
            if ($flag) { $temp->subtract($argument,$object); }
            else       { $temp->subtract($object,$argument); }
            return $temp;
        }
        else
        {
            $object->subtract($object,$argument);
            return($object);
        }
    }
    else
    {
        croak "Math::MatrixReal $name: wrong argument type";
    }
}

sub _exponent
{
    my($matrix, $exp) = @_;
    my($rows,$cols) = ($matrix->[1],$matrix->[2]);

    return $matrix->exponent( $exp );
}
sub _divide
{
	my($matrix,$argument,$flag) = @_;
	# TODO: check dimensions of everything!
	my($mrows,$mcols) = ($matrix->[1],$matrix->[2]);
	my($arows,$acols)=(0,0);
	my($name) = "'/'";
	my $temp = $matrix->clone();
	my $arg;
	my ($inv,$m1);

	if( ref($argument) =~ /Math::MatrixReal/ ){
		$arg =  $argument->clone();
		($arows,$acols)=($arg->[1],$arg->[2]);
	}
	#print "DEBUG: flag= $flag\n";
	#print "DEBUG: arg=$arg\n";
	if( $flag == 1) {
		  #print "DEBUG: ref(arg)= " . ref($arg) . "\n";
		if( ref($argument) =~ /Math::MatrixReal/ ){
			#print "DEBUG: arg is a matrix \n";
			# Matrix Division = A/B = A*B^(-1)
			croak "Math::MatrixReal $name: this operation is defined only for square matrices" unless ($arows == $acols);
			return $temp->multiply(  $arg->inverse() );

		} else {
			#print "DEBUG: Arg is scalar\n";
			#print "DEBUG:arows,acols=$arows,$acols\n";
			#print "DEBGU:mrows,mcols=$mrows,$mcols\n";
			 croak "Math::MatrixReal $name: this operation is defined only for square matrices" unless ($mrows == $mcols);
			$temp->multiply_scalar( $temp , $argument);
			return $temp;
		}
        } else {
	#print "DEBUG: temp=\n";
	#print $temp . "\n";
	#print "DEBUG: ref(arg)= " . ref($arg) . "\n";
	#print "DEBUG: arg=\n";
	#print $arg ."\n";
	if( ref($arg) =~ /Math::MatrixReal/ ){
		#print "DEBUG: matrix division\n";
		if( $arg->is_col_vector() ){
			print "DEBUG: $arg is a col vector\n";
		}
		croak "Math::MatrixReal $name: this operation is defined only for square matrices" unless ($arows == $acols);
		$inv = $arg->inverse();
		return $temp->multiply($inv);
	} else {
		$temp->multiply_scalar($temp,1/$argument);
		return $temp;
	}
    }

}

sub _multiply
{
    my($object,$argument,$flag) = @_;
    my($name) = "'*'";
    my($temp);

    if ((defined $argument) && ref($argument) &&
        (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
    {
        if ((defined $flag) && $flag)
        {
            return( multiply($argument,$object) );
        }
        else
        {
            return( multiply($object,$argument) );
        }
    }
    elsif ((defined $argument) && !(ref($argument)))
    {
        if (defined $flag)
        {
            $temp = $object->new($object->[1],$object->[2]);
            $temp->multiply_scalar($object,$argument);
            return($temp);
        }
        else
        {
            $object->multiply_scalar($object,$argument);
            return($object);
        }
    }
    else
    {
        croak "Math::MatrixReal $name: wrong argument type";
    }
}

sub _assign_add
{
    my($object,$argument) = @_;

    return( &_add($object,$argument,undef) );
}

sub _assign_subtract
{
    my($object,$argument) = @_;

    return( &_subtract($object,$argument,undef) );
}

sub _assign_multiply
{
    my($object,$argument) = @_;

    return( &_multiply($object,$argument,undef) );
}

sub _assign_exponent {
    my($object,$arg) = @_;
    return ( &_exponent($object,$arg,undef) );
}

sub _equal
{
    my($object,$argument,$flag) = @_;
    my($name) = "'=='";
    my($rows,$cols) = ($object->[1],$object->[2]);
    my($i,$j,$result);

    if ((defined $argument) && ref($argument) &&
        (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
    {
        $result = 1;
        EQUAL:
        for ( $i = 0; $i < $rows; $i++ )
        {
            for ( $j = 0; $j < $cols; $j++ )
            {
                if ($object->[0][$i][$j] != $argument->[0][$i][$j])
                {
                    $result = 0;
                    last EQUAL;
                }
            }
        }
        return($result);
    }
    else
    {
        croak "Math::MatrixReal $name: wrong argument type";
    }
}

sub _not_equal
{
    my($object,$argument,$flag) = @_;
    my($name) = "'!='";
    my($rows,$cols) = ($object->[1],$object->[2]);

    if ((defined $argument) && ref($argument) &&
        (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
    {
	my ($r,$c) = $argument->dim;
	return 1 unless ($r == $rows && $c == $cols );
    my $result = 0;
        NOTEQUAL:
        for ( my $i = 0; $i < $rows; $i++ )
        {
            for ( my $j = 0; $j < $cols; $j++ )
            {
                if ($object->[0][$i][$j] != $argument->[0][$i][$j])
                {
                    $result = 1;
                    last NOTEQUAL;
                }
            }
        }
        return $result;
    } else {
        croak "Math::MatrixReal $name: wrong argument type";
    }
}

sub _less_than
{
    my($object,$argument,$flag) = @_;
    my($name) = "'<'";

    if ((defined $argument) && ref($argument) &&
        (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
    {
        if ((defined $flag) && $flag)
        {
            return( $argument->norm_one() < $object->norm_one() );
        } else {
            return( $object->norm_one() < $argument->norm_one() );
        }
    }
    elsif ((defined $argument) && !(ref($argument)))
    {
        if ((defined $flag) && $flag)
        {
            return( abs($argument) < $object->norm_one() );
        } else {
            return( $object->norm_one() < abs($argument) );
        }
    } else {
        croak "Math::MatrixReal $name: wrong argument type";
    }
}

sub _less_than_or_equal
{
    my($object,$argument,$flag) = @_;
    my($name) = "'<='";

    if ((defined $argument) && ref($argument) &&
        (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
    {
        if ((defined $flag) && $flag)
        {
            return( $argument->norm_one() <= $object->norm_one() );
        } else {
            return( $object->norm_one() <= $argument->norm_one() );
        }
    } elsif ((defined $argument) && !(ref($argument))) {
        if ((defined $flag) && $flag)
        {
            return( abs($argument) <= $object->norm_one() );
        } else {
            return( $object->norm_one() <= abs($argument) );
        }
    } else {
        croak "Math::MatrixReal $name: wrong argument type";
    }
}

sub _greater_than
{
    my($object,$argument,$flag) = @_;
    my($name) = "'>'";

    if ((defined $argument) && ref($argument) &&
        (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
    {
        if ((defined $flag) && $flag)
        {
            return( $argument->norm_one() > $object->norm_one() );
        } else {
            return( $object->norm_one() > $argument->norm_one() );
        }
    } elsif ((defined $argument) && !(ref($argument))) {
        if ((defined $flag) && $flag)
        {
            return( abs($argument) > $object->norm_one() );
        } else {
            return( $object->norm_one() > abs($argument) );
        }
    } else {
        croak "Math::MatrixReal $name: wrong argument type";
    }
}

sub _greater_than_or_equal
{
    my($object,$argument,$flag) = @_;
    my($name) = "'>='";

    if ((defined $argument) && ref($argument) &&
        (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
    {
        if ((defined $flag) && $flag)
        {
            return( $argument->norm_one() >= $object->norm_one() );
        } else {
            return( $object->norm_one() >= $argument->norm_one() );
        }
    } elsif ((defined $argument) && !(ref($argument))) {
        if ((defined $flag) && $flag)
        {
            return( abs($argument) >= $object->norm_one() );
        } else {
            return( $object->norm_one() >= abs($argument) );
        }
    } else {
        croak "Math::MatrixReal $name: wrong argument type";
    }
}

sub _clone
{
    my($object) = @_;

    my $temp = $object->new($object->[1],$object->[2]);
    $temp->copy($object);
    $temp->_undo_LR();
    return $temp;
}
{ no warnings; 42 }
__END__