com.vividsolutions.jts.math.Matrix Maven / Gradle / Ivy
The newest version!
/*
* The JTS Topology Suite is a collection of Java classes that
* implement the fundamental operations required to validate a given
* geo-spatial data set to a known topological specification.
*
* Copyright (C) 2001 Vivid Solutions
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with this library; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
* For more information, contact:
*
* Vivid Solutions
* Suite #1A
* 2328 Government Street
* Victoria BC V8T 5G5
* Canada
*
* (250)385-6040
* www.vividsolutions.com
*/
package com.vividsolutions.jts.math;
/**
* Implements some 2D matrix operations
* (in particular, solving systems of linear equations).
*
* @author Martin Davis
*
*/
public class Matrix
{
private static void swapRows(double[][] m, int i, int j)
{
if (i == j) return;
for (int col = 0; col < m[0].length; col++) {
double temp = m[i][col];
m[i][col] = m[j][col];
m[j][col] = temp;
}
}
private static void swapRows(double[] m, int i, int j)
{
if (i == j) return;
double temp = m[i];
m[i] = m[j];
m[j] = temp;
}
/**
* Solves a system of equations using Gaussian Elimination.
* In order to avoid overhead the algorithm runs in-place
* on A - if A should not be modified the client must supply a copy.
*
* @param a an nxn matrix in row/column order )modified by this method)
* @param b a vector of length n
*
* @return a vector containing the solution (if any)
* or null if the system has no or no unique solution
*
* @throws IllegalArgumentException if the matrix is the wrong size
*/
public static double[] solve( double[][] a, double[] b )
{
int n = b.length;
if ( a.length != n || a[0].length != n )
throw new IllegalArgumentException("Matrix A is incorrectly sized");
// Use Gaussian Elimination with partial pivoting.
// Iterate over each row
for (int i = 0; i < n; i++ ) {
// Find the largest pivot in the rows below the current one.
int maxElementRow = i;
for (int j = i + 1; j < n; j++ )
if ( Math.abs( a[j][i] ) > Math.abs( a[maxElementRow][i] ) )
maxElementRow = j;
if ( a[maxElementRow][i] == 0.0 )
return null;
// Exchange current row and maxElementRow in A and b.
swapRows(a, i, maxElementRow );
swapRows(b, i, maxElementRow );
// Eliminate using row i
for (int j = i + 1; j < n; j++ ) {
double rowFactor = a[j][i] / a[i][i];
for (int k = n - 1; k >= i; k-- )
a[j][k] -= a[i][k] * rowFactor;
b[j] -= b[i] * rowFactor;
}
}
/**
* A is now (virtually) in upper-triangular form.
* The solution vector is determined by back-substitution.
*/
double[] solution = new double[n];
for (int j = n - 1; j >= 0; j-- ) {
double t = 0.0;
for (int k = j + 1; k < n; k++ )
t += a[j][k] * solution[k];
solution[j] = ( b[j] - t ) / a[j][j];
}
return solution;
}
}