Parent Directory
| 
Revision Log
Merge the time optimized branch and the main branch. I forgot one subroutine --Mike
1 2 BEGIN { 3 be_strict(); # an alias for use strict. This means that all global variable must contain main:: as a prefix. 4 5 } 6 package Matrix; 7 @Matrix::ISA = qw(MatrixReal1); 8 9 use Carp; 10 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 20 $s = ''; 21 for ( $i = 0; $i < $rows; $i++ ) 22 { 23 $s .= "[ "; 24 for ( $j = 0; $j < $cols; $j++ ) 25 { 26 my $format = (defined($object->rh_options->{display_format})) 27 ? $object->[3]->{display_format} : 28 $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 warn "Can't take trace of non-square matrix " unless $rows == $cols; 49 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 $self * $self ->proj_coeff($vec); 110 } 111 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 my $rows = @{$vec}; 129 my $matrix = new Matrix($rows,1); 130 foreach my $i (1..$rows) { 131 $matrix->assign($i,1,$vec->[$i-1]); 132 } 133 $matrix; 134 } 135 =head4 136 137 This method takes an array of column vectors, or an array of arrays, 138 and converts them to a matrix where each column is one of the previous 139 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 149 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 155 my($i,$j); 156 my $matrix = Matrix->new($rows,$cols); 157 158 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 } 177 178 ###################################################################### 179 # Modifications to MatrixReal.pm which allow use of complex entries 180 ###################################################################### 181 182 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 274 sub decompose_LR 275 { 276 croak "Usage: \$LR_matrix = \$matrix->decompose_LR();" 277 if (@_ != 1); 278 279 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 1;
| aubreyja at gmail dot com | ViewVC Help |
| Powered by ViewVC 1.0.9 |