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