| 629 | } |
| 630 | |
| 631 | @Override |
| 632 | public Matrix[] lup() |
| 633 | { |
| 634 | Matrix[] lup = new Matrix[3]; |
| 635 | |
| 636 | Matrix P = eye(rows()); |
| 637 | DenseMatrix L; |
| 638 | DenseMatrix U = this; |
| 639 | |
| 640 | //Initalization is a little wierd b/c we want to handle rectangular cases as well! |
| 641 | if(rows() > cols())//In this case, we will be changing U before returning it (have to make it smaller, but we can still avoid allocating extra space |
| 642 | L = new DenseMatrix(rows(), cols()); |
| 643 | else |
| 644 | L = new DenseMatrix(rows(), rows()); |
| 645 | |
| 646 | for(int i = 0; i < U.rows(); i++) |
| 647 | { |
| 648 | //If rectangular, we still need to loop through to update ther est of L - even though we wont make many other changes |
| 649 | if(i < U.cols()) |
| 650 | { |
| 651 | //Partial pivoting, find the largest value in this colum and move it to the top! |
| 652 | //Find the largest magintude value in the colum k, row j |
| 653 | int largestRow = i; |
| 654 | double largestVal = Math.abs(U.matrix[i][i]); |
| 655 | for (int j = i + 1; j < U.rows(); j++) |
| 656 | { |
| 657 | double rowJLeadVal = Math.abs(U.matrix[j][i]); |
| 658 | if (rowJLeadVal > largestVal) |
| 659 | { |
| 660 | largestRow = j; |
| 661 | largestVal = rowJLeadVal; |
| 662 | } |
| 663 | } |
| 664 | |
| 665 | //SWAP! |
| 666 | U.swapRows(largestRow, i); |
| 667 | P.swapRows(largestRow, i); |
| 668 | L.swapRows(largestRow, i); |
| 669 | |
| 670 | L.matrix[i][i] = 1; |
| 671 | } |
| 672 | |
| 673 | //Seting up L |
| 674 | for(int k = 0; k < Math.min(i, U.cols()); k++) |
| 675 | { |
| 676 | double tmp = U.matrix[i][k]/U.matrix[k][k]; |
| 677 | L.matrix[i][k] = Double.isNaN(tmp) ? 0.0 : tmp; |
| 678 | U.matrix[i][k] = 0; |
| 679 | |
| 680 | for(int j = k+1; j < U.cols(); j++) |
| 681 | { |
| 682 | U.matrix[i][j] -= L.matrix[i][k]*U.matrix[k][j]; |
| 683 | } |
| 684 | } |
| 685 | } |
| 686 | |
| 687 | |
| 688 | if(rows() > cols())//Clean up! |