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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3360 - (download) (as text) (annotate)
Tue Jul 5 22:13:16 2005 UTC (14 years, 7 months ago) by gage
File size: 11207 byte(s)
Added some pod documentation

    1 
    2 BEGIN {
    3   be_strict(); # an alias for use strict.  This means that all global variable must contain main:: as a prefix.
    4 
    5 }
    6 package Matrix;
    7 @Matrix::ISA = qw(MatrixReal1);
    8 
    9 use Carp;
   10 
   11 $Matrix::DEFAULT_FORMAT = '% #-19.12E ';
   12 # allows specification of the format
   13 
   14 =head4
   15 
   16   Method $matrix->_stringify()
   17   -- overrides MatrixReal1 display mode
   18 
   19 =cut
   20 
   21 
   22 sub _stringify
   23 {
   24     my($object,$argument,$flag) = @_;
   25 #   my($name) = '""'; #&_trace($name,$object,$argument,$flag);
   26     my($rows,$cols) = ($object->[1],$object->[2]);
   27     my($i,$j,$s);
   28 
   29     $s = '';
   30     for ( $i = 0; $i < $rows; $i++ )
   31     {
   32         $s .= "[ ";
   33         for ( $j = 0; $j < $cols; $j++ )
   34         {
   35             my $format = (defined($object->rh_options->{display_format}))
   36                      ?   $object->[3]->{display_format} :
   37                     $Matrix::DEFAULT_FORMAT;
   38             $s .= sprintf($Matrix::DEFAULT_FORMAT, $object->[0][$i][$j]);
   39         }
   40         $s .= "]\n";
   41     }
   42     return($s);
   43 }
   44 
   45 =head4
   46 
   47   Method $matrix->rh_options
   48 
   49 =cut
   50 
   51 sub rh_options {
   52     my $self = shift;
   53     my $last_element = $#$self;
   54     $self->[$last_element] = {} unless defined($self->[3]); # not sure why this needs to be done
   55   $self->[$last_element];     # provides a reference to the options hash MEG
   56 }
   57 
   58 =head4
   59 
   60   Method $matrix->trace
   61 
   62   Returns: scalar which is the trace of the matrix.
   63 
   64 =cut
   65 
   66 sub trace {
   67   my $self = shift;
   68   my $rows = $self->[1];
   69   my $cols = $self->[2];
   70   warn "Can't take trace of non-square matrix " unless $rows == $cols;
   71   my $sum = 0;
   72   for( my $i = 0; $i<$rows;$i++) {
   73     $sum +=$self->[0][$i][$i];
   74   }
   75   $sum;
   76 }
   77 
   78 =head4
   79 
   80   Method $matrix->new_from_array_ref
   81 
   82 =cut
   83 
   84 sub new_from_array_ref {  # this will build a matrix or a row vector from  [a, b, c, ]
   85   my $class = shift;
   86   my $array = shift;
   87   my $rows = @$array;
   88   my $cols = @{$array->[0]};
   89   my $matrix = new Matrix($rows,$cols);
   90   $matrix->[0]=$array;
   91   $matrix;
   92 }
   93 
   94 =head4
   95 
   96   Method $matrix->array_ref
   97 
   98 =cut
   99 
  100 sub array_ref {
  101   my $this = shift;
  102   $this->[0];
  103 }
  104 
  105 =head4
  106 
  107   Method $matrix->list
  108 
  109 =cut
  110 
  111 sub list {           # this is used only for column vectors
  112   my $self = shift;
  113   my @list = ();
  114   warn "This only works with column vectors" unless $self->[2] == 1;
  115   my $rows = $self->[1];
  116   for(my $i=1; $i<=$rows; $i++) {
  117     push(@list, $self->element($i,1) );
  118   }
  119   @list;
  120 }
  121 
  122 =head4
  123 
  124   Method $matrix->new_from_list
  125 
  126 =cut
  127 
  128 sub new_from_list {   # this builds a row vector from an array
  129   my $class = shift;
  130   my @list = @_;
  131   my $cols = @list;
  132   my $rows = 1;
  133   my $matrix = new Matrix($rows, $cols);
  134   my $i=1;
  135   while(@list) {
  136       my $elem = shift(@list);
  137     $matrix->assign($i++,1, $elem);
  138   }
  139   $matrix;
  140 }
  141 
  142 =head4
  143 
  144   Method $matrix->new_row_matrix
  145 
  146 =cut
  147 
  148 sub new_row_matrix {   # this builds a row vector from an array
  149   my $class = shift;
  150   my @list = @_;
  151   my $cols = @list;
  152   my $rows = 1;
  153   my $matrix = new Matrix($rows, $cols);
  154   my $i=1;
  155   while(@list) {
  156       my $elem = shift(@list);
  157     $matrix->assign($i++,1, $elem);
  158   }
  159   $matrix;
  160 }
  161 
  162 =head4
  163 
  164   Method $matrix->proj
  165 
  166 =cut
  167 
  168 sub proj{
  169   my $self = shift;
  170   my ($vec) = @_;
  171   $self * $self ->proj_coeff($vec);
  172 }
  173 
  174 =head4
  175 
  176   Method $matrix->proj_coeff
  177 
  178 =cut
  179 
  180 sub proj_coeff{
  181   my $self= shift;
  182   my ($vec) = @_;
  183   warn 'The vector must be of type Matrix',ref($vec),"|" unless ref($vec) eq 'Matrix';
  184   my $lin_space_tr= ~ $self;
  185   my $matrix = $lin_space_tr * $self;
  186   $vec = $lin_space_tr*$vec;
  187   my $matrix_lr = $matrix->decompose_LR;
  188   my ($dim,$x_vector, $base_matrix) = $matrix_lr->solve_LR($vec);
  189   warn "A unique adapted answer could not be determined.  Possibly the parameters have coefficient zero.<br>  dim = $dim base_matrix is $base_matrix\n" if $dim;  # only print if the dim is not zero.
  190   $x_vector;
  191 }
  192 
  193 =head4
  194 
  195   Method $matrix->new_column_matrix
  196 
  197 =cut
  198 
  199 sub new_column_matrix {
  200   my $class = shift;
  201   my $vec = shift;
  202   warn "The argument to assign column must be a reference to an array" unless ref($vec) =~/ARRAY/;
  203   my $cols = 1;
  204   my $rows = @{$vec};
  205   my $matrix = new Matrix($rows,1);
  206   foreach my $i (1..$rows) {
  207     $matrix->assign($i,1,$vec->[$i-1]);
  208   }
  209   $matrix;
  210 }
  211 
  212 =head4
  213 
  214   This method takes an array of column vectors, or an array of arrays,
  215   and converts them to a matrix where each column is one of the previous
  216   vectors.
  217 
  218   Method $matrix->new_from_col_vecs
  219 
  220 =cut
  221 
  222 sub new_from_col_vecs
  223 {
  224   my $class = shift;
  225   my($vecs) = shift;
  226   my($rows,$cols);
  227 
  228   if(ref($vecs->[0])eq 'Matrix' ){
  229     ($rows,$cols) = (scalar($vecs->[0]->[1]),scalar(@$vecs));
  230   }else{
  231     ($rows,$cols) = (scalar(@{$vecs->[0]}),scalar(@$vecs));
  232   }
  233 
  234   my($i,$j);
  235       my $matrix = Matrix->new($rows,$cols);
  236 
  237     if(ref($vecs->[0])eq 'Matrix' ){
  238         for ( $i = 0; $i < $cols; $i++ )
  239         {
  240           for( $j = 0; $j < $rows; $j++ )
  241       {
  242               $matrix->[0][$j][$i] = $vecs->[$i][0][$j][0];
  243       }
  244         }
  245   }else{
  246     for ( $i = 0; $i < $cols; $i++ )
  247         {
  248           for( $j = 0; $j < $rows; $j++ )
  249       {
  250               $matrix->[0][$j][$i] = $vecs->[$i]->[$j];
  251       }
  252         }
  253   }
  254       return($matrix);
  255 }
  256 
  257 ######################################################################
  258 #  Modifications to MatrixReal.pm which allow use of complex entries
  259 ######################################################################
  260 
  261 =head3
  262 
  263   Overrides of MatrixReal which allow use of complex entries
  264 
  265 =cut
  266 
  267 =head4
  268 
  269   Method $matrix->new_from_col_vecs
  270 
  271 =cut
  272 
  273 sub cp  { # MEG  makes new copies of complex number
  274   my $z = shift;
  275   return $z unless ref($z);
  276   my $w = Complex1::cplx($z->Re,$z->Im);
  277   return $w;
  278 }
  279 
  280 =head4
  281 
  282   Method $matrix->copy
  283 
  284 =cut
  285 
  286 sub copy
  287 {
  288     croak "Usage: \$matrix1->copy(\$matrix2);"
  289       if (@_ != 2);
  290 
  291     my($matrix1,$matrix2) = @_;
  292     my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
  293     my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
  294     my($i,$j);
  295 
  296     croak "MatrixReal1::copy(): matrix size mismatch"
  297       unless (($rows1 == $rows2) && ($cols1 == $cols2));
  298 
  299     for ( $i = 0; $i < $rows1; $i++ )
  300     {
  301   my $r1 = []; # New array ref
  302   my $r2 = $matrix2->[0][$i];
  303   #@$r1 = @$r2; # Copy whole array directly  #MEG
  304   # if the array contains complex objects new objects must be created.
  305   foreach (@$r2) {
  306     push(@$r1,cp($_) );
  307   }
  308   $matrix1->[0][$i] = $r1;
  309     }
  310         $matrix1->[3] = $matrix2->[3]; # sign or option
  311     if (defined $matrix2->[4]) # is an LR decomposition matrix!
  312     {
  313     #    $matrix1->[3] = $matrix2->[3]; # $sign
  314         $matrix1->[4] = $matrix2->[4]; # $perm_row
  315         $matrix1->[5] = $matrix2->[5]; # $perm_col
  316         $matrix1->[6] = $matrix2->[6]; # $option
  317     }
  318 }
  319 ###################################################################
  320 
  321 # MEG added 6/25/03 to accomodate complex entries
  322 
  323 =head4
  324 
  325   Method $matrix->conj
  326 
  327 =cut
  328 
  329 sub conj {
  330   my $elem = shift;
  331     $elem = (ref($elem)) ? ($elem->conjugate) : $elem;
  332     $elem;
  333 }
  334 
  335 =head4
  336 
  337   Method $matrix->transpose
  338 
  339 =cut
  340 
  341 sub transpose
  342 {
  343     croak "Usage: \$matrix1->transpose(\$matrix2);"
  344       if (@_ != 2);
  345 
  346     my($matrix1,$matrix2) = @_;
  347     my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
  348     my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
  349 
  350     croak "MatrixReal1::transpose(): matrix size mismatch"
  351       unless (($rows1 == $cols2) && ($cols1 == $rows2));
  352 
  353     $matrix1->_undo_LR();
  354 
  355     if ($rows1 == $cols1)
  356     {
  357         # more complicated to make in-place possible!
  358         # # conj added by MEG
  359         for (my $i = 0; $i < $rows1; $i++)
  360         {
  361             for (my $j = ($i + 1); $j < $cols1; $j++)
  362             {
  363                 my $swap              = conj($matrix2->[0][$i][$j]);
  364                 $matrix1->[0][$i][$j] = conj($matrix2->[0][$j][$i]);
  365                 $matrix1->[0][$j][$i] = $swap;
  366             }
  367             $matrix1->[0][$i][$i] = conj($matrix2->[0][$i][$i]);
  368         }
  369 
  370     }
  371     else # ($rows1 != $cols1)
  372     {
  373         for (my $i = 0; $i < $rows1; $i++)
  374         {
  375             for (my $j = 0; $j < $cols1; $j++)
  376             {
  377                 $matrix1->[0][$i][$j] = conj($matrix2->[0][$j][$i]);
  378             }
  379         }
  380     }
  381     $matrix1;
  382 }
  383 
  384 =head4
  385 
  386   Method $matrix->decompose_LR
  387 
  388 =cut
  389 
  390 sub decompose_LR
  391 {
  392     croak "Usage: \$LR_matrix = \$matrix->decompose_LR();"
  393       if (@_ != 1);
  394 
  395     my($matrix) = @_;
  396     my($rows,$cols) = ($matrix->[1],$matrix->[2]);
  397     my($perm_row,$perm_col);
  398     my($row,$col,$max);
  399 #    my($i,$j,$k,$n);  #MEG
  400     my($i,$j,$k,);
  401     my($sign) = 1;
  402     my($swap);
  403     my($temp);
  404 #    FIXEME Why won't this work on non-square matrices?
  405 #    croak "MatrixReal1::decompose_LR(): matrix is not quadratic"
  406 #      unless ($rows == $cols);
  407     croak "MatrixReal1::decompose_LR(): matrix has more rows than columns"
  408       unless ($rows <= $cols);
  409 
  410     $temp = $matrix->new($rows,$cols);
  411     $temp->copy($matrix);
  412 #    $n = $rows;
  413     $perm_row = [ ];
  414     $perm_col = [ ];
  415     for ( $i = 0; $i < $rows; $i++ )  #i is a row number
  416     {
  417         $perm_row->[$i] = $i;
  418         $perm_col->[$i] = $i;
  419     }
  420     NONZERO:
  421     for ( $k = 0; $k < $rows; $k++ ) # use Gauss's algorithm:  #k is row number
  422     {
  423         # complete pivot-search:
  424 
  425         $max = 0;
  426         for ( $i = $k; $i < $cols; $i++ )   # i is column number
  427         {
  428             for ( $j = $k; $j < $cols; $j++ )
  429             {
  430                 if (($swap = abs($temp->[0][$i][$j])) > $max)
  431                 {
  432                     $max = $swap;
  433                     $row = $i;
  434                     $col = $j;
  435                 }
  436             }
  437         }
  438         last NONZERO if ($max == 0); # (all remaining elements are zero)
  439         if ($k != $row) # swap row $k and row $row:
  440         {
  441             $sign = -$sign;
  442             $swap             = $perm_row->[$k];
  443             $perm_row->[$k]   = $perm_row->[$row];
  444             $perm_row->[$row] = $swap;
  445             for ( $j = 0; $j < $cols; $j++ )   # j is a column number
  446             {
  447                 # (must run from 0 since L has to be swapped too!)
  448 
  449                 $swap                = $temp->[0][$k][$j];
  450                 $temp->[0][$k][$j]   = $temp->[0][$row][$j];
  451                 $temp->[0][$row][$j] = $swap;
  452             }
  453         }
  454         if ($k != $col) # swap column $k and column $col:
  455         {
  456             $sign = -$sign;
  457             $swap             = $perm_col->[$k];
  458             $perm_col->[$k]   = $perm_col->[$col];
  459             $perm_col->[$col] = $swap;
  460             for ( $i = 0; $i < $rows; $i++ )   #i is a row number
  461             {
  462                 $swap                = $temp->[0][$i][$k];
  463                 $temp->[0][$i][$k]   = $temp->[0][$i][$col];
  464                 $temp->[0][$i][$col] = $swap;
  465             }
  466         }
  467         for ( $i = ($k + 1); $i < $cols; $i++ )   # i is column number
  468         {
  469             # scan the remaining rows, add multiples of row $k to row $i:
  470 
  471             $swap = $temp->[0][$i][$k] / $temp->[0][$k][$k];
  472             if ($swap != 0)
  473             {
  474                 # calculate a row of matrix R:
  475 
  476                 for ( $j = ($k + 1); $j < $cols; $j++ )   #j is also a column number
  477                 {
  478                     $temp->[0][$i][$j] -= $temp->[0][$k][$j] * $swap;
  479                 }
  480 
  481                 # store matrix L in same matrix as R:
  482 
  483                 $temp->[0][$i][$k] = $swap;
  484             }
  485         }
  486     }
  487     my $rh_options = $temp->[3];
  488     $temp->[3] = $sign;
  489     $temp->[4] = $perm_row;
  490     $temp->[5] = $perm_col;
  491     $temp->[6] = $temp->[3];
  492     return($temp);
  493 }
  494 
  495 
  496 
  497 
  498 
  499 1;

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9