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

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9