[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 3372 - (download) (as text) (annotate)
Wed Jul 13 01:22:43 2005 UTC (14 years, 7 months ago) by gage
File size: 106890 byte(s)
Made modifications that allow the use of complex numbers in matrices

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

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9