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

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