[system] / trunk / pg / lib / MatrixReal1.pm Repository:
ViewVC logotype

View of /trunk/pg/lib/MatrixReal1.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1050 - (download) (as text) (annotate)
Fri Jun 6 21:39:42 2003 UTC (16 years, 7 months ago) by sh002i
File size: 106667 byte(s)
moved PG modules and macro files from webwork-modperl to pg
-sam

    1 #  Copyright (c) 1996, 1997 by Steffen Beyer. All rights reserved.
    2 #  Copyright (c) 1999 by Rodolphe Ortalo. All rights reserved.
    3 #  This package is free software; you can redistribute it and/or
    4 #  modify it under the same terms as Perl itself.
    5 
    6 # slightly modified for use in WeBWorK
    7 # modifications by Michael E Gage -- added a reference to options in the object array ($this)
    8 # a better approach would be to rewrite this package so that $this is a hash rather than an array
    9 # grep for MEG to see changes.
   10 
   11 # Changed package name to MatrixReal1 throughout.
   12 package MatrixReal1;
   13 
   14 use strict;
   15 use vars qw(@ISA @EXPORT @EXPORT_OK %EXPORT_TAGS $VERSION);
   16 
   17 require Exporter;
   18 
   19 @ISA = qw(Exporter);
   20 
   21 @EXPORT = qw();
   22 
   23 @EXPORT_OK = qw(min max);
   24 
   25 %EXPORT_TAGS = (all => [@EXPORT_OK]);
   26 
   27 $VERSION = '1.3a5';
   28 
   29 use Carp;
   30 
   31 use overload
   32      'neg' => '_negate',
   33        '~' => '_transpose',
   34     'bool' => '_boolean',
   35        '!' => '_not_boolean',
   36       '""' => '_stringify',
   37      'abs' => '_norm',
   38        '+' => '_add',
   39        '-' => '_subtract',
   40        '*' => '_multiply',
   41       '+=' => '_assign_add',
   42       '-=' => '_assign_subtract',
   43       '*=' => '_assign_multiply',
   44       '==' => '_equal',
   45       '!=' => '_not_equal',
   46        '<' => '_less_than',
   47       '<=' => '_less_than_or_equal',
   48        '>' => '_greater_than',
   49       '>=' => '_greater_than_or_equal',
   50       'eq' => '_equal',
   51       'ne' => '_not_equal',
   52       'lt' => '_less_than',
   53       'le' => '_less_than_or_equal',
   54       'gt' => '_greater_than',
   55       'ge' => '_greater_than_or_equal',
   56        '=' => '_clone',
   57 'fallback' =>   undef;
   58 
   59 sub new
   60 {
   61     croak "Usage: \$new_matrix = MatrixReal1->new(\$rows,\$columns);"
   62       if (@_ != 3);
   63 
   64     my $proto = shift;
   65     my $class = ref($proto) || $proto || 'MatrixReal1';
   66     my $rows = shift;
   67     my $cols = shift;
   68     my($i,$j);
   69     my($this);
   70 
   71     croak "MatrixReal1::new(): number of rows must be > 0"
   72       if ($rows <= 0);
   73 
   74     croak "MatrixReal1::new(): number of columns must be > 0"
   75       if ($cols <= 0);
   76 
   77 #    $this = [ [ ], $rows, $cols ];
   78     $this = [ [ ], $rows, $cols,{} ];  # added a holder for options MEG
   79                                        # see also modifications to LR decomposition
   80 
   81     # Creates first empty row
   82     my $empty = [ ];
   83     $#$empty = $cols - 1; # Lengthens the array
   84     for (my $j = 0; $j < $cols; $j++)
   85     {
   86   $empty->[$j] = 0.0;
   87     }
   88     $this->[0][0] = $empty;
   89     # Creates other rows (by copying)
   90     for (my $i = 1; $i < $rows; $i++)
   91     {
   92   my $arow = [ ];
   93   @$arow = @$empty;
   94   $this->[0][$i] = $arow;
   95     }
   96     bless($this, $class);
   97     return($this);
   98 }
   99 
  100 sub new_from_string
  101 {
  102     croak "Usage: \$new_matrix = MatrixReal1->new_from_string(\$string);"
  103       if (@_ != 2);
  104 
  105     my $proto  = shift;
  106     my $class  = ref($proto) || $proto || 'MatrixReal1';
  107     my $string = shift;
  108     my($line,$values);
  109     my($rows,$cols);
  110     my($row,$col);
  111     my($warn);
  112     my($this);
  113 
  114     $warn = 0;
  115     $rows = 0;
  116     $cols = 0;
  117     $values = [ ];
  118     while ($string =~ m!^\s*
  119   \[ \s+ ( (?: [+-]? \d+ (?: \. \d* )? (?: E [+-]? \d+ )? \s+ )+ ) \] \s*? \n
  120     !x)
  121     {
  122         $line = $1;
  123         $string = $';
  124         $values->[$rows] = [ ];
  125         @{$values->[$rows]} = split(' ', $line);
  126         $col = @{$values->[$rows]};
  127         if ($col != $cols)
  128         {
  129             unless ($cols == 0) { $warn = 1; }
  130             if ($col > $cols) { $cols = $col; }
  131         }
  132         $rows++;
  133     }
  134     if ($string !~ m!^\s*$!)
  135     {
  136         croak "MatrixReal1::new_from_string(): syntax error in input string";
  137     }
  138     if ($rows == 0)
  139     {
  140         croak "MatrixReal1::new_from_string(): empty input string";
  141     }
  142     if ($warn)
  143     {
  144         warn "MatrixReal1::new_from_string(): missing elements will be set to zero!\n";
  145     }
  146     $this = MatrixReal1::new($class,$rows,$cols);
  147     for ( $row = 0; $row < $rows; $row++ )
  148     {
  149         for ( $col = 0; $col < @{$values->[$row]}; $col++ )
  150         {
  151             $this->[0][$row][$col] = $values->[$row][$col];
  152         }
  153     }
  154     return($this);
  155 }
  156 
  157 sub shadow
  158 {
  159     croak "Usage: \$new_matrix = \$some_matrix->shadow();"
  160       if (@_ != 1);
  161 
  162     my($matrix) = @_;
  163     my($temp);
  164 
  165     $temp = $matrix->new($matrix->[1],$matrix->[2]);
  166     return($temp);
  167 }
  168 
  169 
  170 sub copy
  171 {
  172     croak "Usage: \$matrix1->copy(\$matrix2);"
  173       if (@_ != 2);
  174 
  175     my($matrix1,$matrix2) = @_;
  176     my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
  177     my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
  178     my($i,$j);
  179 
  180     croak "MatrixReal1::copy(): matrix size mismatch"
  181       unless (($rows1 == $rows2) && ($cols1 == $cols2));
  182 
  183     for ( $i = 0; $i < $rows1; $i++ )
  184     {
  185   my $r1 = []; # New array ref
  186   my $r2 = $matrix2->[0][$i];
  187   @$r1 = @$r2; # Copy whole array directly
  188   $matrix1->[0][$i] = $r1;
  189     }
  190         $matrix1->[3] = $matrix2->[3]; # sign or option
  191     if (defined $matrix2->[4]) # is an LR decomposition matrix!
  192     {
  193     #    $matrix1->[3] = $matrix2->[3]; # $sign
  194         $matrix1->[4] = $matrix2->[4]; # $perm_row
  195         $matrix1->[5] = $matrix2->[5]; # $perm_col
  196         $matrix1->[6] = $matrix2->[6]; # $option
  197     }
  198 }
  199 
  200 sub clone
  201 {
  202     croak "Usage: \$twin_matrix = \$some_matrix->clone();"
  203       if (@_ != 1);
  204 
  205     my($matrix) = @_;
  206     my($temp);
  207 
  208     $temp = $matrix->new($matrix->[1],$matrix->[2]);
  209     $temp->copy($matrix);
  210     return($temp);
  211 }
  212 
  213 sub row
  214 {
  215     croak "Usage: \$row_vector = \$matrix->row(\$row);"
  216       if (@_ != 2);
  217 
  218     my($matrix,$row) = @_;
  219     my($rows,$cols) = ($matrix->[1],$matrix->[2]);
  220     my($temp);
  221     my($j);
  222 
  223     croak "MatrixReal1::row(): row index out of range"
  224       if (($row < 1) || ($row > $rows));
  225 
  226     $row--;
  227     $temp = $matrix->new(1,$cols);
  228     for ( $j = 0; $j < $cols; $j++ )
  229     {
  230         $temp->[0][0][$j] = $matrix->[0][$row][$j];
  231     }
  232     return($temp);
  233 }
  234 
  235 sub column
  236 {
  237     croak "Usage: \$column_vector = \$matrix->column(\$column);"
  238       if (@_ != 2);
  239 
  240     my($matrix,$col) = @_;
  241     my($rows,$cols) = ($matrix->[1],$matrix->[2]);
  242     my($temp);
  243     my($i);
  244 
  245     croak "MatrixReal1::column(): column index out of range"
  246       if (($col < 1) || ($col > $cols));
  247 
  248     $col--;
  249     $temp = $matrix->new($rows,1);
  250     for ( $i = 0; $i < $rows; $i++ )
  251     {
  252         $temp->[0][$i][0] = $matrix->[0][$i][$col];
  253     }
  254     return($temp);
  255 }
  256 
  257 sub _undo_LR
  258 {
  259     croak "Usage: \$matrix->_undo_LR();"
  260       if (@_ != 1);
  261 
  262     my($this) = @_;
  263     my $rh_options = $this->[6];
  264     undef $this->[3];
  265     undef $this->[4];
  266     undef $this->[5];
  267     undef $this->[6];
  268     $this->[3] = $rh_options;
  269 }
  270 
  271 sub zero
  272 {
  273     croak "Usage: \$matrix->zero();"
  274       if (@_ != 1);
  275 
  276     my($this) = @_;
  277     my($rows,$cols) = ($this->[1],$this->[2]);
  278     my($i,$j);
  279 
  280     $this->_undo_LR();
  281 
  282     # Zero first row
  283     for (my $j = 0; $j < $cols; $j++ )
  284     {
  285   $this->[0][0][$j] = 0.0;
  286     }
  287     # Then propagate to other rows
  288     for (my $i = 0; $i < $rows; $i++)
  289     {
  290   @{$this->[0][$i]} = @{$this->[0][0]};
  291     }
  292 }
  293 
  294 sub one
  295 {
  296     croak "Usage: \$matrix->one();"
  297       if (@_ != 1);
  298 
  299     my($this) = @_;
  300     my($rows,$cols) = ($this->[1],$this->[2]);
  301     my($i,$j);
  302 
  303 # No need for this: done by the 'zero()'
  304 #    $this->_undo_LR();
  305     $this->zero(); # We rely on zero() efficiency
  306     for (my $i = 0; $i < $rows; $i++ )
  307     {
  308         $this->[0][$i][$i] = 1.0;
  309     }
  310 }
  311 
  312 sub assign
  313 {
  314     croak "Usage: \$matrix->assign(\$row,\$column,\$value);"
  315       if (@_ != 4);
  316 
  317     my($this,$row,$col,$value) = @_;
  318     my($rows,$cols) = ($this->[1],$this->[2]);
  319 
  320     croak "MatrixReal1::assign(): row index out of range"
  321       if (($row < 1) || ($row > $rows));
  322 
  323     croak "MatrixReal1::assign(): column index out of range"
  324       if (($col < 1) || ($col > $cols));
  325 
  326     $this->_undo_LR();
  327 
  328     $this->[0][--$row][--$col] = $value;
  329 }
  330 
  331 sub element
  332 {
  333     croak "Usage: \$value = \$matrix->element(\$row,\$column);"
  334       if (@_ != 3);
  335 
  336     my($this,$row,$col) = @_;
  337     my($rows,$cols) = ($this->[1],$this->[2]);
  338 
  339     croak "MatrixReal1::element(): row index out of range"
  340       if (($row < 1) || ($row > $rows));
  341 
  342     croak "MatrixReal1::element(): column index out of range"
  343       if (($col < 1) || ($col > $cols));
  344 
  345     return( $this->[0][--$row][--$col] );
  346 }
  347 
  348 sub dim  #  returns dimensions of a matrix
  349 {
  350     croak "Usage: (\$rows,\$columns) = \$matrix->dim();"
  351       if (@_ != 1);
  352 
  353     my($matrix) = @_;
  354 
  355     return( $matrix->[1], $matrix->[2] );
  356 }
  357 
  358 sub norm_one  #  maximum of sums of each column
  359 {
  360     croak "Usage: \$norm_one = \$matrix->norm_one();"
  361       if (@_ != 1);
  362 
  363     my($this) = @_;
  364     my($rows,$cols) = ($this->[1],$this->[2]);
  365 
  366     my $max = 0.0;
  367     for (my $j = 0; $j < $cols; $j++)
  368     {
  369         my $sum = 0.0;
  370         for (my $i = 0; $i < $rows; $i++)
  371         {
  372             $sum += abs( $this->[0][$i][$j] );
  373         }
  374   $max = $sum if ($sum > $max);
  375     }
  376     return($max);
  377 }
  378 
  379 sub norm_max  #  maximum of sums of each row
  380 {
  381     croak "Usage: \$norm_max = \$matrix->norm_max();"
  382       if (@_ != 1);
  383 
  384     my($this) = @_;
  385     my($rows,$cols) = ($this->[1],$this->[2]);
  386 
  387     my $max = 0.0;
  388     for (my $i = 0; $i < $rows; $i++)
  389     {
  390         my $sum = 0.0;
  391         for (my $j = 0; $j < $cols; $j++)
  392         {
  393             $sum += abs( $this->[0][$i][$j] );
  394         }
  395   $max = $sum if ($sum > $max);
  396     }
  397     return($max);
  398 }
  399 
  400 sub negate
  401 {
  402     croak "Usage: \$matrix1->negate(\$matrix2);"
  403       if (@_ != 2);
  404 
  405     my($matrix1,$matrix2) = @_;
  406     my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
  407     my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
  408 
  409     croak "MatrixReal1::negate(): matrix size mismatch"
  410       unless (($rows1 == $rows2) && ($cols1 == $cols2));
  411 
  412     $matrix1->_undo_LR();
  413 
  414     for (my $i = 0; $i < $rows1; $i++ )
  415     {
  416         for (my $j = 0; $j < $cols1; $j++ )
  417         {
  418             $matrix1->[0][$i][$j] = -($matrix2->[0][$i][$j]);
  419         }
  420     }
  421 }
  422 
  423 sub transpose
  424 {
  425     croak "Usage: \$matrix1->transpose(\$matrix2);"
  426       if (@_ != 2);
  427 
  428     my($matrix1,$matrix2) = @_;
  429     my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
  430     my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
  431 
  432     croak "MatrixReal1::transpose(): matrix size mismatch"
  433       unless (($rows1 == $cols2) && ($cols1 == $rows2));
  434 
  435     $matrix1->_undo_LR();
  436 
  437     if ($rows1 == $cols1)
  438     {
  439         # more complicated to make in-place possible!
  440 
  441         for (my $i = 0; $i < $rows1; $i++)
  442         {
  443             for (my $j = ($i + 1); $j < $cols1; $j++)
  444             {
  445                 my $swap              = $matrix2->[0][$i][$j];
  446                 $matrix1->[0][$i][$j] = $matrix2->[0][$j][$i];
  447                 $matrix1->[0][$j][$i] = $swap;
  448             }
  449             $matrix1->[0][$i][$i] = $matrix2->[0][$i][$i];
  450         }
  451     }
  452     else # ($rows1 != $cols1)
  453     {
  454         for (my $i = 0; $i < $rows1; $i++)
  455         {
  456             for (my $j = 0; $j < $cols1; $j++)
  457             {
  458                 $matrix1->[0][$i][$j] = $matrix2->[0][$j][$i];
  459             }
  460         }
  461     }
  462 }
  463 
  464 sub add
  465 {
  466     croak "Usage: \$matrix1->add(\$matrix2,\$matrix3);"
  467       if (@_ != 3);
  468 
  469     my($matrix1,$matrix2,$matrix3) = @_;
  470     my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
  471     my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
  472     my($rows3,$cols3) = ($matrix3->[1],$matrix3->[2]);
  473     my($i,$j);
  474 
  475     croak "MatrixReal1::add(): matrix size mismatch"
  476       unless (($rows1 == $rows2) && ($rows1 == $rows3) &&
  477               ($cols1 == $cols2) && ($cols1 == $cols3));
  478 
  479     $matrix1->_undo_LR();
  480 
  481     for ( $i = 0; $i < $rows1; $i++ )
  482     {
  483         for ( $j = 0; $j < $cols1; $j++ )
  484         {
  485             $matrix1->[0][$i][$j] =
  486             $matrix2->[0][$i][$j] + $matrix3->[0][$i][$j];
  487         }
  488     }
  489 }
  490 
  491 sub subtract
  492 {
  493     croak "Usage: \$matrix1->subtract(\$matrix2,\$matrix3);"
  494       if (@_ != 3);
  495 
  496     my($matrix1,$matrix2,$matrix3) = @_;
  497     my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
  498     my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
  499     my($rows3,$cols3) = ($matrix3->[1],$matrix3->[2]);
  500     my($i,$j);
  501 
  502     croak "MatrixReal1::subtract(): matrix size mismatch"
  503       unless (($rows1 == $rows2) && ($rows1 == $rows3) &&
  504               ($cols1 == $cols2) && ($cols1 == $cols3));
  505 
  506     $matrix1->_undo_LR();
  507 
  508     for ( $i = 0; $i < $rows1; $i++ )
  509     {
  510         for ( $j = 0; $j < $cols1; $j++ )
  511         {
  512             $matrix1->[0][$i][$j] =
  513             $matrix2->[0][$i][$j] - $matrix3->[0][$i][$j];
  514         }
  515     }
  516 }
  517 
  518 sub multiply_scalar
  519 {
  520     croak "Usage: \$matrix1->multiply_scalar(\$matrix2,\$scalar);"
  521       if (@_ != 3);
  522 
  523     my($matrix1,$matrix2,$scalar) = @_;
  524     my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
  525     my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
  526     my($i,$j);
  527 
  528     croak "MatrixReal1::multiply_scalar(): matrix size mismatch"
  529       unless (($rows1 == $rows2) && ($cols1 == $cols2));
  530 
  531     $matrix1->_undo_LR();
  532 
  533     for ( $i = 0; $i < $rows1; $i++ )
  534     {
  535         for ( $j = 0; $j < $cols1; $j++ )
  536         {
  537             $matrix1->[0][$i][$j] = $matrix2->[0][$i][$j] * $scalar;
  538         }
  539     }
  540 }
  541 
  542 sub multiply
  543 {
  544     croak "Usage: \$product_matrix = \$matrix1->multiply(\$matrix2);"
  545       if (@_ != 2);
  546 
  547     my($matrix1,$matrix2) = @_;
  548     my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
  549     my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
  550     my($temp);
  551 
  552     croak "MatrixReal1::multiply(): matrix size mismatch"
  553       unless ($cols1 == $rows2);
  554 
  555     $temp = $matrix1->new($rows1,$cols2);
  556     for (my $i = 0; $i < $rows1; $i++ )
  557     {
  558         for (my $j = 0; $j < $cols2; $j++ )
  559         {
  560             my $sum = 0.0;
  561             for (my $k = 0; $k < $cols1; $k++ )
  562             {
  563         $sum += ( $matrix1->[0][$i][$k] * $matrix2->[0][$k][$j] );
  564             }
  565             $temp->[0][$i][$j] = $sum;
  566         }
  567     }
  568     return($temp);
  569 }
  570 
  571 sub min
  572 {
  573     croak "Usage: \$minimum = MatrixReal1::min(\$number1,\$number2);"
  574       if (@_ != 2);
  575 
  576     return( $_[0] < $_[1] ? $_[0] : $_[1] );
  577 }
  578 
  579 sub max
  580 {
  581     croak "Usage: \$maximum = MatrixReal1::max(\$number1,\$number2);"
  582       if (@_ != 2);
  583 
  584     return( $_[0] > $_[1] ? $_[0] : $_[1] );
  585 }
  586 
  587 sub kleene
  588 {
  589     croak "Usage: \$minimal_cost_matrix = \$cost_matrix->kleene();"
  590       if (@_ != 1);
  591 
  592     my($matrix) = @_;
  593     my($rows,$cols) = ($matrix->[1],$matrix->[2]);
  594     my($i,$j,$k,$n);
  595     my($temp);
  596 
  597     croak "MatrixReal1::kleene(): matrix is not quadratic"
  598       unless ($rows == $cols);
  599 
  600     $temp = $matrix->new($rows,$cols);
  601     $temp->copy($matrix);
  602     $temp->_undo_LR();
  603     $n = $rows;
  604     for ( $i = 0; $i < $n; $i++ )
  605     {
  606         $temp->[0][$i][$i] = min( $temp->[0][$i][$i] , 0 );
  607     }
  608     for ( $k = 0; $k < $n; $k++ )
  609     {
  610         for ( $i = 0; $i < $n; $i++ )
  611         {
  612             for ( $j = 0; $j < $n; $j++ )
  613             {
  614                 $temp->[0][$i][$j] = min( $temp->[0][$i][$j] ,
  615                                         ( $temp->[0][$i][$k] +
  616                                           $temp->[0][$k][$j] ) );
  617             }
  618         }
  619     }
  620     return($temp);
  621 }
  622 
  623 sub normalize
  624 {
  625     croak "Usage: (\$norm_matrix,\$norm_vector) = \$matrix->normalize(\$vector);"
  626       if (@_ != 2);
  627 
  628     my($matrix,$vector) = @_;
  629     my($rows,$cols) = ($matrix->[1],$matrix->[2]);
  630     my($norm_matrix,$norm_vector);
  631     my($max,$val);
  632     my($i,$j,$n);
  633 
  634     croak "MatrixReal1::normalize(): matrix is not quadratic"
  635       unless ($rows == $cols);
  636 
  637     $n = $rows;
  638 
  639     croak "MatrixReal1::normalize(): vector is not a column vector"
  640       unless ($vector->[2] == 1);
  641 
  642     croak "MatrixReal1::normalize(): matrix and vector size mismatch"
  643       unless ($vector->[1] == $n);
  644 
  645     $norm_matrix = $matrix->new($n,$n);
  646     $norm_vector = $vector->new($n,1);
  647 
  648     $norm_matrix->copy($matrix);
  649     $norm_vector->copy($vector);
  650 
  651     $norm_matrix->_undo_LR();
  652 
  653     for ( $i = 0; $i < $n; $i++ )
  654     {
  655         $max = abs($norm_vector->[0][$i][0]);
  656         for ( $j = 0; $j < $n; $j++ )
  657         {
  658             $val = abs($norm_matrix->[0][$i][$j]);
  659             if ($val > $max) { $max = $val; }
  660         }
  661         if ($max != 0)
  662         {
  663             $norm_vector->[0][$i][0] /= $max;
  664             for ( $j = 0; $j < $n; $j++ )
  665             {
  666                 $norm_matrix->[0][$i][$j] /= $max;
  667             }
  668         }
  669     }
  670     return($norm_matrix,$norm_vector);
  671 }
  672 
  673 sub decompose_LR
  674 {
  675     croak "Usage: \$LR_matrix = \$matrix->decompose_LR();"
  676       if (@_ != 1);
  677 
  678     my($matrix) = @_;
  679     my($rows,$cols) = ($matrix->[1],$matrix->[2]);
  680     my($perm_row,$perm_col);
  681     my($row,$col,$max);
  682     my($i,$j,$k,$n);
  683     my($sign) = 1;
  684     my($swap);
  685     my($temp);
  686 
  687     croak "MatrixReal1::decompose_LR(): matrix is not quadratic"
  688       unless ($rows == $cols);
  689 
  690     $temp = $matrix->new($rows,$cols);
  691     $temp->copy($matrix);
  692     $n = $rows;
  693     $perm_row = [ ];
  694     $perm_col = [ ];
  695     for ( $i = 0; $i < $n; $i++ )
  696     {
  697         $perm_row->[$i] = $i;
  698         $perm_col->[$i] = $i;
  699     }
  700     NONZERO:
  701     for ( $k = 0; $k < $n; $k++ ) # use Gauss's algorithm:
  702     {
  703         # complete pivot-search:
  704 
  705         $max = 0;
  706         for ( $i = $k; $i < $n; $i++ )
  707         {
  708             for ( $j = $k; $j < $n; $j++ )
  709             {
  710                 if (($swap = abs($temp->[0][$i][$j])) > $max)
  711                 {
  712                     $max = $swap;
  713                     $row = $i;
  714                     $col = $j;
  715                 }
  716             }
  717         }
  718         last NONZERO if ($max == 0); # (all remaining elements are zero)
  719         if ($k != $row) # swap row $k and row $row:
  720         {
  721             $sign = -$sign;
  722             $swap             = $perm_row->[$k];
  723             $perm_row->[$k]   = $perm_row->[$row];
  724             $perm_row->[$row] = $swap;
  725             for ( $j = 0; $j < $n; $j++ )
  726             {
  727                 # (must run from 0 since L has to be swapped too!)
  728 
  729                 $swap                = $temp->[0][$k][$j];
  730                 $temp->[0][$k][$j]   = $temp->[0][$row][$j];
  731                 $temp->[0][$row][$j] = $swap;
  732             }
  733         }
  734         if ($k != $col) # swap column $k and column $col:
  735         {
  736             $sign = -$sign;
  737             $swap             = $perm_col->[$k];
  738             $perm_col->[$k]   = $perm_col->[$col];
  739             $perm_col->[$col] = $swap;
  740             for ( $i = 0; $i < $n; $i++ )
  741             {
  742                 $swap                = $temp->[0][$i][$k];
  743                 $temp->[0][$i][$k]   = $temp->[0][$i][$col];
  744                 $temp->[0][$i][$col] = $swap;
  745             }
  746         }
  747         for ( $i = ($k + 1); $i < $n; $i++ )
  748         {
  749             # scan the remaining rows, add multiples of row $k to row $i:
  750 
  751             $swap = $temp->[0][$i][$k] / $temp->[0][$k][$k];
  752             if ($swap != 0)
  753             {
  754                 # calculate a row of matrix R:
  755 
  756                 for ( $j = ($k + 1); $j < $n; $j++ )
  757                 {
  758                     $temp->[0][$i][$j] -= $temp->[0][$k][$j] * $swap;
  759                 }
  760 
  761                 # store matrix L in same matrix as R:
  762 
  763                 $temp->[0][$i][$k] = $swap;
  764             }
  765         }
  766     }
  767     my $rh_options = $temp->[3];
  768     $temp->[3] = $sign;
  769     $temp->[4] = $perm_row;
  770     $temp->[5] = $perm_col;
  771     $temp->[6] = $temp->[3];
  772     return($temp);
  773 }
  774 
  775 sub solve_LR
  776 {
  777     croak "Usage: (\$dimension,\$x_vector,\$base_matrix) = \$LR_matrix->solve_LR(\$b_vector);"
  778       if (@_ != 2);
  779 
  780     my($LR_matrix,$b_vector) = @_;
  781     my($rows,$cols) = ($LR_matrix->[1],$LR_matrix->[2]);
  782     my($dimension,$x_vector,$base_matrix);
  783     my($perm_row,$perm_col);
  784     my($y_vector,$sum);
  785     my($i,$j,$k,$n);
  786 
  787     croak "MatrixReal1::solve_LR(): not an LR decomposition matrix"
  788       unless ((defined $LR_matrix->[4]) && ($rows == $cols));
  789 
  790     $n = $rows;
  791 
  792     croak "MatrixReal1::solve_LR(): vector is not a column vector"
  793       unless ($b_vector->[2] == 1);
  794 
  795     croak "MatrixReal1::solve_LR(): matrix and vector size mismatch"
  796       unless ($b_vector->[1] == $n);
  797 
  798     $perm_row = $LR_matrix->[4];
  799     $perm_col = $LR_matrix->[5];
  800 
  801     $x_vector    =   $b_vector->new($n,1);
  802     $y_vector    =   $b_vector->new($n,1);
  803     $base_matrix = $LR_matrix->new($n,$n);
  804 
  805     # calculate "x" so that LRx = b  ==>  calculate Ly = b, Rx = y:
  806 
  807     for ( $i = 0; $i < $n; $i++ ) # calculate $y_vector:
  808     {
  809         $sum = $b_vector->[0][($perm_row->[$i])][0];
  810         for ( $j = 0; $j < $i; $j++ )
  811         {
  812             $sum -= $LR_matrix->[0][$i][$j] * $y_vector->[0][$j][0];
  813         }
  814         $y_vector->[0][$i][0] = $sum;
  815     }
  816 
  817     $dimension = 0;
  818     for ( $i = ($n - 1); $i >= 0; $i-- ) # calculate $x_vector:
  819     {
  820         if ($LR_matrix->[0][$i][$i] == 0)
  821         {
  822             if ($y_vector->[0][$i][0] != 0)
  823             {
  824                 return(); # a solution does not exist!
  825             }
  826             else
  827             {
  828                 $dimension++;
  829                 $x_vector->[0][($perm_col->[$i])][0] = 0;
  830             }
  831         }
  832         else
  833         {
  834             $sum = $y_vector->[0][$i][0];
  835             for ( $j = ($i + 1); $j < $n; $j++ )
  836             {
  837                 $sum -= $LR_matrix->[0][$i][$j] *
  838                     $x_vector->[0][($perm_col->[$j])][0];
  839             }
  840             $x_vector->[0][($perm_col->[$i])][0] =
  841                 $sum / $LR_matrix->[0][$i][$i];
  842         }
  843     }
  844     if ($dimension)
  845     {
  846         if ($dimension == $n)
  847         {
  848             $base_matrix->one();
  849         }
  850         else
  851         {
  852             for ( $k = 0; $k < $dimension; $k++ )
  853             {
  854                 $base_matrix->[0][($perm_col->[($n-$k-1)])][$k] = 1;
  855                 for ( $i = ($n-$dimension-1); $i >= 0; $i-- )
  856                 {
  857                     $sum = 0;
  858                     for ( $j = ($i + 1); $j < $n; $j++ )
  859                     {
  860                         $sum -= $LR_matrix->[0][$i][$j] *
  861                             $base_matrix->[0][($perm_col->[$j])][$k];
  862                     }
  863                     $base_matrix->[0][($perm_col->[$i])][$k] =
  864                         $sum / $LR_matrix->[0][$i][$i];
  865                 }
  866             }
  867         }
  868     }
  869     return( $dimension, $x_vector, $base_matrix );
  870 }
  871 
  872 sub invert_LR
  873 {
  874     croak "Usage: \$inverse_matrix = \$LR_matrix->invert_LR();"
  875       if (@_ != 1);
  876 
  877     my($matrix) = @_;
  878     my($rows,$cols) = ($matrix->[1],$matrix->[2]);
  879     my($inv_matrix,$x_vector,$y_vector);
  880     my($i,$j,$n);
  881 
  882     croak "MatrixReal1::invert_LR(): not an LR decomposition matrix"
  883       unless ((defined $matrix->[4]) && ($rows == $cols));
  884 
  885     $n = $rows;
  886     if ($matrix->[0][$n-1][$n-1] != 0)
  887     {
  888         $inv_matrix = $matrix->new($n,$n);
  889         $y_vector   = $matrix->new($n,1);
  890         for ( $j = 0; $j < $n; $j++ )
  891         {
  892             if ($j > 0)
  893             {
  894                 $y_vector->[0][$j-1][0] = 0;
  895             }
  896             $y_vector->[0][$j][0] = 1;
  897             if (($rows,$x_vector,$cols) = $matrix->solve_LR($y_vector))
  898             {
  899                 for ( $i = 0; $i < $n; $i++ )
  900                 {
  901                     $inv_matrix->[0][$i][$j] = $x_vector->[0][$i][0];
  902                 }
  903             }
  904             else
  905             {
  906                 die "MatrixReal1::invert_LR(): unexpected error - please inform author!\n";
  907             }
  908         }
  909         return($inv_matrix);
  910     }
  911     else { return(); } # matrix is not invertible!
  912 }
  913 
  914 sub condition
  915 {
  916     # 1st matrix MUST be the inverse of 2nd matrix (or vice-versa)
  917     # for a meaningful result!
  918 
  919     croak "Usage: \$condition = \$matrix->condition(\$inverse_matrix);"
  920       if (@_ != 2);
  921 
  922     my($matrix1,$matrix2) = @_;
  923     my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
  924     my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
  925 
  926     croak "MatrixReal1::condition(): 1st matrix is not quadratic"
  927       unless ($rows1 == $cols1);
  928 
  929     croak "MatrixReal1::condition(): 2nd matrix is not quadratic"
  930       unless ($rows2 == $cols2);
  931 
  932     croak "MatrixReal1::condition(): matrix size mismatch"
  933       unless (($rows1 == $rows2) && ($cols1 == $cols2));
  934 
  935     return( $matrix1->norm_one() * $matrix2->norm_one() );
  936 }
  937 
  938 sub det_LR  #  determinant of LR decomposition matrix
  939 {
  940     croak "Usage: \$determinant = \$LR_matrix->det_LR();"
  941       if (@_ != 1);
  942 
  943     my($matrix) = @_;
  944     my($rows,$cols) = ($matrix->[1],$matrix->[2]);
  945     my($k,$det);
  946 
  947     croak "MatrixReal1::det_LR(): not an LR decomposition matrix"
  948  #   unless ((defined $matrix->[3]) && ($rows == $cols));
  949     unless ((defined $matrix->[4]) && ($rows == $cols)); #options might be in [3] position-- MEG
  950 
  951     $det = $matrix->[3];   # grab the sign from permutation shifts
  952     for ( $k = 0; $k < $rows; $k++ )
  953     {
  954         $det *= $matrix->[0][$k][$k];
  955     }
  956     return($det);
  957 }
  958 
  959 sub order_LR  #  order of LR decomposition matrix (number of non-zero equations)
  960 {
  961     croak "Usage: \$order = \$LR_matrix->order_LR();"
  962       if (@_ != 1);
  963 
  964     my($matrix) = @_;
  965     my($rows,$cols) = ($matrix->[1],$matrix->[2]);
  966     my($order);
  967 
  968     croak "MatrixReal1::order_LR(): not an LR decomposition matrix"
  969       unless ((defined $matrix->[4]) && ($rows == $cols));
  970 
  971     ZERO:
  972     for ( $order = ($rows - 1); $order >= 0; $order-- )
  973     {
  974         last ZERO if ($matrix->[0][$order][$order] != 0);
  975     }
  976     return(++$order);
  977 }
  978 
  979 sub scalar_product
  980 {
  981     croak "Usage: \$scalar_product = \$vector1->scalar_product(\$vector2);"
  982       if (@_ != 2);
  983 
  984     my($vector1,$vector2) = @_;
  985     my($rows1,$cols1) = ($vector1->[1],$vector1->[2]);
  986     my($rows2,$cols2) = ($vector2->[1],$vector2->[2]);
  987     my($k,$sum);
  988 
  989     croak "MatrixReal1::scalar_product(): 1st vector is not a column vector"
  990       unless ($cols1 == 1);
  991 
  992     croak "MatrixReal1::scalar_product(): 2nd vector is not a column vector"
  993       unless ($cols2 == 1);
  994 
  995     croak "MatrixReal1::scalar_product(): vector size mismatch"
  996       unless ($rows1 == $rows2);
  997 
  998     $sum = 0;
  999     for ( $k = 0; $k < $rows1; $k++ )
 1000     {
 1001         $sum += $vector1->[0][$k][0] * $vector2->[0][$k][0];
 1002     }
 1003     return($sum);
 1004 }
 1005 
 1006 sub vector_product
 1007 {
 1008     croak "Usage: \$vector_product = \$vector1->vector_product(\$vector2);"
 1009       if (@_ != 2);
 1010 
 1011     my($vector1,$vector2) = @_;
 1012     my($rows1,$cols1) = ($vector1->[1],$vector1->[2]);
 1013     my($rows2,$cols2) = ($vector2->[1],$vector2->[2]);
 1014     my($temp);
 1015     my($n);
 1016 
 1017     croak "MatrixReal1::vector_product(): 1st vector is not a column vector"
 1018       unless ($cols1 == 1);
 1019 
 1020     croak "MatrixReal1::vector_product(): 2nd vector is not a column vector"
 1021       unless ($cols2 == 1);
 1022 
 1023     croak "MatrixReal1::vector_product(): vector size mismatch"
 1024       unless ($rows1 == $rows2);
 1025 
 1026     $n = $rows1;
 1027 
 1028     croak "MatrixReal1::vector_product(): only defined for 3 dimensions"
 1029       unless ($n == 3);
 1030 
 1031     $temp = $vector1->new($n,1);
 1032     $temp->[0][0][0] = $vector1->[0][1][0] * $vector2->[0][2][0] -
 1033                        $vector1->[0][2][0] * $vector2->[0][1][0];
 1034     $temp->[0][1][0] = $vector1->[0][2][0] * $vector2->[0][0][0] -
 1035                        $vector1->[0][0][0] * $vector2->[0][2][0];
 1036     $temp->[0][2][0] = $vector1->[0][0][0] * $vector2->[0][1][0] -
 1037                        $vector1->[0][1][0] * $vector2->[0][0][0];
 1038     return($temp);
 1039 }
 1040 
 1041 sub length
 1042 {
 1043     croak "Usage: \$length = \$vector->length();"
 1044       if (@_ != 1);
 1045 
 1046     my($vector) = @_;
 1047     my($rows,$cols) = ($vector->[1],$vector->[2]);
 1048     my($k,$comp,$sum);
 1049 
 1050     croak "MatrixReal1::length(): vector is not a column vector"
 1051       unless ($cols == 1);
 1052 
 1053     $sum = 0;
 1054     for ( $k = 0; $k < $rows; $k++ )
 1055     {
 1056         $comp = $vector->[0][$k][0];
 1057         $sum += $comp * $comp;
 1058     }
 1059     return( sqrt( $sum ) );
 1060 }
 1061 
 1062 sub _init_iteration
 1063 {
 1064     croak "Usage: \$which_norm = \$matrix->_init_iteration();"
 1065       if (@_ != 1);
 1066 
 1067     my($matrix) = @_;
 1068     my($rows,$cols) = ($matrix->[1],$matrix->[2]);
 1069     my($ok,$max,$sum,$norm);
 1070     my($i,$j,$n);
 1071 
 1072     croak "MatrixReal1::_init_iteration(): matrix is not quadratic"
 1073       unless ($rows == $cols);
 1074 
 1075     $ok = 1;
 1076     $n = $rows;
 1077     for ( $i = 0; $i < $n; $i++ )
 1078     {
 1079         if ($matrix->[0][$i][$i] == 0) { $ok = 0; }
 1080     }
 1081     if ($ok)
 1082     {
 1083         $norm = 1; # norm_one
 1084         $max = 0;
 1085         for ( $j = 0; $j < $n; $j++ )
 1086         {
 1087             $sum = 0;
 1088             for ( $i = 0; $i < $j; $i++ )
 1089             {
 1090                 $sum += abs($matrix->[0][$i][$j]);
 1091             }
 1092             for ( $i = ($j + 1); $i < $n; $i++ )
 1093             {
 1094                 $sum += abs($matrix->[0][$i][$j]);
 1095             }
 1096             $sum /= abs($matrix->[0][$j][$j]);
 1097             if ($sum > $max) { $max = $sum; }
 1098         }
 1099         $ok = ($max < 1);
 1100         unless ($ok)
 1101         {
 1102             $norm = -1; # norm_max
 1103             $max = 0;
 1104             for ( $i = 0; $i < $n; $i++ )
 1105             {
 1106                 $sum = 0;
 1107                 for ( $j = 0; $j < $i; $j++ )
 1108                 {
 1109                     $sum += abs($matrix->[0][$i][$j]);
 1110                 }
 1111                 for ( $j = ($i + 1); $j < $n; $j++ )
 1112                 {
 1113                     $sum += abs($matrix->[0][$i][$j]);
 1114                 }
 1115                 $sum /= abs($matrix->[0][$i][$i]);
 1116                 if ($sum > $max) { $max = $sum; }
 1117             }
 1118             $ok = ($max < 1)
 1119         }
 1120     }
 1121     if ($ok) { return($norm); }
 1122     else     { return(0); }
 1123 }
 1124 
 1125 sub solve_GSM  #  Global Step Method
 1126 {
 1127     croak "Usage: \$xn_vector = \$matrix->solve_GSM(\$x0_vector,\$b_vector,\$epsilon);"
 1128       if (@_ != 4);
 1129 
 1130     my($matrix,$x0_vector,$b_vector,$epsilon) = @_;
 1131     my($rows1,$cols1) = (   $matrix->[1],   $matrix->[2]);
 1132     my($rows2,$cols2) = ($x0_vector->[1],$x0_vector->[2]);
 1133     my($rows3,$cols3) = ( $b_vector->[1], $b_vector->[2]);
 1134     my($norm,$sum,$diff);
 1135     my($xn_vector);
 1136     my($i,$j,$n);
 1137 
 1138     croak "MatrixReal1::solve_GSM(): matrix is not quadratic"
 1139       unless ($rows1 == $cols1);
 1140 
 1141     $n = $rows1;
 1142 
 1143     croak "MatrixReal1::solve_GSM(): 1st vector is not a column vector"
 1144       unless ($cols2 == 1);
 1145 
 1146     croak "MatrixReal1::solve_GSM(): 2nd vector is not a column vector"
 1147       unless ($cols3 == 1);
 1148 
 1149     croak "MatrixReal1::solve_GSM(): matrix and vector size mismatch"
 1150       unless (($rows2 == $n) && ($rows3 == $n));
 1151 
 1152     return() unless ($norm = $matrix->_init_iteration());
 1153 
 1154     $xn_vector = $x0_vector->new($n,1);
 1155 
 1156     $diff = $epsilon + 1;
 1157     while ($diff >= $epsilon)
 1158     {
 1159         for ( $i = 0; $i < $n; $i++ )
 1160         {
 1161             $sum = $b_vector->[0][$i][0];
 1162             for ( $j = 0; $j < $i; $j++ )
 1163             {
 1164                 $sum -= $matrix->[0][$i][$j] * $x0_vector->[0][$j][0];
 1165             }
 1166             for ( $j = ($i + 1); $j < $n; $j++ )
 1167             {
 1168                 $sum -= $matrix->[0][$i][$j] * $x0_vector->[0][$j][0];
 1169             }
 1170             $xn_vector->[0][$i][0] = $sum / $matrix->[0][$i][$i];
 1171         }
 1172         $x0_vector->subtract($x0_vector,$xn_vector);
 1173         if ($norm > 0) { $diff = $x0_vector->norm_one(); }
 1174         else           { $diff = $x0_vector->norm_max(); }
 1175         for ( $i = 0; $i < $n; $i++ )
 1176         {
 1177             $x0_vector->[0][$i][0] = $xn_vector->[0][$i][0];
 1178         }
 1179     }
 1180     return($xn_vector);
 1181 }
 1182 
 1183 sub solve_SSM  #  Single Step Method
 1184 {
 1185     croak "Usage: \$xn_vector = \$matrix->solve_SSM(\$x0_vector,\$b_vector,\$epsilon);"
 1186       if (@_ != 4);
 1187 
 1188     my($matrix,$x0_vector,$b_vector,$epsilon) = @_;
 1189     my($rows1,$cols1) = (   $matrix->[1],   $matrix->[2]);
 1190     my($rows2,$cols2) = ($x0_vector->[1],$x0_vector->[2]);
 1191     my($rows3,$cols3) = ( $b_vector->[1], $b_vector->[2]);
 1192     my($norm,$sum,$diff);
 1193     my($xn_vector);
 1194     my($i,$j,$n);
 1195 
 1196     croak "MatrixReal1::solve_SSM(): matrix is not quadratic"
 1197       unless ($rows1 == $cols1);
 1198 
 1199     $n = $rows1;
 1200 
 1201     croak "MatrixReal1::solve_SSM(): 1st vector is not a column vector"
 1202       unless ($cols2 == 1);
 1203 
 1204     croak "MatrixReal1::solve_SSM(): 2nd vector is not a column vector"
 1205       unless ($cols3 == 1);
 1206 
 1207     croak "MatrixReal1::solve_SSM(): matrix and vector size mismatch"
 1208       unless (($rows2 == $n) && ($rows3 == $n));
 1209 
 1210     return() unless ($norm = $matrix->_init_iteration());
 1211 
 1212     $xn_vector = $x0_vector->new($n,1);
 1213     $xn_vector->copy($x0_vector);
 1214 
 1215     $diff = $epsilon + 1;
 1216     while ($diff >= $epsilon)
 1217     {
 1218         for ( $i = 0; $i < $n; $i++ )
 1219         {
 1220             $sum = $b_vector->[0][$i][0];
 1221             for ( $j = 0; $j < $i; $j++ )
 1222             {
 1223                 $sum -= $matrix->[0][$i][$j] * $xn_vector->[0][$j][0];
 1224             }
 1225             for ( $j = ($i + 1); $j < $n; $j++ )
 1226             {
 1227                 $sum -= $matrix->[0][$i][$j] * $xn_vector->[0][$j][0];
 1228             }
 1229             $xn_vector->[0][$i][0] = $sum / $matrix->[0][$i][$i];
 1230         }
 1231         $x0_vector->subtract($x0_vector,$xn_vector);
 1232         if ($norm > 0) { $diff = $x0_vector->norm_one(); }
 1233         else           { $diff = $x0_vector->norm_max(); }
 1234         for ( $i = 0; $i < $n; $i++ )
 1235         {
 1236             $x0_vector->[0][$i][0] = $xn_vector->[0][$i][0];
 1237         }
 1238     }
 1239     return($xn_vector);
 1240 }
 1241 
 1242 sub solve_RM  #  Relaxation Method
 1243 {
 1244     croak "Usage: \$xn_vector = \$matrix->solve_RM(\$x0_vector,\$b_vector,\$weight,\$epsilon);"
 1245       if (@_ != 5);
 1246 
 1247     my($matrix,$x0_vector,$b_vector,$weight,$epsilon) = @_;
 1248     my($rows1,$cols1) = (   $matrix->[1],   $matrix->[2]);
 1249     my($rows2,$cols2) = ($x0_vector->[1],$x0_vector->[2]);
 1250     my($rows3,$cols3) = ( $b_vector->[1], $b_vector->[2]);
 1251     my($norm,$sum,$diff);
 1252     my($xn_vector);
 1253     my($i,$j,$n);
 1254 
 1255     croak "MatrixReal1::solve_RM(): matrix is not quadratic"
 1256       unless ($rows1 == $cols1);
 1257 
 1258     $n = $rows1;
 1259 
 1260     croak "MatrixReal1::solve_RM(): 1st vector is not a column vector"
 1261       unless ($cols2 == 1);
 1262 
 1263     croak "MatrixReal1::solve_RM(): 2nd vector is not a column vector"
 1264       unless ($cols3 == 1);
 1265 
 1266     croak "MatrixReal1::solve_RM(): matrix and vector size mismatch"
 1267       unless (($rows2 == $n) && ($rows3 == $n));
 1268 
 1269     return() unless ($norm = $matrix->_init_iteration());
 1270 
 1271     $xn_vector = $x0_vector->new($n,1);
 1272     $xn_vector->copy($x0_vector);
 1273 
 1274     $diff = $epsilon + 1;
 1275     while ($diff >= $epsilon)
 1276     {
 1277         for ( $i = 0; $i < $n; $i++ )
 1278         {
 1279             $sum = $b_vector->[0][$i][0];
 1280             for ( $j = 0; $j < $i; $j++ )
 1281             {
 1282                 $sum -= $matrix->[0][$i][$j] * $xn_vector->[0][$j][0];
 1283             }
 1284             for ( $j = ($i + 1); $j < $n; $j++ )
 1285             {
 1286                 $sum -= $matrix->[0][$i][$j] * $xn_vector->[0][$j][0];
 1287             }
 1288             $xn_vector->[0][$i][0] = $weight * ( $sum / $matrix->[0][$i][$i] )
 1289                                    + (1 - $weight) * $xn_vector->[0][$i][0];
 1290         }
 1291         $x0_vector->subtract($x0_vector,$xn_vector);
 1292         if ($norm > 0) { $diff = $x0_vector->norm_one(); }
 1293         else           { $diff = $x0_vector->norm_max(); }
 1294         for ( $i = 0; $i < $n; $i++ )
 1295         {
 1296             $x0_vector->[0][$i][0] = $xn_vector->[0][$i][0];
 1297         }
 1298     }
 1299     return($xn_vector);
 1300 }
 1301 
 1302 # Core householder reduction routine (when eagenvector
 1303 # are wanted).
 1304 # Adapted from: Numerical Recipes, 2nd edition.
 1305 sub _householder_vectors ($)
 1306 {
 1307     my ($Q) = @_;
 1308     my ($rows, $cols) = ($Q->[1], $Q->[2]);
 1309 
 1310     # Creates tridiagonal
 1311     # Set up tridiagonal needed elements
 1312     my $d = []; # N Diagonal elements 0...n-1
 1313     my $e = []; # N-1 Off-Diagonal elements 0...n-2
 1314 
 1315     my @p = ();
 1316     for (my $i = ($rows-1); $i > 1; $i--)
 1317     {
 1318   my $scale = 0.0;
 1319   # Computes norm of one column (below diagonal)
 1320   for (my $k = 0; $k < $i; $k++)
 1321   {
 1322       $scale += abs($Q->[0][$i][$k]);
 1323   }
 1324   if ($scale == 0.0)
 1325   { # skip the transformation
 1326       $e->[$i-1] = $Q->[0][$i][$i-1];
 1327   }
 1328   else
 1329   {
 1330       my $h = 0.0;
 1331       for (my $k = 0; $k < $i; $k++)
 1332       { # Used scaled Q for transformation
 1333     $Q->[0][$i][$k] /= $scale;
 1334     # Form sigma in h
 1335     $h += $Q->[0][$i][$k] * $Q->[0][$i][$k];
 1336       }
 1337       my $t1 = $Q->[0][$i][$i-1];
 1338       my $t2 = (($t1 >= 0.0) ? -sqrt($h) : sqrt($h));
 1339       $e->[$i-1] = $scale * $t2; # Update off-diagonals
 1340       $h -= $t1 * $t2;
 1341       $Q->[0][$i][$i-1] -= $t2;
 1342       my $f = 0.0;
 1343       for (my $j = 0; $j < $i; $j++)
 1344       {
 1345     $Q->[0][$j][$i] = $Q->[0][$i][$j] / $h;
 1346     my $g = 0.0;
 1347     for (my $k = 0; $k <= $j; $k++)
 1348     {
 1349         $g += $Q->[0][$j][$k] * $Q->[0][$i][$k];
 1350     }
 1351     for (my $k = $j+1; $k < $i; $k++)
 1352     {
 1353         $g += $Q->[0][$k][$j] * $Q->[0][$i][$k];
 1354     }
 1355     # Form elements of P
 1356     $p[$j] = $g / $h;
 1357     $f += $p[$j] * $Q->[0][$i][$j];
 1358       }
 1359       my $hh = $f / ($h + $h);
 1360       for (my $j = 0; $j < $i; $j++)
 1361       {
 1362     my $t3 = $Q->[0][$i][$j];
 1363     my $t4 = $p[$j] - $hh * $t3;
 1364     $p[$j] = $t4;
 1365     for (my $k = 0; $k <= $j; $k++)
 1366     {
 1367         $Q->[0][$j][$k] -= $t3 * $p[$k]
 1368       + $t4 * $Q->[0][$i][$k];
 1369     }
 1370       }
 1371   }
 1372     }
 1373     # Updates for i == 0,1
 1374     $e->[0] = $Q->[0][1][0];
 1375     $d->[0] = $Q->[0][0][0]; # i==0
 1376     $Q->[0][0][0] = 1.0;
 1377     $d->[1] = $Q->[0][1][1]; # i==1
 1378     $Q->[0][1][1] = 1.0;
 1379     $Q->[0][1][0] = $Q->[0][0][1] = 0.0;
 1380     for (my $i = 2; $i < $rows; $i++)
 1381     {
 1382   for (my $j = 0; $j < $i; $j++)
 1383   {
 1384       my $g = 0.0;
 1385       for (my $k = 0; $k < $i; $k++)
 1386       {
 1387     $g += $Q->[0][$i][$k] * $Q->[0][$k][$j];
 1388       }
 1389       for (my $k = 0; $k < $i; $k++)
 1390       {
 1391     $Q->[0][$k][$j] -= $g * $Q->[0][$k][$i];
 1392       }
 1393   }
 1394   $d->[$i] = $Q->[0][$i][$i];
 1395   # Reset row and column of Q for next iteration
 1396   $Q->[0][$i][$i] = 1.0;
 1397   for (my $j = 0; $j < $i; $j++)
 1398   {
 1399       $Q->[0][$i][$j] = $Q->[0][$j][$i] = 0.0;
 1400   }
 1401     }
 1402     return ($d, $e);
 1403 }
 1404 
 1405 # Computes sqrt(a*a + b*b), but more carefully...
 1406 sub _pythag ($$)
 1407 {
 1408     my ($a, $b) = @_;
 1409     my $aa = abs($a);
 1410     my $ab = abs($b);
 1411     if ($aa > $ab)
 1412     {
 1413   # NB: Not needed!: return 0.0 if ($aa == 0.0);
 1414   my $t = $ab / $aa;
 1415   return ($aa * sqrt(1.0 + $t*$t));
 1416     }
 1417     else
 1418     {
 1419   return 0.0 if ($ab == 0.0);
 1420   my $t = $aa / $ab;
 1421   return ($ab * sqrt(1.0 + $t*$t));
 1422     }
 1423 }
 1424 
 1425 # QL algorithm with implicit shifts to determine the eigenvalues
 1426 # of a tridiagonal matrix. Internal routine.
 1427 sub _tridiagonal_QLimplicit
 1428 {
 1429     my ($EV, $d, $e) = @_;
 1430     my ($rows, $cols) = ($EV->[1], $EV->[2]);
 1431 
 1432     $e->[$rows-1] = 0.0;
 1433     # Start real computation
 1434     for (my $l = 0; $l < $rows; $l++)
 1435     {
 1436   my $iter = 0;
 1437   my $m;
 1438   OUTER:
 1439   do {
 1440       for ($m = $l; $m < ($rows - 1); $m++)
 1441       {
 1442     my $dd = abs($d->[$m]) + abs($d->[$m+1]);
 1443     last if ((abs($e->[$m]) + $dd) == $dd);
 1444       }
 1445       if ($m != $l)
 1446       {
 1447     croak("Too many iterations!") if ($iter++ >= 30);
 1448     my $g = ($d->[$l+1] - $d->[$l])
 1449         / (2.0 * $e->[$l]);
 1450     my $r = _pythag($g, 1.0);
 1451     $g = $d->[$m] - $d->[$l]
 1452         + $e->[$l] / ($g + (($g >= 0.0) ? abs($r) : -abs($r)));
 1453     my ($p,$s,$c) = (0.0, 1.0,1.0);
 1454     for (my $i = ($m-1); $i >= $l; $i--)
 1455     {
 1456         my $ii = $i + 1;
 1457         my $f = $s * $e->[$i];
 1458         my $t = _pythag($f, $g);
 1459         $e->[$ii] = $t;
 1460         if ($t == 0.0)
 1461         {
 1462       $d->[$ii] -= $p;
 1463       $e->[$m] = 0.0;
 1464       next OUTER;
 1465         }
 1466         my $b = $c * $e->[$i];
 1467         $s = $f / $t;
 1468         $c = $g / $t;
 1469         $g = $d->[$ii] - $p;
 1470         my $t2 = ($d->[$i] - $g) * $s + 2.0 * $c * $b;
 1471         $p = $s * $t2;
 1472         $d->[$ii] = $g + $p;
 1473         $g = $c * $t2 - $b;
 1474         for (my $k = 0; $k < $rows; $k++)
 1475         {
 1476       my $t1 = $EV->[0][$k][$ii];
 1477       my $t2 = $EV->[0][$k][$i];
 1478       $EV->[0][$k][$ii] = $s * $t2 + $c * $t1;
 1479       $EV->[0][$k][$i] = $c * $t2 - $s * $t1;
 1480         }
 1481     }
 1482     $d->[$l] -= $p;
 1483     $e->[$l] = $g;
 1484     $e->[$m] = 0.0;
 1485       }
 1486   } while ($m != $l);
 1487     }
 1488     return;
 1489 }
 1490 
 1491 # Core householder reduction routine (when eagenvector
 1492 # are NOT wanted).
 1493 sub _householder_values ($)
 1494 {
 1495     my ($Q) = @_; # NB: Q is destroyed on output...
 1496     my ($rows, $cols) = ($Q->[1], $Q->[2]);
 1497 
 1498     # Creates tridiagonal
 1499     # Set up tridiagonal needed elements
 1500     my $d = []; # N Diagonal elements 0...n-1
 1501     my $e = []; # N-1 Off-Diagonal elements 0...n-2
 1502 
 1503     my @p = ();
 1504     for (my $i = ($rows - 1); $i > 1; $i--)
 1505     {
 1506   my $scale = 0.0;
 1507   for (my $k = 0; $k < $i; $k++)
 1508   {
 1509       $scale += abs($Q->[0][$i][$k]);
 1510   }
 1511   if ($scale == 0.0)
 1512   { # skip the transformation
 1513       $e->[$i-1] = $Q->[0][$i][$i-1];
 1514   }
 1515   else
 1516   {
 1517       my $h = 0.0;
 1518       for (my $k = 0; $k < $i; $k++)
 1519       { # Used scaled Q for transformation
 1520     $Q->[0][$i][$k] /= $scale;
 1521     # Form sigma in h
 1522     $h += $Q->[0][$i][$k] * $Q->[0][$i][$k];
 1523       }
 1524       my $t = $Q->[0][$i][$i-1];
 1525       my $t2 = (($t >= 0.0) ? -sqrt($h) : sqrt($h));
 1526       $e->[$i-1] = $scale * $t2; # Updates off-diagonal
 1527       $h -= $t * $t2;
 1528       $Q->[0][$i][$i-1] -= $t2;
 1529       my $f = 0.0;
 1530       for (my $j = 0; $j < $i; $j++)
 1531       {
 1532     my $g = 0.0;
 1533     for (my $k = 0; $k <= $j; $k++)
 1534     {
 1535         $g += $Q->[0][$j][$k] * $Q->[0][$i][$k];
 1536     }
 1537     for (my $k = $j+1; $k < $i; $k++)
 1538     {
 1539         $g += $Q->[0][$k][$j] * $Q->[0][$i][$k];
 1540     }
 1541     # Form elements of P
 1542     $p[$j] = $g / $h;
 1543     $f += $p[$j] * $Q->[0][$i][$j];
 1544       }
 1545       my $hh = $f / ($h + $h);
 1546       for (my $j = 0; $j < $i; $j++)
 1547       {
 1548     my $t = $Q->[0][$i][$j];
 1549     my $g = $p[$j] - $hh * $t;
 1550     $p[$j] = $g;
 1551     for (my $k = 0; $k <= $j; $k++)
 1552     {
 1553         $Q->[0][$j][$k] -= $t * $p[$k]
 1554       + $g * $Q->[0][$i][$k];
 1555     }
 1556       }
 1557   }
 1558     }
 1559     # Updates for i==1
 1560     $e->[0] =  $Q->[0][1][0];
 1561     # Updates diagonal elements
 1562     for (my $i = 0; $i < $rows; $i++)
 1563     {
 1564   $d->[$i] =  $Q->[0][$i][$i];
 1565     }
 1566     return ($d, $e);
 1567 }
 1568 
 1569 # QL algorithm with implicit shifts to determine the
 1570 # eigenvalues ONLY. This is O(N^2) only...
 1571 sub _tridiagonal_QLimplicit_values
 1572 {
 1573     my ($M, $d, $e) = @_; # NB: M is not touched...
 1574     my ($rows, $cols) = ($M->[1], $M->[2]);
 1575 
 1576     $e->[$rows-1] = 0.0;
 1577     # Start real computation
 1578     for (my $l = 0; $l < $rows; $l++)
 1579     {
 1580   my $iter = 0;
 1581   my $m;
 1582   OUTER:
 1583   do {
 1584       for ($m = $l; $m < ($rows - 1); $m++)
 1585       {
 1586     my $dd = abs($d->[$m]) + abs($d->[$m+1]);
 1587     last if ((abs($e->[$m]) + $dd) == $dd);
 1588       }
 1589       if ($m != $l)
 1590       {
 1591     croak("Too many iterations!") if ($iter++ >= 30);
 1592     my $g = ($d->[$l+1] - $d->[$l])
 1593         / (2.0 * $e->[$l]);
 1594     my $r = _pythag($g, 1.0);
 1595     $g = $d->[$m] - $d->[$l]
 1596         + $e->[$l] / ($g + (($g >= 0.0) ? abs($r) : -abs($r)));
 1597     my ($p,$s,$c) = (0.0, 1.0,1.0);
 1598     for (my $i = ($m-1); $i >= $l; $i--)
 1599     {
 1600         my $ii = $i + 1;
 1601         my $f = $s * $e->[$i];
 1602         my $t = _pythag($f, $g);
 1603         $e->[$ii] = $t;
 1604         if ($t == 0.0)
 1605         {
 1606       $d->[$ii] -= $p;
 1607       $e->[$m] = 0.0;
 1608       next OUTER;
 1609         }
 1610         my $b = $c * $e->[$i];
 1611         $s = $f / $t;
 1612         $c = $g / $t;
 1613         $g = $d->[$ii] - $p;
 1614         my $t2 = ($d->[$i] - $g) * $s + 2.0 * $c * $b;
 1615         $p = $s * $t2;
 1616         $d->[$ii] = $g + $p;
 1617         $g = $c * $t2 - $b;
 1618     }
 1619     $d->[$l] -= $p;
 1620     $e->[$l] = $g;
 1621     $e->[$m] = 0.0;
 1622       }
 1623   } while ($m != $l);
 1624     }
 1625     return;
 1626 }
 1627 
 1628 # Householder reduction of a real, symmetric matrix A.
 1629 # Returns a tridiagonal matrix T and the orthogonal matrix
 1630 # Q effecting the transformation between A and T.
 1631 sub householder ($)
 1632 {
 1633     my ($A) = @_;
 1634     my ($rows, $cols) = ($A->[1], $A->[2]);
 1635 
 1636     croak "Matrix is not quadratic"
 1637   unless ($rows = $cols);
 1638     croak "Matrix is not symmetric"
 1639         unless ($A->is_symmetric());
 1640 
 1641     # Copy given matrix TODO: study if we should do in-place modification
 1642     my $Q = $A->clone();
 1643 
 1644     # Do the computation of tridiagonal elements and of
 1645     # transformation matrix
 1646     my ($diag, $offdiag) = $Q->_householder_vectors();
 1647 
 1648     # Creates the tridiagonal matrix
 1649     my $T = $A->shadow();
 1650     for (my $i = 0; $i < $rows; $i++)
 1651     { # Set diagonal
 1652   $T->[0][$i][$i] = $diag->[$i];
 1653     }
 1654     for (my $i = 0; $i < ($rows-1); $i++)
 1655     { # Set off diagonals
 1656   $T->[0][$i+1][$i] = $offdiag->[$i];
 1657   $T->[0][$i][$i+1] = $offdiag->[$i];
 1658     }
 1659     return ($T, $Q);
 1660 }
 1661 
 1662 # QL algorithm with implicit shifts to determine the eigenvalues
 1663 # and eigenvectors of a real tridiagonal matrix - or of a matrix
 1664 # previously reduced to tridiagonal form.
 1665 sub tri_diagonalize ($;$)
 1666 {
 1667   my ($T,$Q) = @_; # Q may be 0 if the original matrix is really tridiagonal
 1668 
 1669   my ($rows, $cols) = ($T->[1], $T->[2]);
 1670 
 1671   croak "Matrix is not quadratic"
 1672     unless ($rows = $cols);
 1673   croak "Matrix is not tridiagonal"
 1674         unless ($T->is_symmetric()); # TODO: Do real tridiag check (not symmetric)!
 1675 
 1676   my $EV;
 1677   # Obtain/Creates the todo eigenvectors matrix
 1678   if ($Q)
 1679     {
 1680       $EV = $Q->clone();
 1681     }
 1682   else
 1683     {
 1684       $EV = $T->shadow();
 1685       $EV->one();
 1686     }
 1687   # Allocates diagonal vector
 1688   my $diag = [ ];
 1689   # Initializes it with T
 1690   for (my $i = 0; $i < $rows; $i++)
 1691     {
 1692       $diag->[$i] = $T->[0][$i][$i];
 1693     }
 1694   # Allocate temporary vector for off-diagonal elements
 1695   my $offdiag = [ ];
 1696   for (my $i = 1; $i < $rows; $i++)
 1697     {
 1698       $offdiag->[$i-1] = $T->[0][$i][$i-1];
 1699     }
 1700 
 1701   # Calls the calculus routine
 1702   $EV->_tridiagonal_QLimplicit($diag, $offdiag);
 1703 
 1704   # Allocate eigenvalues vector
 1705   my $v = MatrixReal1->new($rows,1);
 1706   # Fills it
 1707   for (my $i = 0; $i < $rows; $i++)
 1708     {
 1709   $v->[0][$i][0] = $diag->[$i];
 1710     }
 1711   return ($v, $EV);
 1712 }
 1713 
 1714 # Main routine for diagonalization of a real symmetric
 1715 # matrix M. Operates by transforming M into a tridiagonal
 1716 # matrix and then obtaining the eigenvalues and eigenvectors
 1717 # for that matrix (taking into account the transformation to
 1718 # tridiagonal).
 1719 sub sym_diagonalize ($)
 1720 {
 1721     my ($M) = @_;
 1722     my ($rows, $cols) = ($M->[1], $M->[2]);
 1723 
 1724     croak "Matrix is not quadratic"
 1725   unless ($rows = $cols);
 1726     croak "Matrix is not symmetric"
 1727         unless ($M->is_symmetric());
 1728 
 1729     # Copy initial matrix
 1730     # TODO: study if we should allow in-place modification
 1731     my $VEC = $M->clone();
 1732 
 1733     # Do the computation of tridiagonal elements and of
 1734     # transformation matrix
 1735     my ($diag, $offdiag) = $VEC->_householder_vectors();
 1736     # Calls the calculus routine for diagonalization
 1737     $VEC->_tridiagonal_QLimplicit($diag, $offdiag);
 1738 
 1739     # Allocate eigenvalues vector
 1740     my $val = MatrixReal1->new($rows,1);
 1741     # Fills it
 1742     for (my $i = 0; $i < $rows; $i++)
 1743     {
 1744   $val->[0][$i][0] = $diag->[$i];
 1745     }
 1746     return ($val, $VEC);
 1747 }
 1748 
 1749 # Householder reduction of a real, symmetric matrix A.
 1750 # Returns a tridiagonal matrix T equivalent to A.
 1751 sub householder_tridiagonal ($)
 1752 {
 1753     my ($A) = @_;
 1754     my ($rows, $cols) = ($A->[1], $A->[2]);
 1755 
 1756     croak "Matrix is not quadratic"
 1757   unless ($rows = $cols);
 1758     croak "Matrix is not symmetric"
 1759         unless ($A->is_symmetric());
 1760 
 1761     # Copy given matrix
 1762     my $Q = $A->clone();
 1763 
 1764     # Do the computation of tridiagonal elements and of
 1765     # transformation matrix
 1766     # Q is destroyed after reduction
 1767     my ($diag, $offdiag) = $Q->_householder_values();
 1768 
 1769     # Creates the tridiagonal matrix in Q (avoid allocation)
 1770     my $T = $Q;
 1771     $T->zero();
 1772     for (my $i = 0; $i < $rows; $i++)
 1773     { # Set diagonal
 1774   $T->[0][$i][$i] = $diag->[$i];
 1775     }
 1776     for (my $i = 0; $i < ($rows-1); $i++)
 1777     { # Set off diagonals
 1778   $T->[0][$i+1][$i] = $offdiag->[$i];
 1779   $T->[0][$i][$i+1] = $offdiag->[$i];
 1780     }
 1781     return $T;
 1782 }
 1783 
 1784 # QL algorithm with implicit shifts to determine ONLY
 1785 # the eigenvalues a real tridiagonal matrix - or of a
 1786 # matrix previously reduced to tridiagonal form.
 1787 sub tri_eigenvalues ($;$)
 1788 {
 1789   my ($T) = @_;
 1790   my ($rows, $cols) = ($T->[1], $T->[2]);
 1791 
 1792   croak "Matrix is not quadratic"
 1793     unless ($rows = $cols);
 1794   croak "Matrix is not tridiagonal"
 1795         unless ($T->is_symmetric()); # TODO: Do real tridiag check (not symmetric)!
 1796 
 1797   # Allocates diagonal vector
 1798   my $diag = [ ];
 1799   # Initializes it with T
 1800   for (my $i = 0; $i < $rows; $i++)
 1801     {
 1802       $diag->[$i] = $T->[0][$i][$i];
 1803     }
 1804   # Allocate temporary vector for off-diagonal elements
 1805   my $offdiag = [ ];
 1806   for (my $i = 1; $i < $rows; $i++)
 1807     {
 1808       $offdiag->[$i-1] = $T->[0][$i][$i-1];
 1809     }
 1810 
 1811   # Calls the calculus routine (T is not touched)
 1812   $T->_tridiagonal_QLimplicit_values($diag, $offdiag);
 1813 
 1814   # Allocate eigenvalues vector
 1815   my $v = MatrixReal1->new($rows,1);
 1816   # Fills it
 1817   for (my $i = 0; $i < $rows; $i++)
 1818     {
 1819   $v->[0][$i][0] = $diag->[$i];
 1820     }
 1821   return $v;
 1822 }
 1823 
 1824 # Main routine for diagonalization of a real symmetric
 1825 # matrix M. Operates by transforming M into a tridiagonal
 1826 # matrix and then obtaining the eigenvalues and eigenvectors
 1827 # for that matrix (taking into account the transformation to
 1828 # tridiagonal).
 1829 sub sym_eigenvalues ($)
 1830 {
 1831     my ($M) = @_;
 1832     my ($rows, $cols) = ($M->[1], $M->[2]);
 1833 
 1834     croak "Matrix is not quadratic"
 1835   unless ($rows = $cols);
 1836     croak "Matrix is not symmetric"
 1837         unless ($M->is_symmetric());
 1838 
 1839     # Copy matrix in temporary
 1840     my $A = $M->clone();
 1841     # Do the computation of tridiagonal elements and of
 1842     # transformation matrix. A is destroyed
 1843     my ($diag, $offdiag) = $A->_householder_values();
 1844     # Calls the calculus routine for diagonalization
 1845     # (M is not touched)
 1846     $M->_tridiagonal_QLimplicit_values($diag, $offdiag);
 1847 
 1848     # Allocate eigenvalues vector
 1849     my $val = MatrixReal1->new($rows,1);
 1850     # Fills it
 1851     for (my $i = 0; $i < $rows; $i++)
 1852     {
 1853   $val->[0][$i][0] = $diag->[$i];
 1854     }
 1855     return $val;
 1856 }
 1857 
 1858 # Boolean check routine to see if a matrix is
 1859 # symmetric
 1860 sub is_symmetric ($)
 1861 {
 1862   my ($M) = @_;
 1863   my ($rows, $cols) = ($M->[1], $M->[2]);
 1864   # if it is not quadratic it cannot be symmetric...
 1865   return 0 unless ($rows == $cols);
 1866   for (my $i = 1; $i < $rows; $i++)
 1867     {
 1868       for (my $j = 0; $j < $i; $j++)
 1869   {
 1870     return 0 unless ($M->[0][$i][$j] == $M->[0][$j][$i]);
 1871   }
 1872     }
 1873   return 1;
 1874 }
 1875 
 1876                 ########################################
 1877                 #                                      #
 1878                 # define overloaded operators section: #
 1879                 #                                      #
 1880                 ########################################
 1881 
 1882 sub _negate
 1883 {
 1884     my($object,$argument,$flag) = @_;
 1885 #   my($name) = "neg"; #&_trace($name,$object,$argument,$flag);
 1886     my($temp);
 1887 
 1888     $temp = $object->new($object->[1],$object->[2]);
 1889     $temp->negate($object);
 1890     return($temp);
 1891 }
 1892 
 1893 sub _transpose
 1894 {
 1895     my($object,$argument,$flag) = @_;
 1896 #   my($name) = "'~'"; #&_trace($name,$object,$argument,$flag);
 1897     my($temp);
 1898 
 1899     $temp = $object->new($object->[2],$object->[1]);
 1900     $temp->transpose($object);
 1901     return($temp);
 1902 }
 1903 
 1904 sub _boolean
 1905 {
 1906     my($object,$argument,$flag) = @_;
 1907 #   my($name) = "bool"; #&_trace($name,$object,$argument,$flag);
 1908     my($rows,$cols) = ($object->[1],$object->[2]);
 1909     my($i,$j,$result);
 1910 
 1911     $result = 0;
 1912     BOOL:
 1913     for ( $i = 0; $i < $rows; $i++ )
 1914     {
 1915         for ( $j = 0; $j < $cols; $j++ )
 1916         {
 1917             if ($object->[0][$i][$j] != 0)
 1918             {
 1919                 $result = 1;
 1920                 last BOOL;
 1921             }
 1922         }
 1923     }
 1924     return($result);
 1925 }
 1926 
 1927 sub _not_boolean
 1928 {
 1929     my($object,$argument,$flag) = @_;
 1930 #   my($name) = "'!'"; #&_trace($name,$object,$argument,$flag);
 1931     my($rows,$cols) = ($object->[1],$object->[2]);
 1932     my($i,$j,$result);
 1933 
 1934     $result = 1;
 1935     NOTBOOL:
 1936     for ( $i = 0; $i < $rows; $i++ )
 1937     {
 1938         for ( $j = 0; $j < $cols; $j++ )
 1939         {
 1940             if ($object->[0][$i][$j] != 0)
 1941             {
 1942                 $result = 0;
 1943                 last NOTBOOL;
 1944             }
 1945         }
 1946     }
 1947     return($result);
 1948 }
 1949 
 1950 sub _stringify
 1951 {
 1952     my($object,$argument,$flag) = @_;
 1953 #   my($name) = '""'; #&_trace($name,$object,$argument,$flag);
 1954     my($rows,$cols) = ($object->[1],$object->[2]);
 1955     my($i,$j,$s);
 1956 
 1957     $s = '';
 1958     for ( $i = 0; $i < $rows; $i++ )
 1959     {
 1960         $s .= "[ ";
 1961         for ( $j = 0; $j < $cols; $j++ )
 1962         {
 1963             $s .= sprintf("% #-19.12E ", $object->[0][$i][$j]);
 1964         }
 1965         $s .= "]\n";
 1966     }
 1967     return($s);
 1968 }
 1969 
 1970 sub _norm
 1971 {
 1972     my($object,$argument,$flag) = @_;
 1973 #   my($name) = "abs"; #&_trace($name,$object,$argument,$flag);
 1974 
 1975     return( $object->norm_one() );
 1976 }
 1977 
 1978 sub _add
 1979 {
 1980     my($object,$argument,$flag) = @_;
 1981     my($name) = "'+'"; #&_trace($name,$object,$argument,$flag);
 1982     my($temp);
 1983 
 1984     if ((defined $argument) && ref($argument) &&
 1985         (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
 1986     {
 1987         if (defined $flag)
 1988         {
 1989             $temp = $object->new($object->[1],$object->[2]);
 1990             $temp->add($object,$argument);
 1991             return($temp);
 1992         }
 1993         else
 1994         {
 1995             $object->add($object,$argument);
 1996             return($object);
 1997         }
 1998     }
 1999     else
 2000     {
 2001         croak "MatrixReal1 $name: wrong argument type";
 2002     }
 2003 }
 2004 
 2005 sub _subtract
 2006 {
 2007     my($object,$argument,$flag) = @_;
 2008     my($name) = "'-'"; #&_trace($name,$object,$argument,$flag);
 2009     my($temp);
 2010 
 2011     if ((defined $argument) && ref($argument) &&
 2012         (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
 2013     {
 2014         if (defined $flag)
 2015         {
 2016             $temp = $object->new($object->[1],$object->[2]);
 2017             if ($flag) { $temp->subtract($argument,$object); }
 2018             else       { $temp->subtract($object,$argument); }
 2019             return($temp);
 2020         }
 2021         else
 2022         {
 2023             $object->subtract($object,$argument);
 2024             return($object);
 2025         }
 2026     }
 2027     else
 2028     {
 2029         croak "MatrixReal1 $name: wrong argument type";
 2030     }
 2031 }
 2032 
 2033 sub _multiply
 2034 {
 2035     my($object,$argument,$flag) = @_;
 2036     my($name) = "'*'"; #&_trace($name,$object,$argument,$flag);
 2037     my($temp);
 2038 
 2039     if ((defined $argument) && ref($argument) &&
 2040         (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
 2041     {
 2042         if ((defined $flag) && $flag)
 2043         {
 2044             return( multiply($argument,$object) );
 2045         }
 2046         else
 2047         {
 2048             return( multiply($object,$argument) );
 2049         }
 2050     }
 2051     elsif ((defined $argument) && !(ref($argument)))
 2052     {
 2053         if (defined $flag)
 2054         {
 2055             $temp = $object->new($object->[1],$object->[2]);
 2056             $temp->multiply_scalar($object,$argument);
 2057             return($temp);
 2058         }
 2059         else
 2060         {
 2061             $object->multiply_scalar($object,$argument);
 2062             return($object);
 2063         }
 2064     }
 2065     else
 2066     {
 2067         croak "MatrixReal1 $name: wrong argument type";
 2068     }
 2069 }
 2070 
 2071 sub _assign_add
 2072 {
 2073     my($object,$argument,$flag) = @_;
 2074 #   my($name) = "'+='"; #&_trace($name,$object,$argument,$flag);
 2075 
 2076     return( &_add($object,$argument,undef) );
 2077 }
 2078 
 2079 sub _assign_subtract
 2080 {
 2081     my($object,$argument,$flag) = @_;
 2082 #   my($name) = "'-='"; #&_trace($name,$object,$argument,$flag);
 2083 
 2084     return( &_subtract($object,$argument,undef) );
 2085 }
 2086 
 2087 sub _assign_multiply
 2088 {
 2089     my($object,$argument,$flag) = @_;
 2090 #   my($name) = "'*='"; #&_trace($name,$object,$argument,$flag);
 2091 
 2092     return( &_multiply($object,$argument,undef) );
 2093 }
 2094 
 2095 sub _equal
 2096 {
 2097     my($object,$argument,$flag) = @_;
 2098     my($name) = "'=='"; #&_trace($name,$object,$argument,$flag);
 2099     my($rows,$cols) = ($object->[1],$object->[2]);
 2100     my($i,$j,$result);
 2101 
 2102     if ((defined $argument) && ref($argument) &&
 2103         (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
 2104     {
 2105         $result = 1;
 2106         EQUAL:
 2107         for ( $i = 0; $i < $rows; $i++ )
 2108         {
 2109             for ( $j = 0; $j < $cols; $j++ )
 2110             {
 2111                 if ($object->[0][$i][$j] != $argument->[0][$i][$j])
 2112                 {
 2113                     $result = 0;
 2114                     last EQUAL;
 2115                 }
 2116             }
 2117         }
 2118         return($result);
 2119     }
 2120     else
 2121     {
 2122         croak "MatrixReal1 $name: wrong argument type";
 2123     }
 2124 }
 2125 
 2126 sub _not_equal
 2127 {
 2128     my($object,$argument,$flag) = @_;
 2129     my($name) = "'!='"; #&_trace($name,$object,$argument,$flag);
 2130     my($rows,$cols) = ($object->[1],$object->[2]);
 2131     my($i,$j,$result);
 2132 
 2133     if ((defined $argument) && ref($argument) &&
 2134         (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
 2135     {
 2136         $result = 0;
 2137         NOTEQUAL:
 2138         for ( $i = 0; $i < $rows; $i++ )
 2139         {
 2140             for ( $j = 0; $j < $cols; $j++ )
 2141             {
 2142                 if ($object->[0][$i][$j] != $argument->[0][$i][$j])
 2143                 {
 2144                     $result = 1;
 2145                     last NOTEQUAL;
 2146                 }
 2147             }
 2148         }
 2149         return($result);
 2150     }
 2151     else
 2152     {
 2153         croak "MatrixReal1 $name: wrong argument type";
 2154     }
 2155 }
 2156 
 2157 sub _less_than
 2158 {
 2159     my($object,$argument,$flag) = @_;
 2160     my($name) = "'<'"; #&_trace($name,$object,$argument,$flag);
 2161 
 2162     if ((defined $argument) && ref($argument) &&
 2163         (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
 2164     {
 2165         if ((defined $flag) && $flag)
 2166         {
 2167             return( $argument->norm_one() < $object->norm_one() );
 2168         }
 2169         else
 2170         {
 2171             return( $object->norm_one() < $argument->norm_one() );
 2172         }
 2173     }
 2174     elsif ((defined $argument) && !(ref($argument)))
 2175     {
 2176         if ((defined $flag) && $flag)
 2177         {
 2178             return( abs($argument) < $object->norm_one() );
 2179         }
 2180         else
 2181         {
 2182             return( $object->norm_one() < abs($argument) );
 2183         }
 2184     }
 2185     else
 2186     {
 2187         croak "MatrixReal1 $name: wrong argument type";
 2188     }
 2189 }
 2190 
 2191 sub _less_than_or_equal
 2192 {
 2193     my($object,$argument,$flag) = @_;
 2194     my($name) = "'<='"; #&_trace($name,$object,$argument,$flag);
 2195 
 2196     if ((defined $argument) && ref($argument) &&
 2197         (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
 2198     {
 2199         if ((defined $flag) && $flag)
 2200         {
 2201             return( $argument->norm_one() <= $object->norm_one() );
 2202         }
 2203         else
 2204         {
 2205             return( $object->norm_one() <= $argument->norm_one() );
 2206         }
 2207     }
 2208     elsif ((defined $argument) && !(ref($argument)))
 2209     {
 2210         if ((defined $flag) && $flag)
 2211         {
 2212             return( abs($argument) <= $object->norm_one() );
 2213         }
 2214         else
 2215         {
 2216             return( $object->norm_one() <= abs($argument) );
 2217         }
 2218     }
 2219     else
 2220     {
 2221         croak "MatrixReal1 $name: wrong argument type";
 2222     }
 2223 }
 2224 
 2225 sub _greater_than
 2226 {
 2227     my($object,$argument,$flag) = @_;
 2228     my($name) = "'>'"; #&_trace($name,$object,$argument,$flag);
 2229 
 2230     if ((defined $argument) && ref($argument) &&
 2231         (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
 2232     {
 2233         if ((defined $flag) && $flag)
 2234         {
 2235             return( $argument->norm_one() > $object->norm_one() );
 2236         }
 2237         else
 2238         {
 2239             return( $object->norm_one() > $argument->norm_one() );
 2240         }
 2241     }
 2242     elsif ((defined $argument) && !(ref($argument)))
 2243     {
 2244         if ((defined $flag) && $flag)
 2245         {
 2246             return( abs($argument) > $object->norm_one() );
 2247         }
 2248         else
 2249         {
 2250             return( $object->norm_one() > abs($argument) );
 2251         }
 2252     }
 2253     else
 2254     {
 2255         croak "MatrixReal1 $name: wrong argument type";
 2256     }
 2257 }
 2258 
 2259 sub _greater_than_or_equal
 2260 {
 2261     my($object,$argument,$flag) = @_;
 2262     my($name) = "'>='"; #&_trace($name,$object,$argument,$flag);
 2263 
 2264     if ((defined $argument) && ref($argument) &&
 2265         (ref($argument) !~ /^SCALAR$|^ARRAY$|^HASH$|^CODE$|^REF$/))
 2266     {
 2267         if ((defined $flag) && $flag)
 2268         {
 2269             return( $argument->norm_one() >= $object->norm_one() );
 2270         }
 2271         else
 2272         {
 2273             return( $object->norm_one() >= $argument->norm_one() );
 2274         }
 2275     }
 2276     elsif ((defined $argument) && !(ref($argument)))
 2277     {
 2278         if ((defined $flag) && $flag)
 2279         {
 2280             return( abs($argument) >= $object->norm_one() );
 2281         }
 2282         else
 2283         {
 2284             return( $object->norm_one() >= abs($argument) );
 2285         }
 2286     }
 2287     else
 2288     {
 2289         croak "MatrixReal1 $name: wrong argument type";
 2290     }
 2291 }
 2292 
 2293 sub _clone
 2294 {
 2295     my($object,$argument,$flag) = @_;
 2296 #   my($name) = "'='"; #&_trace($name,$object,$argument,$flag);
 2297     my($temp);
 2298 
 2299     $temp = $object->new($object->[1],$object->[2]);
 2300     $temp->copy($object);
 2301     $temp->_undo_LR();
 2302     return($temp);
 2303 }
 2304 
 2305 sub _trace
 2306 {
 2307     my($text,$object,$argument,$flag) = @_;
 2308 
 2309     unless (defined $object)   { $object   = 'undef'; };
 2310     unless (defined $argument) { $argument = 'undef'; };
 2311     unless (defined $flag)     { $flag     = 'undef'; };
 2312     if (ref($object))   { $object   = ref($object);   }
 2313     if (ref($argument)) { $argument = ref($argument); }
 2314     print "$text: \$obj='$object' \$arg='$argument' \$flag='$flag'\n";
 2315 }
 2316 
 2317 1;
 2318 
 2319 __END__
 2320 
 2321 =head1 NAME
 2322 
 2323 MatrixReal1 - Matrix of Reals
 2324 
 2325 Implements the data type "matrix of reals" (and consequently also
 2326 "vector of reals")
 2327 
 2328 =head1 DESCRIPTION
 2329 
 2330 Implements the data type "matrix of reals", which can be used almost
 2331 like any other basic Perl type thanks to B<OPERATOR OVERLOADING>, i.e.,
 2332 
 2333   $product = $matrix1 * $matrix2;
 2334 
 2335 does what you would like it to do (a matrix multiplication).
 2336 
 2337 Also features many important operations and methods: matrix norm,
 2338 matrix transposition, matrix inverse, determinant of a matrix, order
 2339 and numerical condition of a matrix, scalar product of vectors, vector
 2340 product of vectors, vector length, projection of row and column vectors,
 2341 a comfortable way for reading in a matrix from a file, the keyboard or
 2342 your code, and many more.
 2343 
 2344 Allows to solve linear equation systems using an efficient algorithm
 2345 known as "L-R-decomposition" and several approximative (iterative) methods.
 2346 
 2347 Features an implementation of Kleene's algorithm to compute the minimal
 2348 costs for all paths in a graph with weighted edges (the "weights" being
 2349 the costs associated with each edge).
 2350 
 2351 =head1 SYNOPSIS
 2352 
 2353 =over 2
 2354 
 2355 =item *
 2356 
 2357 C<use MatrixReal1;>
 2358 
 2359 Makes the methods and overloaded operators of this module available
 2360 to your program.
 2361 
 2362 =item *
 2363 
 2364 C<use MatrixReal1 qw(min max);>
 2365 
 2366 =item *
 2367 
 2368 C<use MatrixReal1 qw(:all);>
 2369 
 2370 Use one of these two variants to import (all) the functions which the module
 2371 offers for export; currently these are "min()" and "max()".
 2372 
 2373 =item *
 2374 
 2375 C<$new_matrix = new MatrixReal1($rows,$columns);>
 2376 
 2377 The matrix object constructor method.
 2378 
 2379 Note that this method is implicitly called by many of the other methods
 2380 in this module!
 2381 
 2382 =item *
 2383 
 2384 C<$new_matrix = MatrixReal1-E<gt>>C<new($rows,$columns);>
 2385 
 2386 An alternate way of calling the matrix object constructor method.
 2387 
 2388 =item *
 2389 
 2390 C<$new_matrix = $some_matrix-E<gt>>C<new($rows,$columns);>
 2391 
 2392 Still another way of calling the matrix object constructor method.
 2393 
 2394 Matrix "C<$some_matrix>" is not changed by this in any way.
 2395 
 2396 =item *
 2397 
 2398 C<$new_matrix = MatrixReal1-E<gt>>C<new_from_string($string);>
 2399 
 2400 This method allows you to read in a matrix from a string (for
 2401 instance, from the keyboard, from a file or from your code).
 2402 
 2403 The syntax is simple: each row must start with "C<[ >" and end with
 2404 "C< ]\n>" ("C<\n>" being the newline character and "C< >" a space or
 2405 tab) and contain one or more numbers, all separated from each other
 2406 by spaces or tabs.
 2407 
 2408 Additional spaces or tabs can be added at will, but no comments.
 2409 
 2410 Examples:
 2411 
 2412   $string = "[ 1 2 3 ]\n[ 2 2 -1 ]\n[ 1 1 1 ]\n";
 2413   $matrix = MatrixReal1->new_from_string($string);
 2414   print "$matrix";
 2415 
 2416 By the way, this prints
 2417 
 2418   [  1.000000000000E+00  2.000000000000E+00  3.000000000000E+00 ]
 2419   [  2.000000000000E+00  2.000000000000E+00 -1.000000000000E+00 ]
 2420   [  1.000000000000E+00  1.000000000000E+00  1.000000000000E+00 ]
 2421 
 2422 But you can also do this in a much more comfortable way using the
 2423 shell-like "here-document" syntax:
 2424 
 2425   $matrix = MatrixReal1->new_from_string(<<'MATRIX');
 2426   [  1  0  0  0  0  0  1  ]
 2427   [  0  1  0  0  0  0  0  ]
 2428   [  0  0  1  0  0  0  0  ]
 2429   [  0  0  0  1  0  0  0  ]
 2430   [  0  0  0  0  1  0  0  ]
 2431   [  0  0  0  0  0  1  0  ]
 2432   [  1  0  0  0  0  0 -1  ]
 2433   MATRIX
 2434 
 2435 You can even use variables in the matrix:
 2436 
 2437   $c1 =   2  /  3;
 2438   $c2 =  -2  /  5;
 2439   $c3 =  26  /  9;
 2440 
 2441   $matrix = MatrixReal1->new_from_string(<<"MATRIX");
 2442 
 2443       [   3    2    0   ]
 2444       [   0    3    2   ]
 2445       [  $c1  $c2  $c3  ]
 2446 
 2447   MATRIX
 2448 
 2449 (Remember that you may use spaces and tabs to format the matrix to
 2450 your taste)
 2451 
 2452 Note that this method uses exactly the same representation for a
 2453 matrix as the "stringify" operator "": this means that you can convert
 2454 any matrix into a string with C<$string = "$matrix";> and read it back
 2455 in later (for instance from a file!).
 2456 
 2457 Note however that you may suffer a precision loss in this process
 2458 because only 13 digits are supported in the mantissa when printed!!
 2459 
 2460 If the string you supply (or someone else supplies) does not obey
 2461 the syntax mentioned above, an exception is raised, which can be
 2462 caught by "eval" as follows:
 2463 
 2464   print "Please enter your matrix (in one line): ";
 2465   $string = <STDIN>;
 2466   $string =~ s/\\n/\n/g;
 2467   eval { $matrix = MatrixReal1->new_from_string($string); };
 2468   if ($@)
 2469   {
 2470       print "$@";
 2471       # ...
 2472       # (error handling)
 2473   }
 2474   else
 2475   {
 2476       # continue...
 2477   }
 2478 
 2479 or as follows:
 2480 
 2481   eval { $matrix = MatrixReal1->new_from_string(<<"MATRIX"); };
 2482   [   3    2    0   ]
 2483   [   0    3    2   ]
 2484   [  $c1  $c2  $c3  ]
 2485   MATRIX
 2486   if ($@)
 2487   # ...
 2488 
 2489 Actually, the method shown above for reading a matrix from the keyboard
 2490 is a little awkward, since you have to enter a lot of "\n"'s for the
 2491 newlines.
 2492 
 2493 A better way is shown in this piece of code:
 2494 
 2495   while (1)
 2496   {
 2497       print "\nPlease enter your matrix ";
 2498       print "(multiple lines, <ctrl-D> = done):\n";
 2499       eval { $new_matrix =
 2500           MatrixReal1->new_from_string(join('',<STDIN>)); };
 2501       if ($@)
 2502       {
 2503           $@ =~ s/\s+at\b.*?$//;
 2504           print "${@}Please try again.\n";
 2505       }
 2506       else { last; }
 2507   }
 2508 
 2509 Possible error messages of the "new_from_string()" method are:
 2510 
 2511   MatrixReal1::new_from_string(): syntax error in input string
 2512   MatrixReal1::new_from_string(): empty input string
 2513 
 2514 If the input string has rows with varying numbers of columns,
 2515 the following warning will be printed to STDERR:
 2516 
 2517   MatrixReal1::new_from_string(): missing elements will be set to zero!
 2518 
 2519 If everything is okay, the method returns an object reference to the
 2520 (newly allocated) matrix containing the elements you specified.
 2521 
 2522 =item *
 2523 
 2524 C<$new_matrix = $some_matrix-E<gt>shadow();>
 2525 
 2526 Returns an object reference to a B<NEW> but B<EMPTY> matrix
 2527 (filled with zero's) of the B<SAME SIZE> as matrix "C<$some_matrix>".
 2528 
 2529 Matrix "C<$some_matrix>" is not changed by this in any way.
 2530 
 2531 =item *
 2532 
 2533 C<$matrix1-E<gt>copy($matrix2);>
 2534 
 2535 Copies the contents of matrix "C<$matrix2>" to an B<ALREADY EXISTING>
 2536 matrix "C<$matrix1>" (which must have the same size as matrix "C<$matrix2>"!).
 2537 
 2538 Matrix "C<$matrix2>" is not changed by this in any way.
 2539 
 2540 =item *
 2541 
 2542 C<$twin_matrix = $some_matrix-E<gt>clone();>
 2543 
 2544 Returns an object reference to a B<NEW> matrix of the B<SAME SIZE> as
 2545 matrix "C<$some_matrix>". The contents of matrix "C<$some_matrix>" have
 2546 B<ALREADY BEEN COPIED> to the new matrix "C<$twin_matrix>".
 2547 
 2548 Matrix "C<$some_matrix>" is not changed by this in any way.
 2549 
 2550 =item *
 2551 
 2552 C<$row_vector = $matrix-E<gt>row($row);>
 2553 
 2554 This is a projection method which returns an object reference to
 2555 a B<NEW> matrix (which in fact is a (row) vector since it has only
 2556 one row) to which row number "C<$row>" of matrix "C<$matrix>" has
 2557 already been copied.
 2558 
 2559 Matrix "C<$matrix>" is not changed by this in any way.
 2560 
 2561 =item *
 2562 
 2563 C<$column_vector = $matrix-E<gt>column($column);>
 2564 
 2565 This is a projection method which returns an object reference to
 2566 a B<NEW> matrix (which in fact is a (column) vector since it has
 2567 only one column) to which column number "C<$column>" of matrix
 2568 "C<$matrix>" has already been copied.
 2569 
 2570 Matrix "C<$matrix>" is not changed by this in any way.
 2571 
 2572 =item *
 2573 
 2574 C<$matrix-E<gt>zero();>
 2575 
 2576 Assigns a zero to every element of the matrix "C<$matrix>", i.e.,
 2577 erases all values previously stored there, thereby effectively
 2578 transforming the matrix into a "zero"-matrix or "null"-matrix,
 2579 the neutral element of the addition operation in a Ring.
 2580 
 2581 (For instance the (quadratic) matrices with "n" rows and columns
 2582 and matrix addition and multiplication form a Ring. Most prominent
 2583 characteristic of a Ring is that multiplication is not commutative,
 2584 i.e., in general, "C<matrix1 * matrix2>" is not the same as
 2585 "C<matrix2 * matrix1>"!)
 2586 
 2587 =item *
 2588 
 2589 C<$matrix-E<gt>one();>
 2590 
 2591 Assigns one's to the elements on the main diagonal (elements (1,1),
 2592 (2,2), (3,3) and so on) of matrix "C<$matrix>" and zero's to all others,
 2593 thereby erasing all values previously stored there and transforming the
 2594 matrix into a "one"-matrix, the neutral element of the multiplication
 2595 operation in a Ring.
 2596 
 2597 (If the matrix is quadratic (which this method doesn't require, though),
 2598 then multiplying this matrix with itself yields this same matrix again,
 2599 and multiplying it with some other matrix leaves that other matrix
 2600 unchanged!)
 2601 
 2602 =item *
 2603 
 2604 C<$matrix-E<gt>assign($row,$column,$value);>
 2605 
 2606 Explicitly assigns a value "C<$value>" to a single element of the
 2607 matrix "C<$matrix>", located in row "C<$row>" and column "C<$column>",
 2608 thereby replacing the value previously stored there.
 2609 
 2610 =item *
 2611 
 2612 C<$value = $matrix-E<gt>>C<element($row,$column);>
 2613 
 2614 Returns the value of a specific element of the matrix "C<$matrix>",
 2615 located in row "C<$row>" and column "C<$column>".
 2616 
 2617 =item *
 2618 
 2619 C<($rows,$columns) = $matrix-E<gt>dim();>
 2620 
 2621 Returns a list of two items, representing the number of rows
 2622 and columns the given matrix "C<$matrix>" contains.
 2623 
 2624 =item *
 2625 
 2626 C<$norm_one = $matrix-E<gt>norm_one();>
 2627 
 2628 Returns the "one"-norm of the given matrix "C<$matrix>".
 2629 
 2630 The "one"-norm is defined as follows:
 2631 
 2632 For each column, the sum of the absolute values of the elements in the
 2633 different rows of that column is calculated. Finally, the maximum
 2634 of these sums is returned.
 2635 
 2636 Note that the "one"-norm and the "maximum"-norm are mathematically
 2637 equivalent, although for the same matrix they usually yield a different
 2638 value.
 2639 
 2640 Therefore, you should only compare values that have been calculated
 2641 using the same norm!
 2642 
 2643 Throughout this package, the "one"-norm is (arbitrarily) used
 2644 for all comparisons, for the sake of uniformity and comparability,
 2645 except for the iterative methods "solve_GSM()", "solve_SSM()" and
 2646 "solve_RM()" which use either norm depending on the matrix itself.
 2647 
 2648 =item *
 2649 
 2650 C<$norm_max = $matrix-E<gt>norm_max();>
 2651 
 2652 Returns the "maximum"-norm of the given matrix "C<$matrix>".
 2653 
 2654 The "maximum"-norm is defined as follows:
 2655 
 2656 For each row, the sum of the absolute values of the elements in the
 2657 different columns of that row is calculated. Finally, the maximum
 2658 of these sums is returned.
 2659 
 2660 Note that the "maximum"-norm and the "one"-norm are mathematically
 2661 equivalent, although for the same matrix they usually yield a different
 2662 value.
 2663 
 2664 Therefore, you should only compare values that have been calculated
 2665 using the same norm!
 2666 
 2667 Throughout this package, the "one"-norm is (arbitrarily) used
 2668 for all comparisons, for the sake of uniformity and comparability,
 2669 except for the iterative methods "solve_GSM()", "solve_SSM()" and
 2670 "solve_RM()" which use either norm depending on the matrix itself.
 2671 
 2672 =item *
 2673 
 2674 C<$matrix1-E<gt>negate($matrix2);>
 2675 
 2676 Calculates the negative of matrix "C<$matrix2>" (i.e., multiplies
 2677 all elements with "-1") and stores the result in matrix "C<$matrix1>"
 2678 (which must already exist and have the same size as matrix "C<$matrix2>"!).
 2679 
 2680 This operation can also be carried out "in-place", i.e., input and
 2681 output matrix may be identical.
 2682 
 2683 =item *
 2684 
 2685 C<$matrix1-E<gt>transpose($matrix2);>
 2686 
 2687 Calculates the transposed matrix of matrix "C<$matrix2>" and stores
 2688 the result in matrix "C<$matrix1>" (which must already exist and have
 2689 the same size as matrix "C<$matrix2>"!).
 2690 
 2691 This operation can also be carried out "in-place", i.e., input and
 2692 output matrix may be identical.
 2693 
 2694 Transposition is a symmetry operation: imagine you rotate the matrix
 2695 along the axis of its main diagonal (going through elements (1,1),
 2696 (2,2), (3,3) and so on) by 180 degrees.
 2697 
 2698 Another way of looking at it is to say that rows and columns are
 2699 swapped. In fact the contents of element C<(i,j)> are swapped
 2700 with those of element C<(j,i)>.
 2701 
 2702 Note that (especially for vectors) it makes a big difference if you
 2703 have a row vector, like this:
 2704 
 2705   [ -1 0 1 ]
 2706 
 2707 or a column vector, like this:
 2708 
 2709   [ -1 ]
 2710   [  0 ]
 2711   [  1 ]
 2712 
 2713 the one vector being the transposed of the other!
 2714 
 2715 This is especially true for the matrix product of two vectors:
 2716 
 2717                [ -1 ]
 2718   [ -1 0 1 ] * [  0 ]  =  [ 2 ] ,  whereas
 2719                [  1 ]
 2720 
 2721                              *     [ -1  0  1 ]
 2722   [ -1 ]                                            [  1  0 -1 ]
 2723   [  0 ] * [ -1 0 1 ]  =  [ -1 ]   [  1  0 -1 ]  =  [  0  0  0 ]
 2724   [  1 ]                  [  0 ]   [  0  0  0 ]     [ -1  0  1 ]
 2725                           [  1 ]   [ -1  0  1 ]
 2726 
 2727 So be careful about what you really mean!
 2728 
 2729 Hint: throughout this module, whenever a vector is explicitly required
 2730 for input, a B<COLUMN> vector is expected!
 2731 
 2732 =item *
 2733 
 2734 C<$matrix1-E<gt>add($matrix2,$matrix3);>
 2735 
 2736 Calculates the sum of matrix "C<$matrix2>" and matrix "C<$matrix3>"
 2737 and stores the result in matrix "C<$matrix1>" (which must already exist
 2738 and have the same size as matrix "C<$matrix2>" and matrix "C<$matrix3>"!).
 2739 
 2740 This operation can also be carried out "in-place", i.e., the output and
 2741 one (or both) of the input matrices may be identical.
 2742 
 2743 =item *
 2744 
 2745 C<$matrix1-E<gt>subtract($matrix2,$matrix3);>
 2746 
 2747 Calculates the difference of matrix "C<$matrix2>" minus matrix "C<$matrix3>"
 2748 and stores the result in matrix "C<$matrix1>" (which must already exist
 2749 and have the same size as matrix "C<$matrix2>" and matrix "C<$matrix3>"!).
 2750 
 2751 This operation can also be carried out "in-place", i.e., the output and
 2752 one (or both) of the input matrices may be identical.
 2753 
 2754 Note that this operation is the same as
 2755 C<$matrix1-E<gt>add($matrix2,-$matrix3);>, although the latter is
 2756 a little less efficient.
 2757 
 2758 =item *
 2759 
 2760 C<$matrix1-E<gt>multiply_scalar($matrix2,$scalar);>
 2761 
 2762 Calculates the product of matrix "C<$matrix2>" and the number "C<$scalar>"
 2763 (i.e., multiplies each element of matrix "C<$matrix2>" with the factor
 2764 "C<$scalar>") and stores the result in matrix "C<$matrix1>" (which must
 2765 already exist and have the same size as matrix "C<$matrix2>"!).
 2766 
 2767 This operation can also be carried out "in-place", i.e., input and
 2768 output matrix may be identical.
 2769 
 2770 =item *
 2771 
 2772 C<$product_matrix = $matrix1-E<gt>multiply($matrix2);>
 2773 
 2774 Calculates the product of matrix "C<$matrix1>" and matrix "C<$matrix2>"
 2775 and returns an object reference to a new matrix "C<$product_matrix>" in
 2776 which the result of this operation has been stored.
 2777 
 2778 Note that the dimensions of the two matrices "C<$matrix1>" and "C<$matrix2>"
 2779 (i.e., their numbers of rows and columns) must harmonize in the following
 2780 way (example):
 2781 
 2782                           [ 2 2 ]
 2783                           [ 2 2 ]
 2784                           [ 2 2 ]
 2785 
 2786               [ 1 1 1 ]   [ * * ]
 2787               [ 1 1 1 ]   [ * * ]
 2788               [ 1 1 1 ]   [ * * ]
 2789               [ 1 1 1 ]   [ * * ]
 2790 
 2791 I.e., the number of columns of matrix "C<$matrix1>" has to be the same
 2792 as the number of rows of matrix "C<$matrix2>".
 2793 
 2794 The number of rows and columns of the resulting matrix "C<$product_matrix>"
 2795 is determined by the number of rows of matrix "C<$matrix1>" and the number
 2796 of columns of matrix "C<$matrix2>", respectively.
 2797 
 2798 =item *
 2799 
 2800 C<$minimum = MatrixReal1::min($number1,$number2);>
 2801 
 2802 Returns the minimum of the two numbers "C<number1>" and "C<number2>".
 2803 
 2804 =item *
 2805 
 2806 C<$minimum = MatrixReal1::max($number1,$number2);>
 2807 
 2808 Returns the maximum of the two numbers "C<number1>" and "C<number2>".
 2809 
 2810 =item *
 2811 
 2812 C<$minimal_cost_matrix = $cost_matrix-E<gt>kleene();>
 2813 
 2814 Copies the matrix "C<$cost_matrix>" (which has to be quadratic!) to
 2815 a new matrix of the same size (i.e., "clones" the input matrix) and
 2816 applies Kleene's algorithm to it.
 2817 
 2818 See L<Math::Kleene(3)> for more details about this algorithm!
 2819 
 2820 The method returns an object reference to the new matrix.
 2821 
 2822 Matrix "C<$cost_matrix>" is not changed by this method in any way.
 2823 
 2824 =item *
 2825 
 2826 C<($norm_matrix,$norm_vector) = $matrix-E<gt>normalize($vector);>
 2827 
 2828 This method is used to improve the numerical stability when solving
 2829 linear equation systems.
 2830 
 2831 Suppose you have a matrix "A" and a vector "b" and you want to find
 2832 out a vector "x" so that C<A * x = b>, i.e., the vector "x" which
 2833 solves the equation system represented by the matrix "A" and the
 2834 vector "b".
 2835 
 2836 Applying this method to the pair (A,b) yields a pair (A',b') where
 2837 each row has been divided by (the absolute value of) the greatest
 2838 coefficient appearing in that row. So this coefficient becomes equal
 2839 to "1" (or "-1") in the new pair (A',b') (all others become smaller
 2840 than one and greater than minus one).
 2841 
 2842 Note that this operation does not change the equation system itself
 2843 because the same division is carried out on either side of the equation
 2844 sign!
 2845 
 2846 The method requires a quadratic (!) matrix "C<$matrix>" and a vector
 2847 "C<$vector>" for input (the vector must be a column vector with the same
 2848 number of rows as the input matrix) and returns a list of two items
 2849 which are object references to a new matrix and a new vector, in this
 2850 order.
 2851 
 2852 The output matrix and vector are clones of the input matrix and vector
 2853 to which the operation explained above has been applied.
 2854 
 2855 The input matrix and vector are not changed by this in any way.
 2856 
 2857 Example of how this method can affect the result of the methods to solve
 2858 equation systems (explained immediately below following this method):
 2859 
 2860 Consider the following little program:
 2861 
 2862   #!perl -w
 2863 
 2864   use MatrixReal1 qw(new_from_string);
 2865 
 2866   $A = MatrixReal1->new_from_string(<<"MATRIX");
 2867   [  1   2   3  ]
 2868   [  5   7  11  ]
 2869   [ 23  19  13  ]
 2870   MATRIX
 2871 
 2872   $b = MatrixReal1->new_from_string(<<"MATRIX");
 2873   [   0   ]
 2874   [   1   ]
 2875   [  29   ]
 2876   MATRIX
 2877 
 2878   $LR = $A->decompose_LR();
 2879   if (($dim,$x,$B) = $LR->solve_LR($b))
 2880   {
 2881       $test = $A * $x;
 2882       print "x = \n$x";
 2883       print "A * x = \n$test";
 2884   }
 2885 
 2886   ($A_,$b_) = $A->normalize($b);
 2887 
 2888   $LR = $A_->decompose_LR();
 2889   if (($dim,$x,$B) = $LR->solve_LR($b_))
 2890   {
 2891       $test = $A * $x;
 2892       print "x = \n$x";
 2893       print "A * x = \n$test";
 2894   }
 2895 
 2896 This will print:
 2897 
 2898   x =
 2899   [  1.000000000000E+00 ]
 2900   [  1.000000000000E+00 ]
 2901   [ -1.000000000000E+00 ]
 2902   A * x =
 2903   [  4.440892098501E-16 ]
 2904   [  1.000000000000E+00 ]
 2905   [  2.900000000000E+01 ]
 2906   x =
 2907   [  1.000000000000E+00 ]
 2908   [  1.000000000000E+00 ]
 2909   [ -1.000000000000E+00 ]
 2910   A * x =
 2911   [  0.000000000000E+00 ]
 2912   [  1.000000000000E+00 ]
 2913   [  2.900000000000E+01 ]
 2914 
 2915 You can see that in the second example (where "normalize()" has been used),
 2916 the result is "better", i.e., more accurate!
 2917 
 2918 =item *
 2919 
 2920 C<$LR_matrix = $matrix-E<gt>decompose_LR();>
 2921 
 2922 This method is needed to solve linear equation systems.
 2923 
 2924 Suppose you have a matrix "A" and a vector "b" and you want to find
 2925 out a vector "x" so that C<A * x = b>, i.e., the vector "x" which
 2926 solves the equation system represented by the matrix "A" and the
 2927 vector "b".
 2928 
 2929 You might also have a matrix "A" and a whole bunch of different
 2930 vectors "b1".."bk" for which you need to find vectors "x1".."xk"
 2931 so that C<A * xi = bi>, for C<i=1..k>.
 2932 
 2933 Using Gaussian transformations (multiplying a row or column with
 2934 a factor, swapping two rows or two columns and adding a multiple
 2935 of one row or column to another), it is possible to decompose any
 2936 matrix "A" into two triangular matrices, called "L" and "R" (for
 2937 "Left" and "Right").
 2938 
 2939 "L" has one's on the main diagonal (the elements (1,1), (2,2), (3,3)
 2940 and so so), non-zero values to the left and below of the main diagonal
 2941 and all zero's in the upper right half of the matrix.
 2942 
 2943 "R" has non-zero values on the main diagonal as well as to the right
 2944 and above of the main diagonal and all zero's in the lower left half
 2945 of the matrix, as follows:
 2946 
 2947           [ 1 0 0 0 0 ]      [ x x x x x ]
 2948           [ x 1 0 0 0 ]      [ 0 x x x x ]
 2949       L = [ x x 1 0 0 ]  R = [ 0 0 x x x ]
 2950           [ x x x 1 0 ]      [ 0 0 0 x x ]
 2951           [ x x x x 1 ]      [ 0 0 0 0 x ]
 2952 
 2953 Note that "C<L * R>" is equivalent to matrix "A" in the sense that
 2954 C<L * R * x = b  E<lt>==E<gt>  A * x = b> for all vectors "x", leaving
 2955 out of account permutations of the rows and columns (these are taken
 2956 care of "magically" by this module!) and numerical errors.
 2957 
 2958 Trick:
 2959 
 2960 Because we know that "L" has one's on its main diagonal, we can
 2961 store both matrices together in the same array without information
 2962 loss! I.e.,
 2963 
 2964                  [ R R R R R ]
 2965                  [ L R R R R ]
 2966             LR = [ L L R R R ]
 2967                  [ L L L R R ]
 2968                  [ L L L L R ]
 2969 
 2970 Beware, though, that "LR" and "C<L * R>" are not the same!!!
 2971 
 2972 Note also that for the same reason, you cannot apply the method "normalize()"
 2973 to an "LR" decomposition matrix. Trying to do so will yield meaningless
 2974 rubbish!
 2975 
 2976 (You need to apply "normalize()" to each pair (Ai,bi) B<BEFORE> decomposing
 2977 the matrix "Ai'"!)
 2978 
 2979 Now what does all this help us in solving linear equation systems?
 2980 
 2981 It helps us because a triangular matrix is the next best thing
 2982 that can happen to us besides a diagonal matrix (a matrix that
 2983 has non-zero values only on its main diagonal - in which case
 2984 the solution is trivial, simply divide "C<b[i]>" by "C<A[i,i]>"
 2985 to get "C<x[i]>"!).
 2986 
 2987 To find the solution to our problem "C<A * x = b>", we divide this
 2988 problem in parts: instead of solving C<A * x = b> directly, we first
 2989 decompose "A" into "L" and "R" and then solve "C<L * y = b>" and
 2990 finally "C<R * x = y>" (motto: divide and rule!).
 2991 
 2992 From the illustration above it is clear that solving "C<L * y = b>"
 2993 and "C<R * x = y>" is straightforward: we immediately know that
 2994 C<y[1] = b[1]>. We then deduce swiftly that
 2995 
 2996   y[2] = b[2] - L[2,1] * y[1]
 2997 
 2998 (and we know "C<y[1]>" by now!), that
 2999 
 3000   y[3] = b[3] - L[3,1] * y[1] - L[3,2] * y[2]
 3001 
 3002 and so on.
 3003 
 3004 Having effortlessly calculated the vector "y", we now proceed to
 3005 calculate the vector "x" in a similar fashion: we see immediately
 3006 that C<x[n] = y[n] / R[n,n]>. It follows that
 3007 
 3008   x[n-1] = ( y[n-1] - R[n-1,n] * x[n] ) / R[n-1,n-1]
 3009 
 3010 and
 3011 
 3012   x[n-2] = ( y[n-2] - R[n-2,n-1] * x[n-1] - R[n-2,n] * x[n] )
 3013            / R[n-2,n-2]
 3014 
 3015 and so on.
 3016 
 3017 You can see that - especially when you have many vectors "b1".."bk"
 3018 for which you are searching solutions to C<A * xi = bi> - this scheme
 3019 is much more efficient than a straightforward, "brute force" approach.
 3020 
 3021 This method requires a quadratic matrix as its input matrix.
 3022 
 3023 If you don't have that many equations, fill up with zero's (i.e., do
 3024 nothing to fill the superfluous rows if it's a "fresh" matrix, i.e.,
 3025 a matrix that has been created with "new()" or "shadow()").
 3026 
 3027 The method returns an object reference to a new matrix containing the
 3028 matrices "L" and "R".
 3029 
 3030 The input matrix is not changed by this method in any way.
 3031 
 3032 Note that you can "copy()" or "clone()" the result of this method without
 3033 losing its "magical" properties (for instance concerning the hidden
 3034 permutations of its rows and columns).
 3035 
 3036 However, as soon as you are applying any method that alters the contents
 3037 of the matrix, its "magical" properties are stripped off, and the matrix
 3038 immediately reverts to an "ordinary" matrix (with the values it just happens
 3039 to contain at that moment, be they meaningful as an ordinary matrix or not!).
 3040 
 3041 =item *
 3042 
 3043 C<($dimension,$x_vector,$base_matrix) = $LR_matrix>C<-E<gt>>C<solve_LR($b_vector);>
 3044 
 3045 Use this method to actually solve an equation system.
 3046 
 3047 Matrix "C<$LR_matrix>" must be a (quadratic) matrix returned by the
 3048 method "decompose_LR()", the LR decomposition matrix of the matrix
 3049 "A" of your equation system C<A * x = b>.
 3050 
 3051 The input vector "C<$b_vector>" is the vector "b" in your equation system
 3052 C<A * x = b>, which must be a column vector and have the same number of
 3053 rows as the input matrix "C<$LR_matrix>".
 3054 
 3055 The method returns a list of three items if a solution exists or an
 3056 empty list otherwise (!).
 3057 
 3058 Therefore, you should always use this method like this:
 3059 
 3060   if ( ($dim,$x_vec,$base) = $LR->solve_LR($b_vec) )
 3061   {
 3062       # do something with the solution...
 3063   }
 3064   else
 3065   {
 3066       # do something with the fact that there is no solution...
 3067   }
 3068 
 3069 The three items returned are: the dimension "C<$dimension>" of the solution
 3070 space (which is zero if only one solution exists, one if the solution is
 3071 a straight line, two if the solution is a plane, and so on), the solution
 3072 vector "C<$x_vector>" (which is the vector "x" of your equation system
 3073 C<A * x = b>) and a matrix "C<$base_matrix>" representing a base of the
 3074 solution space (a set of vectors which put up the solution space like
 3075 the spokes of an umbrella).
 3076 
 3077 Only the first "C<$dimension>" columns of this base matrix actually
 3078 contain entries, the remaining columns are all zero.
 3079 
 3080 Now what is all this stuff with that "base" good for?
 3081 
 3082 The output vector "x" is B<ALWAYS> a solution of your equation system
 3083 C<A * x = b>.
 3084 
 3085 But also any vector "C<$vector>"
 3086 
 3087   $vector = $x_vector->clone();
 3088 
 3089   $machine_infinity = 1E+99; # or something like that
 3090 
 3091   for ( $i = 1; $i <= $dimension; $i++ )
 3092   {
 3093       $vector += rand($machine_infinity) * $base_matrix->column($i);
 3094   }
 3095 
 3096 is a solution to your problem C<A * x = b>, i.e., if "C<$A_matrix>" contains
 3097 your matrix "A", then
 3098 
 3099   print abs( $A_matrix * $vector - $b_vector ), "\n";
 3100 
 3101 should print a number around 1E-16 or so!
 3102 
 3103 By the way, note that you can actually calculate those vectors "C<$vector>"
 3104 a little more efficient as follows:
 3105 
 3106   $rand_vector = $x_vector->shadow();
 3107 
 3108   $machine_infinity = 1E+99; # or something like that
 3109 
 3110   for ( $i = 1; $i <= $dimension; $i++ )
 3111   {
 3112       $rand_vector->assign($i,1, rand($machine_infinity) );
 3113   }
 3114 
 3115   $vector = $x_vector + ( $base_matrix * $rand_vector );
 3116 
 3117 Note that the input matrix and vector are not changed by this method
 3118 in any way.
 3119 
 3120 =item *
 3121 
 3122 C<$inverse_matrix = $LR_matrix-E<gt>invert_LR();>
 3123 
 3124 Use this method to calculate the inverse of a given matrix "C<$LR_matrix>",
 3125 which must be a (quadratic) matrix returned by the method "decompose_LR()".
 3126 
 3127 The method returns an object reference to a new matrix of the same size as
 3128 the input matrix containing the inverse of the matrix that you initially
 3129 fed into "decompose_LR()" B<IF THE INVERSE EXISTS>, or an empty list
 3130 otherwise.
 3131 
 3132 Therefore, you should always use this method in the following way:
 3133 
 3134   if ( $inverse_matrix = $LR->invert_LR() )
 3135   {
 3136       # do something with the inverse matrix...
 3137   }
 3138   else
 3139   {
 3140       # do something with the fact that there is no inverse matrix...
 3141   }
 3142 
 3143 Note that by definition (disregarding numerical errors), the product
 3144 of the initial matrix and its inverse (or vice-versa) is always a matrix
 3145 containing one's on the main diagonal (elements (1,1), (2,2), (3,3) and
 3146 so on) and zero's elsewhere.
 3147 
 3148 The input matrix is not changed by this method in any way.
 3149 
 3150 =item *
 3151 
 3152 C<$condition = $matrix-E<gt>condition($inverse_matrix);>
 3153 
 3154 In fact this method is just a shortcut for
 3155 
 3156   abs($matrix) * abs($inverse_matrix)
 3157 
 3158 Both input matrices must be quadratic and have the same size, and the result
 3159 is meaningful only if one of them is the inverse of the other (for instance,
 3160 as returned by the method "invert_LR()").
 3161 
 3162 The number returned is a measure of the "condition" of the given matrix
 3163 "C<$matrix>", i.e., a measure of the numerical stability of the matrix.
 3164 
 3165 This number is always positive, and the smaller its value, the better the
 3166 condition of the matrix (the better the stability of all subsequent
 3167 computations carried out using this matrix).
 3168 
 3169 Numerical stability means for example that if
 3170 
 3171   abs( $vec_correct - $vec_with_error ) < $epsilon
 3172 
 3173 holds, there must be a "C<$delta>" which doesn't depend on the vector
 3174 "C<$vec_correct>" (nor "C<$vec_with_error>", by the way) so that
 3175 
 3176   abs( $matrix * $vec_correct - $matrix * $vec_with_error ) < $delta
 3177 
 3178 also holds.
 3179 
 3180 =item *
 3181 
 3182 C<$determinant = $LR_matrix-E<gt>det_LR();>
 3183 
 3184 Calculates the determinant of a matrix, whose LR decomposition matrix
 3185 "C<$LR_matrix>" must be given (which must be a (quadratic) matrix
 3186 returned by the method "decompose_LR()").
 3187 
 3188 In fact the determinant is a by-product of the LR decomposition: It is
 3189 (in principle, that is, except for the sign) simply the product of the
 3190 elements on the main diagonal (elements (1,1), (2,2), (3,3) and so on)
 3191 of the LR decomposition matrix.
 3192 
 3193 (The sign is taken care of "magically" by this module)
 3194 
 3195 =item *
 3196 
 3197 C<$order = $LR_matrix-E<gt>order_LR();>
 3198 
 3199 Calculates the order (called "Rang" in German) of a matrix, whose
 3200 LR decomposition matrix "C<$LR_matrix>" must be given (which must
 3201 be a (quadratic) matrix returned by the method "decompose_LR()").
 3202 
 3203 This number is a measure of the number of linear independent row
 3204 and column vectors (= number of linear independent equations in
 3205 the case of a matrix representing an equation system) of the
 3206 matrix that was initially fed into "decompose_LR()".
 3207 
 3208 If "n" is the number of rows and columns of the (quadratic!) matrix,
 3209 then "n - order" is the dimension of the solution space of the
 3210 associated equation system.
 3211 
 3212 =item *
 3213 
 3214 C<$scalar_product = $vector1-E<gt>scalar_product($vector2);>
 3215 
 3216 Returns the scalar product of vector "C<$vector1>" and vector "C<$vector2>".
 3217 
 3218 Both vectors must be column vectors (i.e., a matrix having
 3219 several rows but only one column).
 3220 
 3221 This is a (more efficient!) shortcut for
 3222 
 3223   $temp           = ~$vector1 * $vector2;
 3224   $scalar_product =  $temp->element(1,1);
 3225 
 3226 or the sum C<i=1..n> of the products C<vector1[i] * vector2[i]>.
 3227 
 3228 Provided none of the two input vectors is the null vector, then
 3229 the two vectors are orthogonal, i.e., have an angle of 90 degrees
 3230 between them, exactly when their scalar product is zero, and
 3231 vice-versa.
 3232 
 3233 =item *
 3234 
 3235 C<$vector_product = $vector1-E<gt>vector_product($vector2);>
 3236 
 3237 Returns the vector product of vector "C<$vector1>" and vector "C<$vector2>".
 3238 
 3239 Both vectors must be column vectors (i.e., a matrix having several rows
 3240 but only one column).
 3241 
 3242 Currently, the vector product is only defined for 3 dimensions (i.e.,
 3243 vectors with 3 rows); all other vectors trigger an error message.
 3244 
 3245 In 3 dimensions, the vector product of two vectors "x" and "y"
 3246 is defined as
 3247 
 3248               |  x[1]  y[1]  e[1]  |
 3249   determinant |  x[2]  y[2]  e[2]  |
 3250               |  x[3]  y[3]  e[3]  |
 3251 
 3252 where the "C<x[i]>" and "C<y[i]>" are the components of the two vectors
 3253 "x" and "y", respectively, and the "C<e[i]>" are unity vectors (i.e.,
 3254 vectors with a length equal to one) with a one in row "i" and zero's
 3255 elsewhere (this means that you have numbers and vectors as elements
 3256 in this matrix!).
 3257 
 3258 This determinant evaluates to the rather simple formula
 3259 
 3260   z[1] = x[2] * y[3] - x[3] * y[2]
 3261   z[2] = x[3] * y[1] - x[1] * y[3]
 3262   z[3] = x[1] * y[2] - x[2] * y[1]
 3263 
 3264 A characteristic property of the vector product is that the resulting
 3265 vector is orthogonal to both of the input vectors (if neither of both
 3266 is the null vector, otherwise this is trivial), i.e., the scalar product
 3267 of each of the input vectors with the resulting vector is always zero.
 3268 
 3269 =item *
 3270 
 3271 C<$length = $vector-E<gt>length();>
 3272 
 3273 This is actually a shortcut for
 3274 
 3275   $length = sqrt( $vector->scalar_product($vector) );
 3276 
 3277 and returns the length of a given (column!) vector "C<$vector>".
 3278 
 3279 Note that the "length" calculated by this method is in fact the
 3280 "two"-norm of a vector "C<$vector>"!
 3281 
 3282 The general definition for norms of vectors is the following:
 3283 
 3284   sub vector_norm
 3285   {
 3286       croak "Usage: \$norm = \$vector->vector_norm(\$n);"
 3287         if (@_ != 2);
 3288 
 3289       my($vector,$n) = @_;
 3290       my($rows,$cols) = ($vector->[1],$vector->[2]);
 3291       my($k,$comp,$sum);
 3292 
 3293       croak "MatrixReal1::vector_norm(): vector is not a column vector"
 3294         unless ($cols == 1);
 3295 
 3296       croak "MatrixReal1::vector_norm(): norm index must be > 0"
 3297         unless ($n > 0);
 3298 
 3299       croak "MatrixReal1::vector_norm(): norm index must be integer"
 3300         unless ($n == int($n));
 3301 
 3302       $sum = 0;
 3303       for ( $k = 0; $k < $rows; $k++ )
 3304       {
 3305           $comp = abs( $vector->[0][$k][0] );
 3306           $sum += $comp ** $n;
 3307       }
 3308       return( $sum ** (1 / $n) );
 3309   }
 3310 
 3311 Note that the case "n = 1" is the "one"-norm for matrices applied to a
 3312 vector, the case "n = 2" is the euclidian norm or length of a vector,
 3313 and if "n" goes to infinity, you have the "infinity"- or "maximum"-norm
 3314 for matrices applied to a vector!
 3315 
 3316 =item *
 3317 
 3318 C<$xn_vector = $matrix-E<gt>>C<solve_GSM($x0_vector,$b_vector,$epsilon);>
 3319 
 3320 =item *
 3321 
 3322 C<$xn_vector = $matrix-E<gt>>C<solve_SSM($x0_vector,$b_vector,$epsilon);>
 3323 
 3324 =item *
 3325 
 3326 C<$xn_vector = $matrix-E<gt>>C<solve_RM($x0_vector,$b_vector,$weight,$epsilon);>
 3327 
 3328 In some cases it might not be practical or desirable to solve an
 3329 equation system "C<A * x = b>" using an analytical algorithm like
 3330 the "decompose_LR()" and "solve_LR()" method pair.
 3331 
 3332 In fact in some cases, due to the numerical properties (the "condition")
 3333 of the matrix "A", the numerical error of the obtained result can be
 3334 greater than by using an approximative (iterative) algorithm like one
 3335 of the three implemented here.
 3336 
 3337 All three methods, GSM ("Global Step Method" or "Gesamtschrittverfahren"),
 3338 SSM ("Single Step Method" or "Einzelschrittverfahren") and RM ("Relaxation
 3339 Method" or "Relaxationsverfahren"), are fix-point iterations, that is, can
 3340 be described by an iteration function "C<x(t+1) = Phi( x(t) )>" which has
 3341 the property:
 3342 
 3343   Phi(x)  =  x    <==>    A * x  =  b
 3344 
 3345 We can define "C<Phi(x)>" as follows:
 3346 
 3347   Phi(x)  :=  ( En - A ) * x  +  b
 3348 
 3349 where "En" is a matrix of the same size as "A" ("n" rows and columns)
 3350 with one's on its main diagonal and zero's elsewhere.
 3351 
 3352 This function has the required property.
 3353 
 3354 Proof:
 3355 
 3356            A * x        =   b
 3357 
 3358   <==>  -( A * x )      =  -b
 3359 
 3360   <==>  -( A * x ) + x  =  -b + x
 3361 
 3362   <==>  -( A * x ) + x + b  =  x
 3363 
 3364   <==>  x - ( A * x ) + b  =  x
 3365 
 3366   <==>  ( En - A ) * x + b  =  x
 3367 
 3368 This last step is true because
 3369 
 3370   x[i] - ( a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] ) + b[i]
 3371 
 3372 is the same as
 3373 
 3374   ( -a[i,1] x[1] + ... + (1 - a[i,i]) x[i] + ... + -a[i,n] x[n] ) + b[i]
 3375 
 3376 qed
 3377 
 3378 Note that actually solving the equation system "C<A * x = b>" means
 3379 to calculate
 3380 
 3381         a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n]  =  b[i]
 3382 
 3383   <==>  a[i,i] x[i]  =
 3384         b[i]
 3385         - ( a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] )
 3386         + a[i,i] x[i]
 3387 
 3388   <==>  x[i]  =
 3389         ( b[i]
 3390             - ( a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] )
 3391             + a[i,i] x[i]
 3392         ) / a[i,i]
 3393 
 3394   <==>  x[i]  =
 3395         ( b[i] -
 3396             ( a[i,1] x[1] + ... + a[i,i-1] x[i-1] +
 3397               a[i,i+1] x[i+1] + ... + a[i,n] x[n] )
 3398         ) / a[i,i]
 3399 
 3400 There is one major restriction, though: a fix-point iteration is
 3401 guaranteed to converge only if the first derivative of the iteration
 3402 function has an absolute value less than one in an area around the
 3403 point "C<x(*)>" for which "C<Phi( x(*) ) = x(*)>" is to be true, and
 3404 if the start vector "C<x(0)>" lies within that area!
 3405 
 3406 This is best verified grafically, which unfortunately is impossible
 3407 to do in this textual documentation!
 3408 
 3409 See literature on Numerical Analysis for details!
 3410 
 3411 In our case, this restriction translates to the following three conditions:
 3412 
 3413 There must exist a norm so that the norm of the matrix of the iteration
 3414 function, C<( En - A )>, has a value less than one, the matrix "A" may
 3415 not have any zero value on its main diagonal and the initial vector
 3416 "C<x(0)>" must be "good enough", i.e., "close enough" to the solution
 3417 "C<x(*)>".
 3418 
 3419 (Remember school math: the first derivative of a straight line given by
 3420 "C<y = a * x + b>" is "a"!)
 3421 
 3422 The three methods expect a (quadratic!) matrix "C<$matrix>" as their
 3423 first argument, a start vector "C<$x0_vector>", a vector "C<$b_vector>"
 3424 (which is the vector "b" in your equation system "C<A * x = b>"), in the
 3425 case of the "Relaxation Method" ("RM"), a real number "C<$weight>" best
 3426 between zero and two, and finally an error limit (real number) "C<$epsilon>".
 3427 
 3428 (Note that the weight "C<$weight>" used by the "Relaxation Method" ("RM")
 3429 is B<NOT> checked to lie within any reasonable range!)
 3430 
 3431 The three methods first test the first two conditions of the three
 3432 conditions listed above and return an empty list if these conditions
 3433 are not fulfilled.
 3434 
 3435 Therefore, you should always test their return value using some
 3436 code like:
 3437 
 3438   if ( $xn_vector = $A_matrix->solve_GSM($x0_vector,$b_vector,1E-12) )
 3439   {
 3440       # do something with the solution...
 3441   }
 3442   else
 3443   {
 3444       # do something with the fact that there is no solution...
 3445   }
 3446 
 3447 Otherwise, they iterate until C<abs( Phi(x) - x ) E<lt> epsilon>.
 3448 
 3449 (Beware that theoretically, infinite loops might result if the starting
 3450 vector is too far "off" the solution! In practice, this shouldn't be
 3451 a problem. Anyway, you can always press <ctrl-C> if you think that the
 3452 iteration takes too long!)
 3453 
 3454 The difference between the three methods is the following:
 3455 
 3456 In the "Global Step Method" ("GSM"), the new vector "C<x(t+1)>"
 3457 (called "y" here) is calculated from the vector "C<x(t)>"
 3458 (called "x" here) according to the formula:
 3459 
 3460   y[i] =
 3461   ( b[i]
 3462       - ( a[i,1] x[1] + ... + a[i,i-1] x[i-1] +
 3463           a[i,i+1] x[i+1] + ... + a[i,n] x[n] )
 3464   ) / a[i,i]
 3465 
 3466 In the "Single Step Method" ("SSM"), the components of the vector
 3467 "C<x(t+1)>" which have already been calculated are used to calculate
 3468 the remaining components, i.e.
 3469 
 3470   y[i] =
 3471   ( b[i]
 3472       - ( a[i,1] y[1] + ... + a[i,i-1] y[i-1] +  # note the "y[]"!
 3473           a[i,i+1] x[i+1] + ... + a[i,n] x[n] )  # note the "x[]"!
 3474   ) / a[i,i]
 3475 
 3476 In the "Relaxation method" ("RM"), the components of the vector
 3477 "C<x(t+1)>" are calculated by "mixing" old and new value (like
 3478 cold and hot water), and the weight "C<$weight>" determines the
 3479 "aperture" of both the "hot water tap" as well as of the "cold
 3480 water tap", according to the formula:
 3481 
 3482   y[i] =
 3483   ( b[i]
 3484       - ( a[i,1] y[1] + ... + a[i,i-1] y[i-1] +  # note the "y[]"!
 3485           a[i,i+1] x[i+1] + ... + a[i,n] x[n] )  # note the "x[]"!
 3486   ) / a[i,i]
 3487   y[i] = weight * y[i] + (1 - weight) * x[i]
 3488 
 3489 Note that the weight "C<$weight>" should be greater than zero and
 3490 less than two (!).
 3491 
 3492 The three methods are supposed to be of different efficiency.
 3493 Experiment!
 3494 
 3495 Remember that in most cases, it is probably advantageous to first
 3496 "normalize()" your equation system prior to solving it!
 3497 
 3498 =back
 3499 
 3500 =head2 Eigensystems
 3501 
 3502 =over 2
 3503 
 3504 =item *
 3505 
 3506 C<$matrix-E<gt>is_symmetric();>
 3507 
 3508 Returns a boolean value indicating if the given matrix is
 3509 symmetric (B<M>[I<i>,I<j>]=B<M>[I<j>,I<i>]). This is equivalent to
 3510 C<($matrix == ~$matrix)> but without memory allocation.
 3511 
 3512 =item *
 3513 
 3514 C<($l, $V) = $matrix-E<gt>sym_diagonalize();>
 3515 
 3516 This method performs the diagonalization of the quadratic
 3517 I<symmetric> matrix B<M> stored in $matrix.
 3518 On output, B<l> is a column vector containing all the eigenvalues
 3519 of B<M> and B<V> is an orthogonal matrix which columns are the
 3520 corresponding normalized eigenvectors.
 3521 The primary property of an eigenvalue I<l> and an eigenvector
 3522 B<x> is of course that: B<M> * B<x> = I<l> * B<x>.
 3523 
 3524 The method uses a Householder reduction to tridiagonal form
 3525 followed by a QL algoritm with implicit shifts on this
 3526 tridiagonal. (The tridiagonal matrix is kept internally
 3527 in a compact form in this routine to save memory.)
 3528 In fact, this routine wraps the householder() and
 3529 tri_diagonalize() methods described below when their
 3530 intermediate results are not desired.
 3531 The overall algorithmic complexity of this technique
 3532 is O(N^3). According to several books, the coefficient
 3533 hidden by the 'O' is one of the best possible for general
 3534 (symmetric) matrixes.
 3535 
 3536 =item *
 3537 
 3538 C<($T, $Q) = $matrix-E<gt>householder();>
 3539 
 3540 This method performs the Householder algorithm which reduces
 3541 the I<n> by I<n> real I<symmetric> matrix B<M> contained
 3542 in $matrix to tridiagonal form.
 3543 On output, B<T> is a symmetric tridiagonal matrix (only
 3544 diagonal and off-diagonal elements are non-zero) and B<Q>
 3545 is an I<orthogonal> matrix performing the tranformation
 3546 between B<M> and B<T> (C<$M == $Q * $T * ~$Q>).
 3547 
 3548 =item *
 3549 
 3550 C<($l, $V) = $T-E<gt>tri_diagonalize([$Q]);>
 3551 
 3552 This method diagonalizes the symmetric tridiagonal
 3553 matrix B<T>. On output, $l and $V are similar to the
 3554 output values described for sym_diagonalize().
 3555 
 3556 The optional argument $Q corresponds to an orthogonal
 3557 transformation matrix B<Q> that should be used additionally
 3558 during B<V> (eigenvectors) computation. It should be supplied
 3559 if the desired eigenvectors correspond to a more general
 3560 symmetric matrix B<M> previously reduced by the
 3561 householder() method, not a mere tridiagonal. If B<T> is
 3562 really a tridiagonal matrix, B<Q> can be omitted (it
 3563 will be internally created in fact as an identity matrix).
 3564 The method uses a QL algorithm (with implicit shifts).
 3565 
 3566 =item *
 3567 
 3568 C<$l = $matrix-E<gt>sym_eigenvalues();>
 3569 
 3570 This method computes the eigenvalues of the quadratic
 3571 I<symmetric> matrix B<M> stored in $matrix.
 3572 On output, B<l> is a column vector containing all the eigenvalues
 3573 of B<M>. Eigenvectors are not computed (on the contrary of
 3574 C<sym_diagonalize()>) and this method is more efficient
 3575 (even though it uses a similar algorithm with two phases).
 3576 However, understand that the algorithmic complexity of this
 3577 technique is still also O(N^3). But the coefficient hidden
 3578 by the 'O' is better by a factor of..., well, see your
 3579 benchmark, it's wiser.
 3580 
 3581 This routine wraps the householder_tridiagonal() and
 3582 tri_eigenvalues() methods described below when the
 3583 intermediate tridiagonal matrix is not needed.
 3584 
 3585 =item *
 3586 
 3587 C<$T = $matrix-E<gt>householder_tridiagonal();>
 3588 
 3589 This method performs the Householder algorithm which reduces
 3590 the I<n> by I<n> real I<symmetric> matrix B<M> contained
 3591 in $matrix to tridiagonal form.
 3592 On output, B<T> is the obtained symmetric tridiagonal matrix
 3593 (only diagonal and off-diagonal elements are non-zero). The
 3594 operation is similar to the householder() method, but potentially
 3595 a little more efficient as the transformation matrix is not
 3596 computed.
 3597 
 3598 =item *
 3599 
 3600 C<$l = $T-E<gt>tri_eigenvalues();>
 3601 
 3602 This method compute the eigenvalues of the symmetric
 3603 tridiagonal matrix B<T>. On output, $l is a vector
 3604 containing the eigenvalues (similar to C<sym_eigenvalues()>).
 3605 This method is much more efficient than tri_diagonalize()
 3606 when eigenvectors are not needed.
 3607 
 3608 =back
 3609 
 3610 =head1 OVERLOADED OPERATORS
 3611 
 3612 =head2 SYNOPSIS
 3613 
 3614 =over 2
 3615 
 3616 =item *
 3617 
 3618 Unary operators:
 3619 
 3620 "C<->", "C<~>", "C<abs>", C<test>, "C<!>", 'C<"">'
 3621 
 3622 =item *
 3623 
 3624 Binary (arithmetic) operators:
 3625 
 3626 "C<+>", "C<->", "C<*>"
 3627 
 3628 =item *
 3629 
 3630 Binary (relational) operators:
 3631 
 3632 "C<==>", "C<!=>", "C<E<lt>>", "C<E<lt>=>", "C<E<gt>>", "C<E<gt>=>"
 3633 
 3634 "C<eq>", "C<ne>", "C<lt>", "C<le>", "C<gt>", "C<ge>"
 3635 
 3636 Note that the latter ("C<eq>", "C<ne>", ... ) are just synonyms
 3637 of the former ("C<==>", "C<!=>", ... ), defined for convenience
 3638 only.
 3639 
 3640 =back
 3641 
 3642 =head2 DESCRIPTION
 3643 
 3644 =over 5
 3645 
 3646 =item '-'
 3647 
 3648 Unary minus
 3649 
 3650 Returns the negative of the given matrix, i.e., the matrix with
 3651 all elements multiplied with the factor "-1".
 3652 
 3653 Example:
 3654 
 3655     $matrix = -$matrix;
 3656 
 3657 =item '~'
 3658 
 3659 Transposition
 3660 
 3661 Returns the transposed of the given matrix.
 3662 
 3663 Examples:
 3664 
 3665     $temp = ~$vector * $vector;
 3666     $length = sqrt( $temp->element(1,1) );
 3667 
 3668     if (~$matrix == $matrix) { # matrix is symmetric ... }
 3669 
 3670 =item abs
 3671 
 3672 Norm
 3673 
 3674 Returns the "one"-Norm of the given matrix.
 3675 
 3676 Example:
 3677 
 3678     $error = abs( $A * $x - $b );
 3679 
 3680 =item test
 3681 
 3682 Boolean test
 3683 
 3684 Tests wether there is at least one non-zero element in the matrix.
 3685 
 3686 Example:
 3687 
 3688     if ($xn_vector) { # result of iteration is not zero ... }
 3689 
 3690 =item '!'
 3691 
 3692 Negated boolean test
 3693 
 3694 Tests wether the matrix contains only zero's.
 3695 
 3696 Examples:
 3697 
 3698     if (! $b_vector) { # heterogenous equation system ... }
 3699     else             { # homogenous equation system ... }
 3700 
 3701     unless ($x_vector) { # not the null-vector! }
 3702 
 3703 =item '""""'
 3704 
 3705 "Stringify" operator
 3706 
 3707 Converts the given matrix into a string.
 3708 
 3709 Uses scientific representation to keep precision loss to a minimum in case
 3710 you want to read this string back in again later with "new_from_string()".
 3711 
 3712 Uses a 13-digit mantissa and a 20-character field for each element so that
 3713 lines will wrap nicely on an 80-column screen.
 3714 
 3715 Examples:
 3716 
 3717     $matrix = MatrixReal1->new_from_string(<<"MATRIX");
 3718     [ 1  0 ]
 3719     [ 0 -1 ]
 3720     MATRIX
 3721     print "$matrix";
 3722 
 3723     [  1.000000000000E+00  0.000000000000E+00 ]
 3724     [  0.000000000000E+00 -1.000000000000E+00 ]
 3725 
 3726     $string = "$matrix";
 3727     $test = MatrixReal1->new_from_string($string);
 3728     if ($test == $matrix) { print ":-)\n"; } else { print ":-(\n"; }
 3729 
 3730 =item '+'
 3731 
 3732 Addition
 3733 
 3734 Returns the sum of the two given matrices.
 3735 
 3736 Examples:
 3737 
 3738     $matrix_S = $matrix_A + $matrix_B;
 3739 
 3740     $matrix_A += $matrix_B;
 3741 
 3742 =item '-'
 3743 
 3744 Subtraction
 3745 
 3746 Returns the difference of the two given matrices.
 3747 
 3748 Examples:
 3749 
 3750     $matrix_D = $matrix_A - $matrix_B;
 3751 
 3752     $matrix_A -= $matrix_B;
 3753 
 3754 Note that this is the same as:
 3755 
 3756     $matrix_S = $matrix_A + -$matrix_B;
 3757 
 3758     $matrix_A += -$matrix_B;
 3759 
 3760 (The latter are less efficient, though)
 3761 
 3762 =item '*'
 3763 
 3764 Multiplication
 3765 
 3766 Returns the matrix product of the two given matrices or
 3767 the product of the given matrix and scalar factor.
 3768 
 3769 Examples:
 3770 
 3771     $matrix_P = $matrix_A * $matrix_B;
 3772 
 3773     $matrix_A *= $matrix_B;
 3774 
 3775     $vector_b = $matrix_A * $vector_x;
 3776 
 3777     $matrix_B = -1 * $matrix_A;
 3778 
 3779     $matrix_B = $matrix_A * -1;
 3780 
 3781     $matrix_A *= -1;
 3782 
 3783 =item '=='
 3784 
 3785 Equality
 3786 
 3787 Tests two matrices for equality.
 3788 
 3789 Example:
 3790 
 3791     if ( $A * $x == $b ) { print "EUREKA!\n"; }
 3792 
 3793 Note that in most cases, due to numerical errors (due to the finite
 3794 precision of computer arithmetics), it is a bad idea to compare two
 3795 matrices or vectors this way.
 3796 
 3797 Better use the norm of the difference of the two matrices you want
 3798 to compare and compare that norm with a small number, like this:
 3799 
 3800     if ( abs( $A * $x - $b ) < 1E-12 ) { print "BINGO!\n"; }
 3801 
 3802 =item '!='
 3803 
 3804 Inequality
 3805 
 3806 Tests two matrices for inequality.
 3807 
 3808 Example:
 3809 
 3810     while ($x0_vector != $xn_vector) { # proceed with iteration ... }
 3811 
 3812 (Stops when the iteration becomes stationary)
 3813 
 3814 Note that (just like with the '==' operator), it is usually a bad idea
 3815 to compare matrices or vectors this way. Compare the norm of the difference
 3816 of the two matrices with a small number instead.
 3817 
 3818 =item 'E<lt>'
 3819 
 3820 Less than
 3821 
 3822 Examples:
 3823 
 3824     if ( $matrix1 < $matrix2 ) { # ... }
 3825 
 3826     if ( $vector < $epsilon ) { # ... }
 3827 
 3828     if ( 1E-12 < $vector ) { # ... }
 3829 
 3830     if ( $A * $x - $b < 1E-12 ) { # ... }
 3831 
 3832 These are just shortcuts for saying:
 3833 
 3834     if ( abs($matrix1) < abs($matrix2) ) { # ... }
 3835 
 3836     if ( abs($vector) < abs($epsilon) ) { # ... }
 3837 
 3838     if ( abs(1E-12) < abs($vector) ) { # ... }
 3839 
 3840     if ( abs( $A * $x - $b ) < abs(1E-12) ) { # ... }
 3841 
 3842 Uses the "one"-norm for matrices and Perl's built-in "abs()" for scalars.
 3843 
 3844 =item 'E<lt>='
 3845 
 3846 Less than or equal
 3847 
 3848 As with the '<' operator, this is just a shortcut for the same expression
 3849 with "abs()" around all arguments.
 3850 
 3851 Example:
 3852 
 3853     if ( $A * $x - $b <= 1E-12 ) { # ... }
 3854 
 3855 which in fact is the same as:
 3856 
 3857     if ( abs( $A * $x - $b ) <= abs(1E-12) ) { # ... }
 3858 
 3859 Uses the "one"-norm for matrices and Perl's built-in "abs()" for scalars.
 3860 
 3861 =item 'E<gt>'
 3862 
 3863 Greater than
 3864 
 3865 As with the '<' and '<=' operator, this
 3866 
 3867     if ( $xn - $x0 > 1E-12 ) { # ... }
 3868 
 3869 is just a shortcut for:
 3870 
 3871     if ( abs( $xn - $x0 ) > abs(1E-12) ) { # ... }
 3872 
 3873 Uses the "one"-norm for matrices and Perl's built-in "abs()" for scalars.
 3874 
 3875 =item 'E<gt>='
 3876 
 3877 Greater than or equal
 3878 
 3879 As with the '<', '<=' and '>' operator, the following
 3880 
 3881     if ( $LR >= $A ) { # ... }
 3882 
 3883 is simply a shortcut for:
 3884 
 3885     if ( abs($LR) >= abs($A) ) { # ... }
 3886 
 3887 Uses the "one"-norm for matrices and Perl's built-in "abs()" for scalars.
 3888 
 3889 =back
 3890 
 3891 =head1 SEE ALSO
 3892 
 3893 Math::MatrixBool(3), DFA::Kleene(3), Math::Kleene(3),
 3894 Set::IntegerRange(3), Set::IntegerFast(3).
 3895 
 3896 =head1 VERSION
 3897 
 3898 This man page documents MatrixReal1 version 1.3.
 3899 
 3900 =head1 AUTHORS
 3901 
 3902 Steffen Beyer <sb@sdm.de>, Rodolphe Ortalo <ortalo@laas.fr>.
 3903 
 3904 =head1 CREDITS
 3905 
 3906 Many thanks to Prof. Pahlings for stoking the fire of my enthusiasm for
 3907 Algebra and Linear Algebra at the university (RWTH Aachen, Germany), and
 3908 to Prof. Esser and his assistant, Mr. Jarausch, for their fascinating
 3909 lectures in Numerical Analysis!
 3910 
 3911 =head1 COPYRIGHT
 3912 
 3913 Copyright (c) 1996, 1997, 1999 by Steffen Beyer and Rodolphe Ortalo.
 3914 All rights reserved.
 3915 
 3916 =head1 LICENSE AGREEMENT
 3917 
 3918 This package is free software; you can redistribute it and/or
 3919 modify it under the same terms as Perl itself.
 3920 

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9