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

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9