Parent Directory
| 
Revision Log
The beginnings of making decompose_LR work for non-square matrices --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 # Modifications to MatrixReal.pm which allow use of complex entries 179 ###################################################################### 180 181 sub cp { # MEG makes new copies of complex number 182 my $z = shift; 183 return $z unless ref($z); 184 my $w = Complex1::cplx($z->Re,$z->Im); 185 return $w; 186 } 187 sub copy 188 { 189 croak "Usage: \$matrix1->copy(\$matrix2);" 190 if (@_ != 2); 191 192 my($matrix1,$matrix2) = @_; 193 my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]); 194 my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]); 195 my($i,$j); 196 197 croak "MatrixReal1::copy(): matrix size mismatch" 198 unless (($rows1 == $rows2) && ($cols1 == $cols2)); 199 200 for ( $i = 0; $i < $rows1; $i++ ) 201 { 202 my $r1 = []; # New array ref 203 my $r2 = $matrix2->[0][$i]; 204 #@$r1 = @$r2; # Copy whole array directly #MEG 205 # if the array contains complex objects new objects must be created. 206 foreach (@$r2) { 207 push(@$r1,cp($_) ); 208 } 209 $matrix1->[0][$i] = $r1; 210 } 211 $matrix1->[3] = $matrix2->[3]; # sign or option 212 if (defined $matrix2->[4]) # is an LR decomposition matrix! 213 { 214 # $matrix1->[3] = $matrix2->[3]; # $sign 215 $matrix1->[4] = $matrix2->[4]; # $perm_row 216 $matrix1->[5] = $matrix2->[5]; # $perm_col 217 $matrix1->[6] = $matrix2->[6]; # $option 218 } 219 } 220 ################################################################### 221 222 # MEG added 6/25/03 to accomodate complex entries 223 sub conj { 224 my $elem = shift; 225 $elem = (ref($elem)) ? ($elem->conjugate) : $elem; 226 $elem; 227 } 228 sub transpose 229 { 230 croak "Usage: \$matrix1->transpose(\$matrix2);" 231 if (@_ != 2); 232 233 my($matrix1,$matrix2) = @_; 234 my($rows1,$cols1) = ($matrix1->[1],$matrix1->[2]); 235 my($rows2,$cols2) = ($matrix2->[1],$matrix2->[2]); 236 237 croak "MatrixReal1::transpose(): matrix size mismatch" 238 unless (($rows1 == $cols2) && ($cols1 == $rows2)); 239 240 $matrix1->_undo_LR(); 241 242 if ($rows1 == $cols1) 243 { 244 # more complicated to make in-place possible! 245 # # conj added by MEG 246 for (my $i = 0; $i < $rows1; $i++) 247 { 248 for (my $j = ($i + 1); $j < $cols1; $j++) 249 { 250 my $swap = conj($matrix2->[0][$i][$j]); 251 $matrix1->[0][$i][$j] = conj($matrix2->[0][$j][$i]); 252 $matrix1->[0][$j][$i] = $swap; 253 } 254 $matrix1->[0][$i][$i] = conj($matrix2->[0][$i][$i]); 255 } 256 257 } 258 else # ($rows1 != $cols1) 259 { 260 for (my $i = 0; $i < $rows1; $i++) 261 { 262 for (my $j = 0; $j < $cols1; $j++) 263 { 264 $matrix1->[0][$i][$j] = conj($matrix2->[0][$j][$i]); 265 } 266 } 267 } 268 $matrix1; 269 } 270 271 sub decompose_LR 272 { 273 croak "Usage: \$LR_matrix = \$matrix->decompose_LR();" 274 if (@_ != 1); 275 276 my($matrix) = @_; 277 my($rows,$cols) = ($matrix->[1],$matrix->[2]); 278 my($perm_row,$perm_col); 279 my($row,$col,$max); 280 # my($i,$j,$k,$n); #MEG 281 my($i,$j,$k,); 282 my($sign) = 1; 283 my($swap); 284 my($temp); 285 # Why won't this work on non-square matrices? 286 # croak "MatrixReal1::decompose_LR(): matrix is not quadratic" 287 # unless ($rows == $cols); 288 croak "MatrixReal1::decompose_LR(): matrix has more rows than columns" 289 unless ($rows <= $cols); 290 291 $temp = $matrix->new($rows,$cols); 292 $temp->copy($matrix); 293 # $n = $rows; 294 $perm_row = [ ]; 295 $perm_col = [ ]; 296 for ( $i = 0; $i < $rows; $i++ ) #i is a row number 297 { 298 $perm_row->[$i] = $i; 299 $perm_col->[$i] = $i; 300 } 301 NONZERO: 302 for ( $k = 0; $k < $rows; $k++ ) # use Gauss's algorithm: #k is row number 303 { 304 # complete pivot-search: 305 306 $max = 0; 307 for ( $i = $k; $i < $cols; $i++ ) # i is column number 308 { 309 for ( $j = $k; $j < $cols; $j++ ) 310 { 311 if (($swap = abs($temp->[0][$i][$j])) > $max) 312 { 313 $max = $swap; 314 $row = $i; 315 $col = $j; 316 } 317 } 318 } 319 last NONZERO if ($max == 0); # (all remaining elements are zero) 320 if ($k != $row) # swap row $k and row $row: 321 { 322 $sign = -$sign; 323 $swap = $perm_row->[$k]; 324 $perm_row->[$k] = $perm_row->[$row]; 325 $perm_row->[$row] = $swap; 326 for ( $j = 0; $j < $cols; $j++ ) # j is a column number 327 { 328 # (must run from 0 since L has to be swapped too!) 329 330 $swap = $temp->[0][$k][$j]; 331 $temp->[0][$k][$j] = $temp->[0][$row][$j]; 332 $temp->[0][$row][$j] = $swap; 333 } 334 } 335 if ($k != $col) # swap column $k and column $col: 336 { 337 $sign = -$sign; 338 $swap = $perm_col->[$k]; 339 $perm_col->[$k] = $perm_col->[$col]; 340 $perm_col->[$col] = $swap; 341 for ( $i = 0; $i < $rows; $i++ ) #i is a row number 342 { 343 $swap = $temp->[0][$i][$k]; 344 $temp->[0][$i][$k] = $temp->[0][$i][$col]; 345 $temp->[0][$i][$col] = $swap; 346 } 347 } 348 for ( $i = ($k + 1); $i < $cols; $i++ ) # i is column number 349 { 350 # scan the remaining rows, add multiples of row $k to row $i: 351 352 $swap = $temp->[0][$i][$k] / $temp->[0][$k][$k]; 353 if ($swap != 0) 354 { 355 # calculate a row of matrix R: 356 357 for ( $j = ($k + 1); $j < $cols; $j++ ) #j is also a column number 358 { 359 $temp->[0][$i][$j] -= $temp->[0][$k][$j] * $swap; 360 } 361 362 # store matrix L in same matrix as R: 363 364 $temp->[0][$i][$k] = $swap; 365 } 366 } 367 } 368 my $rh_options = $temp->[3]; 369 $temp->[3] = $sign; 370 $temp->[4] = $perm_row; 371 $temp->[5] = $perm_col; 372 $temp->[6] = $temp->[3]; 373 return($temp); 374 } 375 376 377 1;
| aubreyja at gmail dot com | ViewVC Help |
| Powered by ViewVC 1.0.9 |