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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3373 - (view) (download) (as text)

1 : gage 3363 =head1 NAME
2 : apizer 1079
3 : gage 3363 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 : gage 3373 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 : sh002i 1050 package Matrix;
25 :     @Matrix::ISA = qw(MatrixReal1);
26 :    
27 : malsyned 1280 use Carp;
28 : sh002i 1050
29 :     $Matrix::DEFAULT_FORMAT = '% #-19.12E ';
30 :     # allows specification of the format
31 : gage 3360
32 :     =head4
33 :    
34 :     Method $matrix->_stringify()
35 :     -- overrides MatrixReal1 display mode
36 :    
37 :     =cut
38 :    
39 :    
40 : gage 3373 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 : sh002i 1050 my($rows,$cols) = ($object->[1],$object->[2]);
50 :     my($i,$j,$s);
51 : gage 3373
52 : sh002i 1050 $s = '';
53 :     for ( $i = 0; $i < $rows; $i++ )
54 :     {
55 :     $s .= "[ ";
56 :     for ( $j = 0; $j < $cols; $j++ )
57 : gage 3373 { #warn " i $i j $j ",$object->rh_options;
58 : apizer 1079 my $format = (defined($object->rh_options->{display_format}))
59 : gage 3373 ? $object->rh_options->{display_format} :
60 : sh002i 1050 $Matrix::DEFAULT_FORMAT;
61 : gage 3373 $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 : sh002i 1050 }
65 :     $s .= "]\n";
66 :     }
67 :     return($s);
68 :     }
69 :    
70 : gage 3373 # 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 : gage 3360 =head4
122 :    
123 :     Method $matrix->rh_options
124 :    
125 :     =cut
126 :    
127 : sh002i 1050 sub rh_options {
128 :     my $self = shift;
129 : gage 3373 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 : sh002i 1050 }
133 :    
134 : gage 3360 =head4
135 : sh002i 1050
136 : gage 3360 Method $matrix->trace
137 :    
138 :     Returns: scalar which is the trace of the matrix.
139 :    
140 :     =cut
141 :    
142 : sh002i 1050 sub trace {
143 :     my $self = shift;
144 :     my $rows = $self->[1];
145 :     my $cols = $self->[2];
146 : gage 1070 warn "Can't take trace of non-square matrix " unless $rows == $cols;
147 : sh002i 1050 my $sum = 0;
148 :     for( my $i = 0; $i<$rows;$i++) {
149 :     $sum +=$self->[0][$i][$i];
150 :     }
151 :     $sum;
152 :     }
153 : gage 3360
154 :     =head4
155 :    
156 :     Method $matrix->new_from_array_ref
157 :    
158 :     =cut
159 :    
160 : sh002i 1050 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 : gage 3360 =head4
171 :    
172 :     Method $matrix->array_ref
173 :    
174 :     =cut
175 :    
176 : sh002i 1050 sub array_ref {
177 :     my $this = shift;
178 :     $this->[0];
179 :     }
180 :    
181 : gage 3360 =head4
182 :    
183 :     Method $matrix->list
184 :    
185 :     =cut
186 :    
187 : sh002i 1050 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 : gage 3360
198 :     =head4
199 :    
200 :     Method $matrix->new_from_list
201 :    
202 :     =cut
203 :    
204 : sh002i 1050 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 : gage 3360
218 :     =head4
219 :    
220 :     Method $matrix->new_row_matrix
221 :    
222 :     =cut
223 :    
224 : sh002i 1050 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 : gage 3360
238 :     =head4
239 :    
240 :     Method $matrix->proj
241 :    
242 :     =cut
243 :    
244 : sh002i 1050 sub proj{
245 :     my $self = shift;
246 :     my ($vec) = @_;
247 : apizer 1079 $self * $self ->proj_coeff($vec);
248 :     }
249 : gage 3360
250 :     =head4
251 :    
252 :     Method $matrix->proj_coeff
253 :    
254 :     =cut
255 :    
256 : sh002i 1050 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 : gage 3360
269 :     =head4
270 :    
271 :     Method $matrix->new_column_matrix
272 :    
273 :     =cut
274 :    
275 : sh002i 1050 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 : apizer 1079 my $rows = @{$vec};
281 : sh002i 1050 my $matrix = new Matrix($rows,1);
282 :     foreach my $i (1..$rows) {
283 :     $matrix->assign($i,1,$vec->[$i-1]);
284 :     }
285 : gage 1070 $matrix;
286 :     }
287 : gage 3360
288 : gage 1070 =head4
289 :    
290 :     This method takes an array of column vectors, or an array of arrays,
291 : apizer 1079 and converts them to a matrix where each column is one of the previous
292 : gage 1070 vectors.
293 : gage 3360
294 :     Method $matrix->new_from_col_vecs
295 : gage 1070
296 :     =cut
297 :    
298 :     sub new_from_col_vecs
299 :     {
300 :     my $class = shift;
301 :     my($vecs) = shift;
302 :     my($rows,$cols);
303 : apizer 1079
304 : gage 1070 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 : apizer 1079
310 : gage 1070 my($i,$j);
311 : apizer 1079 my $matrix = Matrix->new($rows,$cols);
312 :    
313 : gage 1070 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 : apizer 1079 }
332 : gage 1285
333 : malsyned 1280 ######################################################################
334 :     # Modifications to MatrixReal.pm which allow use of complex entries
335 :     ######################################################################
336 : gage 1070
337 : gage 3360 =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 : malsyned 1280 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 : gage 3360
356 :     =head4
357 :    
358 :     Method $matrix->copy
359 :    
360 :     =cut
361 :    
362 : malsyned 1280 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 : gage 3360
399 :     =head4
400 :    
401 :     Method $matrix->conj
402 :    
403 :     =cut
404 :    
405 : malsyned 1280 sub conj {
406 :     my $elem = shift;
407 :     $elem = (ref($elem)) ? ($elem->conjugate) : $elem;
408 :     $elem;
409 :     }
410 : gage 3360
411 :     =head4
412 :    
413 :     Method $matrix->transpose
414 :    
415 :     =cut
416 :    
417 : malsyned 1280 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 : gage 3360 =head4
461 : malsyned 1280
462 : gage 3360 Method $matrix->decompose_LR
463 : gage 1285
464 : gage 3360 =cut
465 :    
466 : gage 1290 sub decompose_LR
467 :     {
468 :     croak "Usage: \$LR_matrix = \$matrix->decompose_LR();"
469 :     if (@_ != 1);
470 : gage 1285
471 : gage 1290 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 : gage 3373 my $rh_options = $matrix->[$MatrixReal1::OPTION_ENTRY];
480 : gage 3360 # FIXEME Why won't this work on non-square matrices?
481 : gage 1290 # croak "MatrixReal1::decompose_LR(): matrix is not quadratic"
482 :     # unless ($rows == $cols);
483 : gage 3373 # croak "MatrixReal1::decompose_LR(): matrix has more rows than columns"
484 :     # unless ($rows <= $cols);
485 : gage 1290
486 :     $temp = $matrix->new($rows,$cols);
487 :     $temp->copy($matrix);
488 :     # $n = $rows;
489 :     $perm_row = [ ];
490 :     $perm_col = [ ];
491 : gage 3373 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 : gage 1290 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 : gage 3373 for ( $i = $k; $i < $rows; $i++ ) # i is row number
500 : gage 1290 {
501 : gage 3373 for ( $j = $k; $j < $cols; $j++ ) #j is a col number
502 : gage 1290 {
503 :     if (($swap = abs($temp->[0][$i][$j])) > $max)
504 :     {
505 :     $max = $swap;
506 :     $row = $i;
507 :     $col = $j;
508 :     }
509 :     }
510 :     }
511 : gage 3373 # warn "max is $max row is $row and col is $col and k is $k";
512 : gage 1290 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 : gage 3373 { my $swap; # localize variable MEG
530 : gage 1290 $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 : gage 3373 for (my $i = ($k + 1); $i < $rows; $i++ ) # i is row number
542 : gage 1290 {
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 : gage 3373 for (my $j = ($k + 1); $j < $cols; $j++ ) #j is a column number
551 : gage 1290 {
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 : gage 3373 #my $rh_options = $temp->[3];
561 : gage 1290 $temp->[3] = $sign;
562 :     $temp->[4] = $perm_row;
563 :     $temp->[5] = $perm_col;
564 : gage 3373 $temp->[$MatrixReal1::OPTION_ENTRY] = $rh_options;
565 : gage 1290 return($temp);
566 :     }
567 :    
568 :    
569 :    
570 :    
571 :    
572 : gage 1070 1;

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9