[system] / trunk / pg / lib / Matrix.pm Repository: Repository Listing bbplugincoursesdistsnplrochestersystemwww

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

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