------------------------------------------------------------------------------ Factor linear equations Ax = b using LU decomposition A = LU where L is lower triangular matrix and U is upper triangular matrix. Input is square matrix A, integer array of pivot indices index[0->n-1], and size of square matrix n. Output factorization LU is in matrix A. If error is found, method returns 0.
| 483 | // of square matrix n. Output factorization LU is in matrix A. If error is |
| 484 | // found, method returns 0. |
| 485 | vtkTypeBool vtkMath::LUFactorLinearSystem(double** A, int* index, int size) |
| 486 | { |
| 487 | double scratch[VTK_MAX_SCRATCH_ARRAY_SIZE]; |
| 488 | double* scale = (size <= VTK_MAX_SCRATCH_ARRAY_SIZE ? scratch : new double[size]); |
| 489 | |
| 490 | int i, j, k; |
| 491 | int maxI = 0; |
| 492 | double largest, temp1, temp2, sum; |
| 493 | |
| 494 | // |
| 495 | // Loop over rows to get implicit scaling information |
| 496 | // |
| 497 | for (i = 0; i < size; ++i) |
| 498 | { |
| 499 | for (largest = 0.0, j = 0; j < size; ++j) |
| 500 | { |
| 501 | if ((temp2 = std::abs(A[i][j])) > largest) |
| 502 | { |
| 503 | largest = temp2; |
| 504 | } |
| 505 | } |
| 506 | |
| 507 | if (largest == 0.0) |
| 508 | { |
| 509 | vtkGenericWarningMacro(<< "Unable to factor linear system"); |
| 510 | if (size > VTK_MAX_SCRATCH_ARRAY_SIZE) |
| 511 | { |
| 512 | delete[] scale; |
| 513 | } |
| 514 | return 0; |
| 515 | } |
| 516 | scale[i] = 1.0 / largest; |
| 517 | } |
| 518 | // |
| 519 | // Loop over all columns using Crout's method |
| 520 | // |
| 521 | for (j = 0; j < size; ++j) |
| 522 | { |
| 523 | for (i = 0; i < j; ++i) |
| 524 | { |
| 525 | sum = A[i][j]; |
| 526 | for (k = 0; k < i; ++k) |
| 527 | { |
| 528 | sum -= A[i][k] * A[k][j]; |
| 529 | } |
| 530 | A[i][j] = sum; |
| 531 | } |
| 532 | // |
| 533 | // Begin search for largest pivot element |
| 534 | // |
| 535 | for (largest = 0.0, i = j; i < size; ++i) |
| 536 | { |
| 537 | sum = A[i][j]; |
| 538 | for (k = 0; k < j; ++k) |
| 539 | { |
| 540 | sum -= A[i][k] * A[k][j]; |
| 541 | } |
| 542 | A[i][j] = sum; |