01: package com.vividsolutions.jump.geom;
02:
03: import com.vividsolutions.jump.I18N;
04:
05: /**
06: * Implements some 2D matrix operations
07: * (in particular, solving systems of linear equations).
08: *
09: * @author Martin Davis
10: *
11: */
12: public class Matrix {
13: private static void swapRows(double[][] m, int i, int j) {
14: if (i == j)
15: return;
16: for (int col = 0; col < m[0].length; col++) {
17: double temp = m[i][col];
18: m[i][col] = m[j][col];
19: m[j][col] = temp;
20: }
21: }
22:
23: private static void swapRows(double[] m, int i, int j) {
24: if (i == j)
25: return;
26: double temp = m[i];
27: m[i] = m[j];
28: m[j] = temp;
29: }
30:
31: /**
32: * Solves a system of equations using Gaussian Elimination.
33: * In order to avoid overhead the algorithm runs in-place
34: * on A - if A should not be modified the client must supply a copy.
35: *
36: * @param a an nxn matrix in row/column order )modified by this method)
37: * @param b a vector of length n
38: *
39: * @return a vector containing the solution (if any)
40: * @return null if the system has no or no unique solution
41: *
42: * @throws IllegalArgumentException if the matrix is the wrong size
43: */
44: public static double[] solve(double[][] a, double[] b) {
45: int n = b.length;
46: if (a.length != n || a[0].length != n)
47: throw new IllegalArgumentException(
48: I18N
49: .get("jump.geom.Matrix.Matrix-A-is-incorrectly-sized"));
50:
51: // Use Gaussian Elimination with partial pivoting.
52: // Iterate over each row
53: for (int i = 0; i < n; i++) {
54: // Find the largest pivot in the rows below the current one.
55: int maxElementRow = i;
56: for (int j = i + 1; j < n; j++)
57: if (Math.abs(a[j][i]) > Math.abs(a[maxElementRow][i]))
58: maxElementRow = j;
59:
60: if (a[maxElementRow][i] == 0.0)
61: return null;
62:
63: // Exchange current row and maxElementRow in A and b.
64: swapRows(a, i, maxElementRow);
65: swapRows(b, i, maxElementRow);
66:
67: // Eliminate using row i
68: for (int j = i + 1; j < n; j++) {
69: double rowFactor = a[j][i] / a[i][i];
70: for (int k = n - 1; k >= i; k--)
71: a[j][k] -= a[i][k] * rowFactor;
72: b[j] -= b[i] * rowFactor;
73: }
74: }
75:
76: /**
77: * A is now (virtually) in upper-triangular form.
78: * The solution vector is determined by back-substitution.
79: */
80: double[] solution = new double[n];
81: for (int j = n - 1; j >= 0; j--) {
82: double t = 0.0;
83: for (int k = j + 1; k < n; k++)
84: t += a[j][k] * solution[k];
85: solution[j] = (b[j] - t) / a[j][j];
86: }
87: return solution;
88: }
89:
90: }
|