Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
tree: ebbcd265ff
Fetching contributors…

Cannot retrieve contributors at this time

5390 lines (4199 sloc) 157.952 kb
# Copyright (c) 1996, 1997 by Steffen Beyer. All rights reserved.
# Copyright (c) 1999 by Rodolphe Ortalo. All rights reserved.
# Copyright (c) 2001-2011 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.09';
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 (msouth@fulcrum.org)
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 (msouth@fulcrum.org)
sub new_from_columns {
my $self = shift;
$self->new_from_cols(@_);
}
# from Math::MatrixReal::Ext1 (msouth@fulcrum.org)
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 );
}
# from Math::MatrixReal::Ext1 (msouth@fulcrum.org)
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 _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 = <<LATEX;
\\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;
}
########################################
# #
# 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__
=head1 FUNCTIONS
=head2 Constructor Methods And Such
=over 4
=item * use Math::MatrixReal;
Makes the methods and overloaded operators of this module available
to your program.
=item * $new_matrix = new Math::MatrixReal($rows,$columns);
The matrix object constructor method. A new matrix of size $rows by $columns
will be created, with the value C<0.0> for all elements.
Note that this method is implicitly called by many of the other methods
in this module.
=item * $new_matrix = $some_matrix-E<gt>new($rows,$columns);
Another way of calling the matrix object constructor method.
Matrix $some_matrix is not changed by this in any way.
=item * $new_matrix = $matrix-E<gt>new_from_cols( [ $column_vector|$array_ref|$string, ... ] )
Creates a new matrix given a reference to an array of any of the following:
=over 4
=item * column vectors ( n by 1 Math::MatrixReal matrices )
=item * references to arrays
=item * strings properly formatted to create a column with Math::MatrixReal's
new_from_string command
=back
You may mix and match these as you wish. However, all must be of the
same dimension--no padding happens automatically. Example:
my $matrix = Math::MatrixReal->new_from_cols( [ [1,2], [3,4] ] );
print $matrix;
will print
[ 1.000000000000E+00 3.000000000000E+00 ]
[ 2.000000000000E+00 4.000000000000E+00 ]
=item * new_from_rows( [ $row_vector|$array_ref|$string, ... ] )
Creates a new matrix given a reference to an array of any of the following:
=over 4
=item * row vectors ( 1 by n Math::MatrixReal matrices )
=item * references to arrays
=item * strings properly formatted to create a row with Math::MatrixReal's new_from_string command
=back
You may mix and match these as you wish. However, all must be of the
same dimension--no padding happens automatically. Example:
my $matrix = Math::MatrixReal->new_from_rows( [ [1,2], [3,4] ] );
print $matrix;
will print
[ 1.000000000000E+00 2.000000000000E+00 ]
[ 3.000000000000E+00 4.000000000000E+00 ]
=item * $new_matrix = Math::MatrixReal-E<gt>new_random($rows, $cols, %options );
This method allows you to create a random matrix with various properties controlled
by the %options matrix, which is optional. The default values of the %options matrix
are { integer => 0, symmetric => 0, tridiagonal => 0, diagonal => 0, bounded_by => [0,10] } .
Example:
$matrix = Math::MatrixReal->new_random(4, { diagonal => 1, integer => 1 } );
print $matrix;
will print a 4x4 random diagonal matrix with integer entries between zero and ten, something like
[ 5.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 ]
[ 0.000000000000E+00 2.000000000000E+00 0.000000000000E+00 0.000000000000E+00 ]
[ 0.000000000000E+00 0.000000000000E+00 1.000000000000E+00 0.000000000000E+00 ]
[ 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 8.000000000000E+00 ]
=item * $new_matrix = Math::MatrixReal-E<gt>new_diag( $array_ref );
This method allows you to create a diagonal matrix by only specifying
the diagonal elements. Example:
$matrix = Math::MatrixReal->new_diag( [ 1,2,3,4 ] );
print $matrix;
will print
[ 1.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 ]
[ 0.000000000000E+00 2.000000000000E+00 0.000000000000E+00 0.000000000000E+00 ]
[ 0.000000000000E+00 0.000000000000E+00 3.000000000000E+00 0.000000000000E+00 ]
[ 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 4.000000000000E+00 ]
=item * $new_matrix = Math::MatrixReal-E<gt>new_tridiag( $lower, $diag, $upper );
This method allows you to create a tridiagonal matrix by only specifying
the lower diagonal, diagonal and upper diagonal, respectively.
$matrix = Math::MatrixReal->new_tridiag( [ 6, 4, 2 ], [1,2,3,4], [1, 8, 9] );
print $matrix;
will print
[ 1.000000000000E+00 1.000000000000E+00 0.000000000000E+00 0.000000000000E+00 ]
[ 6.000000000000E+00 2.000000000000E+00 8.000000000000E+00 0.000000000000E+00 ]
[ 0.000000000000E+00 4.000000000000E+00 3.000000000000E+00 9.000000000000E+00 ]
[ 0.000000000000E+00 0.000000000000E+00 2.000000000000E+00 4.000000000000E+00 ]
=item * $new_matrix = Math::MatrixReal-E<gt>new_from_string($string);
This method allows you to read in a matrix from a string (for
instance, from the keyboard, from a file or from your code).
The syntax is simple: each row must start with "C<[ >" and end with
"C< ]\n>" ("C<\n>" being the newline character and "C< >" a space or
tab) and contain one or more numbers, all separated from each other
by spaces or tabs.
Additional spaces or tabs can be added at will, but no comments.
Examples:
$string = "[ 1 2 3 ]\n[ 2 2 -1 ]\n[ 1 1 1 ]\n";
$matrix = Math::MatrixReal->new_from_string($string);
print "$matrix";
By the way, this prints
[ 1.000000000000E+00 2.000000000000E+00 3.000000000000E+00 ]
[ 2.000000000000E+00 2.000000000000E+00 -1.000000000000E+00 ]
[ 1.000000000000E+00 1.000000000000E+00 1.000000000000E+00 ]
But you can also do this in a much more comfortable way using the
shell-like "here-document" syntax:
$matrix = Math::MatrixReal->new_from_string(<<'MATRIX');
[ 1 0 0 0 0 0 1 ]
[ 0 1 0 0 0 0 0 ]
[ 0 0 1 0 0 0 0 ]
[ 0 0 0 1 0 0 0 ]
[ 0 0 0 0 1 0 0 ]
[ 0 0 0 0 0 1 0 ]
[ 1 0 0 0 0 0 -1 ]
MATRIX
You can even use variables in the matrix:
$c1 = 2 / 3;
$c2 = -2 / 5;
$c3 = 26 / 9;
$matrix = Math::MatrixReal->new_from_string(<<"MATRIX");
[ 3 2 0 ]
[ 0 3 2 ]
[ $c1 $c2 $c3 ]
MATRIX
(Remember that you may use spaces and tabs to format the matrix to
your taste)
Note that this method uses exactly the same representation for a
matrix as the "stringify" operator "": this means that you can convert
any matrix into a string with C<$string = "$matrix";> and read it back
in later (for instance from a file!).
Note however that you may suffer a precision loss in this process
because only 13 digits are supported in the mantissa when printed!!
If the string you supply (or someone else supplies) does not obey
the syntax mentioned above, an exception is raised, which can be
caught by "eval" as follows:
print "Please enter your matrix (in one line): ";
$string = <STDIN>;
$string =~ s/\\n/\n/g;
eval { $matrix = Math::MatrixReal->new_from_string($string); };
if ($@)
{
print "$@";
# ...
# (error handling)
}
else
{
# continue...
}
or as follows:
eval { $matrix = Math::MatrixReal->new_from_string(<<"MATRIX"); };
[ 3 2 0 ]
[ 0 3 2 ]
[ $c1 $c2 $c3 ]
MATRIX
if ($@)
# ...
Actually, the method shown above for reading a matrix from the keyboard
is a little awkward, since you have to enter a lot of "\n"'s for the
newlines.
A better way is shown in this piece of code:
while (1)
{
print "\nPlease enter your matrix ";
print "(multiple lines, <ctrl-D> = done):\n";
eval { $new_matrix =
Math::MatrixReal->new_from_string(join('',<STDIN>)); };
if ($@)
{
$@ =~ s/\s+at\b.*?$//;
print "${@}Please try again.\n";
}
else { last; }
}
Possible error messages of the "new_from_string()" method are:
Math::MatrixReal::new_from_string(): syntax error in input string
Math::MatrixReal::new_from_string(): empty input string
If the input string has rows with varying numbers of columns,
the following warning will be printed to STDERR:
Math::MatrixReal::new_from_string(): missing elements will be set to zero!
If everything is okay, the method returns an object reference to the
(newly allocated) matrix containing the elements you specified.
=item * $new_matrix = $some_matrix-E<gt>shadow();
Returns an object reference to a B<NEW> but B<EMPTY> matrix
(filled with zero's) of the B<SAME SIZE> as matrix "C<$some_matrix>".
Matrix "C<$some_matrix>" is not changed by this in any way.
=item * $matrix1-E<gt>copy($matrix2);
Copies the contents of matrix "C<$matrix2>" to an B<ALREADY EXISTING>
matrix "C<$matrix1>" (which must have the same size as matrix "C<$matrix2>"!).
Matrix "C<$matrix2>" is not changed by this in any way.
=item * $twin_matrix = $some_matrix-E<gt>clone();
Returns an object reference to a B<NEW> matrix of the B<SAME SIZE> as
matrix "C<$some_matrix>". The contents of matrix "C<$some_matrix>" have
B<ALREADY BEEN COPIED> to the new matrix "C<$twin_matrix>". This
is the method that the operator "=" is overloaded to when you type
C<$a = $b>, when C<$a> and C<$b> are matrices.
Matrix "C<$some_matrix>" is not changed by this in any way.
=back
=head2 Matrix Row, Column and Element operations
=over 4
=item * $row_vector = $matrix-E<gt>row($row);
This is a projection method which returns an object reference to
a B<NEW> matrix (which in fact is a (row) vector since it has only
one row) to which row number "C<$row>" of matrix "C<$matrix>" has
already been copied.
Matrix "C<$matrix>" is not changed by this in any way.
=item * $column_vector = $matrix-E<gt>column($column);
This is a projection method which returns an object reference to
a B<NEW> matrix (which in fact is a (column) vector since it has
only one column) to which column number "C<$column>" of matrix
"C<$matrix>" has already been copied.
Matrix "C<$matrix>" is not changed by this in any way.
=item * $matrix-E<gt>assign($row,$column,$value);
Explicitly assigns a value "C<$value>" to a single element of the
matrix "C<$matrix>", located in row "C<$row>" and column "C<$column>",
thereby replacing the value previously stored there.
=item * $value = $matrix-E<gt>element($row,$column);
Returns the value of a specific element of the matrix "C<$matrix>",
located in row "C<$row>" and column "C<$column>".
=item * $new_matrix = $matrix-E<gt>each( \&function );
Creates a new matrix by evaluating a code reference on each element of the
given matrix. The function is passed the element, the row index and the column
index, in that order. The value the function returns ( or the value of the last
executed statement ) is the value given to the corresponding element in $new_matrix.
Example:
# add 1 to every element in the matrix
$matrix = $matrix->each ( sub { (shift) + 1 } );
Example:
my $cofactor = $matrix->each( sub { my(undef,$i,$j) = @_;
($i+$j) % 2 == 0 ? $matrix->minor($i,$j)->det()
: -1*$matrix->minor($i,$j)->det();
} );
This code needs some explanation. For each element of $matrix, it throws away the actual value
and stores the row and column indexes in $i and $j. Then it sets element [$i,$j] in $cofactor
to the determinant of C<$matrix-E<gt>minor($i,$j)> if it is an "even" element, or C<-1*$matrix-E<gt>minor($i,$j)>
if it is an "odd" element.
=item * $new_matrix = $matrix-E<gt>each_diag( \&function );
Creates a new matrix by evaluating a code reference on each diagonal element of the
given matrix. The function is passed the element, the row index and the column
index, in that order. The value the function returns ( or the value of the last
executed statement ) is the value given to the corresponding element in $new_matrix.
=item * $matrix-E<gt>swap_col( $col1, $col2 );
This method takes two one-based column numbers and swaps the values of each element in each column.
C<$matrix-E<gt>swap_col(2,3)> would replace column 2 in $matrix with column 3, and replace column
3 with column 2.
=item * $matrix-E<gt>swap_row( $row1, $row2 );
This method takes two one-based row numbers and swaps the values of each element in each row.
C<$matrix-E<gt>swap_row(2,3)> would replace row 2 in $matrix with row 3, and replace row
3 with row 2.
=item * $matrix-E<gt>assign_row( $row_number , $new_row_vector );
This method takes a one-based row number and assigns row $row_number of $matrix
with $new_row_vector and returns the resulting matrix.
C<$matrix-E<gt>assign_row(5, $x)> would replace row 2 in $matrix with the row vector $x.
=back
=head2 Matrix Operations
=over 4
=item * $det = $matrix-E<gt>det();
Returns the determinant of the matrix, without going through
the rigamarole of computing a LR decomposition. This method should
be much faster than LR decomposition if the matrix is diagonal or
triangular. Otherwise, it is just a wrapper for
C<$matrix-E<gt>decompose_LR-E<gt>det_LR>. If the determinant is zero,
there is no inverse and vice-versa. Only quadratic matrices have
determinants.
=item *
C<$inverse = $matrix-E<gt>inverse();>
Returns the inverse of a matrix, without going through the
rigamarole of computing a LR decomposition. If no inverse exists,
undef is returned and an error is printed via C<carp()>.
This is nothing but a wrapper for C<$matrix-E<gt>decompose_LR-E<gt>invert_LR>.
=item *
C<($rows,$columns) = $matrix-E<gt>dim();>
Returns a list of two items, representing the number of rows
and columns the given matrix "C<$matrix>" contains.
=item *
C<$norm_one = $matrix-E<gt>norm_one();>
Returns the "one"-norm of the given matrix "C<$matrix>".
The "one"-norm is defined as follows:
For each column, the sum of the absolute values of the elements in the
different rows of that column is calculated. Finally, the maximum
of these sums is returned.
Note that the "one"-norm and the "maximum"-norm are mathematically
equivalent, although for the same matrix they usually yield a different
value.
Therefore, you should only compare values that have been calculated
using the same norm!
Throughout this package, the "one"-norm is (arbitrarily) used
for all comparisons, for the sake of uniformity and comparability,
except for the iterative methods "solve_GSM()", "solve_SSM()" and
"solve_RM()" which use either norm depending on the matrix itself.
=item *
C<$norm_max = $matrix-E<gt>norm_max();>
Returns the "maximum"-norm of the given matrix $matrix.
The "maximum"-norm is defined as follows: