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

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9