| 1 | package de.uka.ipd.sdq.reliability.solver; |
| 2 | |
| 3 | /////////////////////////////////////////////////////////////////////////// |
| 4 | // // |
| 5 | // Program file name: Inverse.java // |
| 6 | // // |
| 7 | // � Tao Pang 2006 // |
| 8 | // // |
| 9 | // Last modified: January 18, 2006 // |
| 10 | // // |
| 11 | // (1) This Java program is part of the book, "An Introduction to // |
| 12 | // Computational Physics, 2nd Edition," written by Tao Pang and // |
| 13 | // published by Cambridge University Press on January 19, 2006. // |
| 14 | // // |
| 15 | // (2) No warranties, express or implied, are made for this program. // |
| 16 | // // |
| 17 | /////////////////////////////////////////////////////////////////////////// |
| 18 | |
| 19 | // An example of performing matrix inversion through the |
| 20 | // partial-pivoting Gaussian elimination. |
| 21 | |
| 22 | public class Inverse { |
| 23 | |
| 24 | public static void gaussian(double a[][], int index[]) { |
| 25 | int n = index.length; |
| 26 | double c[] = new double[n]; |
| 27 | |
| 28 | // Initialize the index |
| 29 | for (int i = 0; i < n; ++i) |
| 30 | index[i] = i; |
| 31 | |
| 32 | // Find the rescaling factors, one from each row |
| 33 | for (int i = 0; i < n; ++i) { |
| 34 | double c1 = 0; |
| 35 | for (int j = 0; j < n; ++j) { |
| 36 | double c0 = Math.abs(a[i][j]); |
| 37 | if (c0 > c1) |
| 38 | c1 = c0; |
| 39 | } |
| 40 | c[i] = c1; |
| 41 | } |
| 42 | |
| 43 | // Search the pivoting element from each column |
| 44 | int k = 0; |
| 45 | for (int j = 0; j < n - 1; ++j) { |
| 46 | double pi1 = 0; |
| 47 | for (int i = j; i < n; ++i) { |
| 48 | double pi0 = Math.abs(a[index[i]][j]); |
| 49 | pi0 /= c[index[i]]; |
| 50 | if (pi0 > pi1) { |
| 51 | pi1 = pi0; |
| 52 | k = i; |
| 53 | } |
| 54 | } |
| 55 | |
| 56 | // Interchange rows according to the pivoting order |
| 57 | int itmp = index[j]; |
| 58 | index[j] = index[k]; |
| 59 | index[k] = itmp; |
| 60 | for (int i = j + 1; i < n; ++i) { |
| 61 | double pj = a[index[i]][j] / a[index[j]][j]; |
| 62 | |
| 63 | // Record pivoting ratios below the diagonal |
| 64 | a[index[i]][j] = pj; |
| 65 | |
| 66 | // Modify other elements accordingly |
| 67 | for (int l = j + 1; l < n; ++l) |
| 68 | a[index[i]][l] -= pj * a[index[j]][l]; |
| 69 | } |
| 70 | } |
| 71 | } |
| 72 | |
| 73 | // Method to carry out the partial-pivoting Gaussian |
| 74 | // elimination. Here index[] stores pivoting order. |
| 75 | |
| 76 | public static double[][] invert(double a[][]) { |
| 77 | int n = a.length; |
| 78 | double x[][] = new double[n][n]; |
| 79 | double b[][] = new double[n][n]; |
| 80 | int index[] = new int[n]; |
| 81 | for (int i = 0; i < n; ++i) |
| 82 | b[i][i] = 1; |
| 83 | |
| 84 | // Transform the matrix into an upper triangle |
| 85 | gaussian(a, index); |
| 86 | |
| 87 | // Update the matrix b[i][j] with the ratios stored |
| 88 | for (int i = 0; i < n - 1; ++i) |
| 89 | for (int j = i + 1; j < n; ++j) |
| 90 | for (int k = 0; k < n; ++k) |
| 91 | b[index[j]][k] -= a[index[j]][i] * b[index[i]][k]; |
| 92 | |
| 93 | // Perform backward substitutions |
| 94 | for (int i = 0; i < n; ++i) { |
| 95 | x[n - 1][i] = b[index[n - 1]][i] / a[index[n - 1]][n - 1]; |
| 96 | for (int j = n - 2; j >= 0; --j) { |
| 97 | x[j][i] = b[index[j]][i]; |
| 98 | for (int k = j + 1; k < n; ++k) { |
| 99 | x[j][i] -= a[index[j]][k] * x[k][i]; |
| 100 | } |
| 101 | x[j][i] /= a[index[j]][j]; |
| 102 | } |
| 103 | } |
| 104 | return x; |
| 105 | } |
| 106 | } |