Doc Classic

Class Matrix

 1 /*
 2  * Copyright (c) 2016 Martin Davis.
 3  *
 4  * All rights reserved. This program and the accompanying materials
 5  * are made available under the terms of the Eclipse Public License 2.0
 6  * and Eclipse Distribution License v. 1.0 which accompanies this distribution.
 7  * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html
 8  * and the Eclipse Distribution License is available at
 9  *
10  * http://www.eclipse.org/org/documents/edl-v10.php.
11  */
12 package org.locationtech.jts.math;
13  
14 /**
15  * Implements some 2D matrix operations 
16  * (in particular, solving systems of linear equations).
17  * 
18  * @author Martin Davis
19  *
20  */
21 public class Matrix
22 {
23   private static void swapRows(double[][] m, int i, int j)
24   {
25     if (i == j) return;
26     for (int col = 0; col < m[0].length; col++) {
27       double temp = m[i][col];
28       m[i][col] = m[j][col];
29       m[j][col] = temp;
30     }
31   }
32   
33   private static void swapRows(double[] m, int i, int j)
34   {
35     if (i == j) return;
36     double temp = m[i];
37     m[i] = m[j];
38     m[j] = temp;
39   }
40   
41   /**
42    * Solves a system of equations using Gaussian Elimination.
43    * In order to avoid overhead the algorithm runs in-place
44    * on A - if A should not be modified the client must supply a copy.
45    * 
46    * @param a an nxn matrix in row/column order )modified by this method)
47    * @param b a vector of length n
48    * 
49    * @return a vector containing the solution (if any)
50    * or null if the system has no or no unique solution
51    * 
52    * @throws IllegalArgumentException if the matrix is the wrong size 
53    */
54   public static double[] solve( double[][] a, double[] b )
55   {
56     int n = b.length;
57     if ( a.length != n || a[0].length != n )
58       throw new IllegalArgumentException("Matrix A is incorrectly sized");
59     
60     // Use Gaussian Elimination with partial pivoting.
61     // Iterate over each row
62     for (int i = 0; i < n; i++ ) {
63       // Find the largest pivot in the rows below the current one.
64       int maxElementRow = i;
65       for (int j = i + 1; j < n; j++ )
66         if ( Math.abs( a[j][i] ) > Math.abs( a[maxElementRow][i] ) )
67           maxElementRow = j;
68         
69       if ( a[maxElementRow][i] == 0.0 )
70         return null;
71       
72       // Exchange current row and maxElementRow in A and b.
73       swapRows(a, i, maxElementRow );
74       swapRows(b, i, maxElementRow );
75       
76       // Eliminate using row i
77       for (int j = i + 1; j < n; j++ ) {
78         double rowFactor = a[j][i] / a[i][i];
79         for (int k = n - 1; k >= i; k-- )
80           a[j][k] -= a[i][k] * rowFactor;
81         b[j] -= b[i] * rowFactor;
82       }
83     }
84     
85     /**
86      * A is now (virtually) in upper-triangular form.
87      * The solution vector is determined by back-substitution.
88      */
89     double[] solution = new double[n];
90     for (int j = n - 1; j >= 0; j-- ) {
91       double t = 0.0;
92       for (int k = j + 1; k < n; k++ )
93         t += a[j][k] * solution[k];
94       solution[j] = ( b[j] - t ) / a[j][j];
95     }
96     return solution;
97   }
98   
99 }
100