[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 3360 - (view) (download) (as text)

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

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9