[system] / trunk / pg / lib / Matrix.pm Repository:
ViewVC logotype

Diff of /trunk/pg/lib/Matrix.pm

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

Revision 1285 Revision 1290
269 $matrix1; 269 $matrix1;
270} 270}
271 271
272 272
273 273
274sub 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
274 382
2751; 3831;

Legend:
Removed from v.1285  
changed lines
  Added in v.1290

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9