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