[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 1290 - (download) (as text) (annotate)
Thu Jun 26 22:43:15 2003 UTC (16 years, 6 months ago) by gage
File size: 10295 byte(s)
Merge the time optimized branch and the main branch.  I forgot one subroutine
--Mike

    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 sub _stringify
   14 {
   15     my($object,$argument,$flag) = @_;
   16 #   my($name) = '""'; #&_trace($name,$object,$argument,$flag);
   17     my($rows,$cols) = ($object->[1],$object->[2]);
   18     my($i,$j,$s);
   19 
   20     $s = '';
   21     for ( $i = 0; $i < $rows; $i++ )
   22     {
   23         $s .= "[ ";
   24         for ( $j = 0; $j < $cols; $j++ )
   25         {
   26             my $format = (defined($object->rh_options->{display_format}))
   27                      ?   $object->[3]->{display_format} :
   28                     $Matrix::DEFAULT_FORMAT;
   29             $s .= sprintf($Matrix::DEFAULT_FORMAT, $object->[0][$i][$j]);
   30         }
   31         $s .= "]\n";
   32     }
   33     return($s);
   34 }
   35 
   36 sub rh_options {
   37     my $self = shift;
   38     my $last_element = $#$self;
   39     $self->[$last_element] = {} unless defined($self->[3]); # not sure why this needs to be done
   40   $self->[$last_element];     # provides a reference to the options hash MEG
   41 }
   42 
   43 
   44 sub trace {
   45   my $self = shift;
   46   my $rows = $self->[1];
   47   my $cols = $self->[2];
   48   warn "Can't take trace of non-square matrix " unless $rows == $cols;
   49   my $sum = 0;
   50   for( my $i = 0; $i<$rows;$i++) {
   51     $sum +=$self->[0][$i][$i];
   52   }
   53   $sum;
   54 }
   55 sub new_from_array_ref {  # this will build a matrix or a row vector from  [a, b, c, ]
   56   my $class = shift;
   57   my $array = shift;
   58   my $rows = @$array;
   59   my $cols = @{$array->[0]};
   60   my $matrix = new Matrix($rows,$cols);
   61   $matrix->[0]=$array;
   62   $matrix;
   63 }
   64 
   65 sub array_ref {
   66   my $this = shift;
   67   $this->[0];
   68 }
   69 
   70 sub list {           # this is used only for column vectors
   71   my $self = shift;
   72   my @list = ();
   73   warn "This only works with column vectors" unless $self->[2] == 1;
   74   my $rows = $self->[1];
   75   for(my $i=1; $i<=$rows; $i++) {
   76     push(@list, $self->element($i,1) );
   77   }
   78   @list;
   79 }
   80 sub new_from_list {   # this builds a row vector from an array
   81   my $class = shift;
   82   my @list = @_;
   83   my $cols = @list;
   84   my $rows = 1;
   85   my $matrix = new Matrix($rows, $cols);
   86   my $i=1;
   87   while(@list) {
   88       my $elem = shift(@list);
   89     $matrix->assign($i++,1, $elem);
   90   }
   91   $matrix;
   92 }
   93 sub new_row_matrix {   # this builds a row vector from an array
   94   my $class = shift;
   95   my @list = @_;
   96   my $cols = @list;
   97   my $rows = 1;
   98   my $matrix = new Matrix($rows, $cols);
   99   my $i=1;
  100   while(@list) {
  101       my $elem = shift(@list);
  102     $matrix->assign($i++,1, $elem);
  103   }
  104   $matrix;
  105 }
  106 sub proj{
  107   my $self = shift;
  108   my ($vec) = @_;
  109   $self * $self ->proj_coeff($vec);
  110 }
  111 sub proj_coeff{
  112   my $self= shift;
  113   my ($vec) = @_;
  114   warn 'The vector must be of type Matrix',ref($vec),"|" unless ref($vec) eq 'Matrix';
  115   my $lin_space_tr= ~ $self;
  116   my $matrix = $lin_space_tr * $self;
  117   $vec = $lin_space_tr*$vec;
  118   my $matrix_lr = $matrix->decompose_LR;
  119   my ($dim,$x_vector, $base_matrix) = $matrix_lr->solve_LR($vec);
  120   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.
  121   $x_vector;
  122 }
  123 sub new_column_matrix {
  124   my $class = shift;
  125   my $vec = shift;
  126   warn "The argument to assign column must be a reference to an array" unless ref($vec) =~/ARRAY/;
  127   my $cols = 1;
  128   my $rows = @{$vec};
  129   my $matrix = new Matrix($rows,1);
  130   foreach my $i (1..$rows) {
  131     $matrix->assign($i,1,$vec->[$i-1]);
  132   }
  133   $matrix;
  134 }
  135 =head4
  136 
  137   This method takes an array of column vectors, or an array of arrays,
  138   and converts them to a matrix where each column is one of the previous
  139   vectors.
  140 
  141 =cut
  142 
  143 sub new_from_col_vecs
  144 {
  145   my $class = shift;
  146   my($vecs) = shift;
  147   my($rows,$cols);
  148 
  149   if(ref($vecs->[0])eq 'Matrix' ){
  150     ($rows,$cols) = (scalar($vecs->[0]->[1]),scalar(@$vecs));
  151   }else{
  152     ($rows,$cols) = (scalar(@{$vecs->[0]}),scalar(@$vecs));
  153   }
  154 
  155   my($i,$j);
  156       my $matrix = Matrix->new($rows,$cols);
  157 
  158     if(ref($vecs->[0])eq 'Matrix' ){
  159         for ( $i = 0; $i < $cols; $i++ )
  160         {
  161           for( $j = 0; $j < $rows; $j++ )
  162       {
  163               $matrix->[0][$j][$i] = $vecs->[$i][0][$j][0];
  164       }
  165         }
  166   }else{
  167     for ( $i = 0; $i < $cols; $i++ )
  168         {
  169           for( $j = 0; $j < $rows; $j++ )
  170       {
  171               $matrix->[0][$j][$i] = $vecs->[$i]->[$j];
  172       }
  173         }
  174   }
  175       return($matrix);
  176 }
  177 
  178 ######################################################################
  179 #  Modifications to MatrixReal.pm which allow use of complex entries
  180 ######################################################################
  181 
  182 sub cp  { # MEG  makes new copies of complex number
  183   my $z = shift;
  184   return $z unless ref($z);
  185   my $w = Complex1::cplx($z->Re,$z->Im);
  186   return $w;
  187 }
  188 sub copy
  189 {
  190     croak "Usage: \$matrix1->copy(\$matrix2);"
  191       if (@_ != 2);
  192 
  193     my($matrix1,$matrix2) = @_;
  194     my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
  195     my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
  196     my($i,$j);
  197 
  198     croak "MatrixReal1::copy(): matrix size mismatch"
  199       unless (($rows1 == $rows2) && ($cols1 == $cols2));
  200 
  201     for ( $i = 0; $i < $rows1; $i++ )
  202     {
  203   my $r1 = []; # New array ref
  204   my $r2 = $matrix2->[0][$i];
  205   #@$r1 = @$r2; # Copy whole array directly  #MEG
  206   # if the array contains complex objects new objects must be created.
  207   foreach (@$r2) {
  208     push(@$r1,cp($_) );
  209   }
  210   $matrix1->[0][$i] = $r1;
  211     }
  212         $matrix1->[3] = $matrix2->[3]; # sign or option
  213     if (defined $matrix2->[4]) # is an LR decomposition matrix!
  214     {
  215     #    $matrix1->[3] = $matrix2->[3]; # $sign
  216         $matrix1->[4] = $matrix2->[4]; # $perm_row
  217         $matrix1->[5] = $matrix2->[5]; # $perm_col
  218         $matrix1->[6] = $matrix2->[6]; # $option
  219     }
  220 }
  221 ###################################################################
  222 
  223 # MEG added 6/25/03 to accomodate complex entries
  224 sub conj {
  225   my $elem = shift;
  226     $elem = (ref($elem)) ? ($elem->conjugate) : $elem;
  227     $elem;
  228 }
  229 sub transpose
  230 {
  231     croak "Usage: \$matrix1->transpose(\$matrix2);"
  232       if (@_ != 2);
  233 
  234     my($matrix1,$matrix2) = @_;
  235     my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]);
  236     my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]);
  237 
  238     croak "MatrixReal1::transpose(): matrix size mismatch"
  239       unless (($rows1 == $cols2) && ($cols1 == $rows2));
  240 
  241     $matrix1->_undo_LR();
  242 
  243     if ($rows1 == $cols1)
  244     {
  245         # more complicated to make in-place possible!
  246         # # conj added by MEG
  247         for (my $i = 0; $i < $rows1; $i++)
  248         {
  249             for (my $j = ($i + 1); $j < $cols1; $j++)
  250             {
  251                 my $swap              = conj($matrix2->[0][$i][$j]);
  252                 $matrix1->[0][$i][$j] = conj($matrix2->[0][$j][$i]);
  253                 $matrix1->[0][$j][$i] = $swap;
  254             }
  255             $matrix1->[0][$i][$i] = conj($matrix2->[0][$i][$i]);
  256         }
  257 
  258     }
  259     else # ($rows1 != $cols1)
  260     {
  261         for (my $i = 0; $i < $rows1; $i++)
  262         {
  263             for (my $j = 0; $j < $cols1; $j++)
  264             {
  265                 $matrix1->[0][$i][$j] = conj($matrix2->[0][$j][$i]);
  266             }
  267         }
  268     }
  269     $matrix1;
  270 }
  271 
  272 
  273 
  274 sub decompose_LR
  275 {
  276     croak "Usage: \$LR_matrix = \$matrix->decompose_LR();"
  277       if (@_ != 1);
  278 
  279     my($matrix) = @_;
  280     my($rows,$cols) = ($matrix->[1],$matrix->[2]);
  281     my($perm_row,$perm_col);
  282     my($row,$col,$max);
  283 #    my($i,$j,$k,$n);  #MEG
  284     my($i,$j,$k,);
  285     my($sign) = 1;
  286     my($swap);
  287     my($temp);
  288 #    Why won't this work on non-square matrices?
  289 #    croak "MatrixReal1::decompose_LR(): matrix is not quadratic"
  290 #      unless ($rows == $cols);
  291     croak "MatrixReal1::decompose_LR(): matrix has more rows than columns"
  292       unless ($rows <= $cols);
  293 
  294     $temp = $matrix->new($rows,$cols);
  295     $temp->copy($matrix);
  296 #    $n = $rows;
  297     $perm_row = [ ];
  298     $perm_col = [ ];
  299     for ( $i = 0; $i < $rows; $i++ )  #i is a row number
  300     {
  301         $perm_row->[$i] = $i;
  302         $perm_col->[$i] = $i;
  303     }
  304     NONZERO:
  305     for ( $k = 0; $k < $rows; $k++ ) # use Gauss's algorithm:  #k is row number
  306     {
  307         # complete pivot-search:
  308 
  309         $max = 0;
  310         for ( $i = $k; $i < $cols; $i++ )   # i is column number
  311         {
  312             for ( $j = $k; $j < $cols; $j++ )
  313             {
  314                 if (($swap = abs($temp->[0][$i][$j])) > $max)
  315                 {
  316                     $max = $swap;
  317                     $row = $i;
  318                     $col = $j;
  319                 }
  320             }
  321         }
  322         last NONZERO if ($max == 0); # (all remaining elements are zero)
  323         if ($k != $row) # swap row $k and row $row:
  324         {
  325             $sign = -$sign;
  326             $swap             = $perm_row->[$k];
  327             $perm_row->[$k]   = $perm_row->[$row];
  328             $perm_row->[$row] = $swap;
  329             for ( $j = 0; $j < $cols; $j++ )   # j is a column number
  330             {
  331                 # (must run from 0 since L has to be swapped too!)
  332 
  333                 $swap                = $temp->[0][$k][$j];
  334                 $temp->[0][$k][$j]   = $temp->[0][$row][$j];
  335                 $temp->[0][$row][$j] = $swap;
  336             }
  337         }
  338         if ($k != $col) # swap column $k and column $col:
  339         {
  340             $sign = -$sign;
  341             $swap             = $perm_col->[$k];
  342             $perm_col->[$k]   = $perm_col->[$col];
  343             $perm_col->[$col] = $swap;
  344             for ( $i = 0; $i < $rows; $i++ )   #i is a row number
  345             {
  346                 $swap                = $temp->[0][$i][$k];
  347                 $temp->[0][$i][$k]   = $temp->[0][$i][$col];
  348                 $temp->[0][$i][$col] = $swap;
  349             }
  350         }
  351         for ( $i = ($k + 1); $i < $cols; $i++ )   # i is column number
  352         {
  353             # scan the remaining rows, add multiples of row $k to row $i:
  354 
  355             $swap = $temp->[0][$i][$k] / $temp->[0][$k][$k];
  356             if ($swap != 0)
  357             {
  358                 # calculate a row of matrix R:
  359 
  360                 for ( $j = ($k + 1); $j < $cols; $j++ )   #j is also a column number
  361                 {
  362                     $temp->[0][$i][$j] -= $temp->[0][$k][$j] * $swap;
  363                 }
  364 
  365                 # store matrix L in same matrix as R:
  366 
  367                 $temp->[0][$i][$k] = $swap;
  368             }
  369         }
  370     }
  371     my $rh_options = $temp->[3];
  372     $temp->[3] = $sign;
  373     $temp->[4] = $perm_row;
  374     $temp->[5] = $perm_col;
  375     $temp->[6] = $temp->[3];
  376     return($temp);
  377 }
  378 
  379 
  380 
  381 
  382 
  383 1;

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9