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