1 /*
2 * Copyright (C) 2015 Alberto Irurueta Carro (alberto@irurueta.com)
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16 package com.irurueta.algebra;
17
18 /**
19 * Computes Gauss-Jordan elimination for provided matrix using full pivoting,
20 * which provides greater stability.
21 * Gauss-Jordan elimination can be used to compute matrix inversion or to
22 * solve linear systems of equations of the form A * x = b, where Gauss-Jordan
23 * elimination both inverts matrix A and finds solution x at the same time.
24 */
25 public class GaussJordanElimination {
26
27 /**
28 * Constructor.
29 * Prevents instantiation.
30 */
31 private GaussJordanElimination() {
32 }
33
34 /**
35 * Computes Gauss-Jordan elimination by attempting to solve linear system
36 * of equations a * x = b. This method computes inverse of matrix "a", and
37 * modifies provided matrix so that its inverse is stored in it after
38 * execution of this method. Likewise, this method modifies b so that
39 * solution x is stored on it after execution of this method.
40 * This method can only be used on squared a matrices.
41 *
42 * @param a linear system of equations matrix. Will contain its inverse
43 * after execution of this method.
44 * @param b linear system of equations parameters. Will contain the solution
45 * after execution of this method. If null is provided, solution is not
46 * stored but matrix inverse is computed anyway. Each column of b is
47 * considered a new linear system of equations and its solution x is
48 * computed on the corresponding column of b.
49 * @throws SingularMatrixException if provided matrix "a" is found to be
50 * singular.
51 * @throws WrongSizeException if provided matrix "a" is not square.
52 */
53 public static void process(final Matrix a, final Matrix b) throws SingularMatrixException, WrongSizeException {
54 if (a.getRows() != a.getColumns()) {
55 throw new WrongSizeException();
56 }
57 if (b != null && b.getRows() != a.getRows()) {
58 throw new WrongSizeException();
59 }
60
61 int i;
62 var icol = 0;
63 var irow = 0;
64 int j;
65 int k;
66 int l;
67 int ll;
68 final var n = a.getRows();
69 final var m = b != null ? b.getColumns() : 0;
70
71 double big;
72 double dum;
73 double pivinv;
74 double value;
75 final var indxr = new int[n];
76 final var indxc = new int[n];
77 final var ipiv = new int[n];
78
79 for (j = 0; j < n; j++) {
80 ipiv[j] = 0;
81 }
82 for (i = 0; i < n; i++) {
83 big = 0.0;
84 for (j = 0; j < n; j++) {
85 if (ipiv[j] != 1) {
86 for (k = 0; k < n; k++) {
87 if (ipiv[k] == 0) {
88 value = Math.abs(a.getElementAt(j, k));
89 if (value >= big) {
90 big = value;
91 irow = j;
92 icol = k;
93 }
94 }
95 }
96 }
97 }
98 ++(ipiv[icol]);
99 if (irow != icol) {
100 for (l = 0; l < n; l++) {
101 swap(a.getBuffer(), a.getBuffer(), a.getIndex(irow, l), a.getIndex(icol, l));
102 }
103 for (l = 0; l < m; l++) {
104 swap(b.getBuffer(), b.getBuffer(), b.getIndex(irow, l), b.getIndex(icol, l));
105 }
106 }
107 indxr[i] = irow;
108 indxc[i] = icol;
109 value = a.getElementAt(icol, icol);
110 if (value == 0.0) {
111 throw new SingularMatrixException();
112 }
113 pivinv = 1.0 / value;
114 a.setElementAt(icol, icol, 1.0);
115 for (l = 0; l < n; l++) {
116 a.setElementAt(icol, l, a.getElementAt(icol, l) * pivinv);
117 }
118 for (l = 0; l < m; l++) {
119 b.setElementAt(icol, l, b.getElementAt(icol, l) * pivinv);
120 }
121 for (ll = 0; ll < n; ll++) {
122 if (ll != icol) {
123 dum = a.getElementAt(ll, icol);
124 a.setElementAt(ll, icol, 0.0);
125 for (l = 0; l < n; l++) {
126 a.setElementAt(ll, l, a.getElementAt(ll, l) - a.getElementAt(icol, l) * dum);
127 }
128 for (l = 0; l < m; l++) {
129 b.setElementAt(ll, l, b.getElementAt(ll, l) - b.getElementAt(icol, l) * dum);
130 }
131 }
132 }
133 }
134 for (l = n - 1; l >= 0; l--) {
135 if (indxr[l] != indxc[l]) {
136 for (k = 0; k < n; k++) {
137 swap(a.getBuffer(), a.getBuffer(), a.getIndex(k, indxr[l]), a.getIndex(k, indxc[l]));
138 }
139 }
140 }
141 }
142
143 /**
144 * Computes Gauss-Jordan elimination by attempting to solve linear system
145 * of equations a * x = b. This method computes inverse of matrix "a", and
146 * modifies provided matrix so that its inverse is stored in it after
147 * execution of this method. Likewise, this method modifies b so that
148 * solution x is stored on it after execution of this method.
149 * This method can only be used on squared a matrices.
150 *
151 * @param a linear system of equations matrix. Will contain its inverse
152 * after execution of this method.
153 * @param b linear system of equations parameters. Will contain the solution
154 * after execution of this method. If null is provided, solution is not
155 * stored but matrix inverse is computed anyway.
156 * @throws SingularMatrixException if provided matrix "a" is found to be
157 * singular.
158 * @throws WrongSizeException if provided matrix "a" is not square.
159 */
160 public static void process(final Matrix a, final double[] b) throws SingularMatrixException, WrongSizeException {
161 final var mb = b != null ? Matrix.newFromArray(b) : null;
162 process(a, mb);
163 if (mb != null) {
164 final var buffer = mb.getBuffer();
165 System.arraycopy(buffer, 0, b, 0, b.length);
166 }
167 }
168
169 /**
170 * Computes inverse of matrix "a". No solution of a linear system of equations
171 * is computed. This method modifies provided matrix storing the inverse
172 * on it after execution of this method
173 *
174 * @param a matrix to be inverted.
175 * @throws SingularMatrixException if provided matrix is found to be
176 * singular
177 * @throws WrongSizeException if provided matrix is not square
178 */
179 public static void inverse(final Matrix a) throws SingularMatrixException, WrongSizeException {
180 process(a, (Matrix) null);
181 }
182
183 /**
184 * Swaps values in arrays at provided positions
185 *
186 * @param array1 1st array
187 * @param array2 2nd array
188 * @param pos1 1st position to be swapped
189 * @param pos2 2nd position to be swapped
190 */
191 private static void swap(final double[] array1, final double[] array2,
192 final int pos1, final int pos2) {
193 final var value1 = array1[pos1];
194 final var value2 = array2[pos2];
195
196 array1[pos1] = value2;
197 array2[pos1] = value1;
198 }
199 }