Gaussian elimination

In linear algebra, Gaussian elimination is an algorithm  for solving systems of linear equations, finding the rank of a matrix, and calculating the inverse of an invertible square matrix.

Gauss-Jordan Elimination is a variant of Gaussian Elimination.

package org.javaresources.math.matrix;
/**
 * Created by IntelliJ IDEA.
 * User: glen
 * Date: 2010-3-29
 * Time: 20:10:28
 *
 */
public class MatrixUtil {
    /**
     *   Ax=b
     * @param a A  matrix
     * @param b b  vector
     */
    public static void gaussJordan(double a[][], double b[]){
        int i, j, k, m;
        double temp;
        int numRows = a.length;
        int numCols = a[0].length;
        int index[][] = new int[numRows][2];
        partialPivot(a, b, index);
        for (i = 0; i < numRows; ++i) {
            temp = a[i][i];
            for (j = 0; j < numCols; ++j) {
                a[i][j] /= temp;
            }
            b[i] /= temp;
            a[i][i] = 1.0 / temp;
            for (k = 0; k < numRows; ++k) {
                if (k != i) {
                    temp = a[k][i];
                    for (j = 0; j < numCols; ++j) {
                        a[k][j] -= temp * a[i][j];
                    }
                    b[k] -= temp * b[i];
                    a[k][i] = -temp * a[i][i];
                }
            }
        }
        for (j = numCols - 1; j >= 0; --j) {
            k = index[j][0];
            m = index[j][1];
            if (k != m) {
                for (i = 0; i < numRows; ++i) {
                    temp = a[i][m];
                    a[i][m] = a[i][k];
                    a[i][k] = temp;
                }
            }
        }
        return;
    }
    /**
     *
     * @param mt  matrix
     * @param v     vector
     * @param index
     */
    private static void partialPivot(
            double mt[][], double v[], int index[][]) {
        double temp;
        double tempRow[];
        int i, j, m;
        int numRows = mt.length;
        int numCols = mt[0].length;
        double scale[] = new double[numRows];
        for (i = 0; i < numRows; ++i) {
            index[i][0] = i;
            index[i][1] = i;
            for (j = 0; j < numCols; ++j) {
                scale[i] = Math.max(scale[i], Math.abs(mt[i][j]));
            }
        }
        for (j = 0; j < numCols - 1; ++j) {
            m = j;
            for (i = j + 1; i < numRows; ++i) {
                if (Math.abs(mt[i][j]) / scale[i] >
                        Math.abs(mt[m][j]) / scale[m]) {
                    m = i;
                }
            }
            if (m != j) {
                index[j][0] = j;
                index[j][1] = m;
                tempRow = mt[j];
                mt[j] = mt[m];
                mt[m] = tempRow;
                temp = v[j];
                v[j] = v[m];
                v[m] = temp;
                temp = scale[j];
                scale[j] = scale[m];
                scale[m] = temp;
            }
        }
    }

    /**
     *    Inverse matrix
     * @param a
     */
    public static void invertMatrix(double a[][]) {
        int i, j, k, m;
        double temp;

        int numRows = a.length;
        int numCols = a[0].length;
        int index[][] = new int[numRows][2];
        partialPivot(a, new double[numRows], index);
        for (i = 0; i < numRows; ++i) {
            temp = a[i][i];
            for (j = 0; j < numCols; ++j) {
                a[i][j] /= temp;
            }
            a[i][i] = 1.0 / temp;
            for (k = 0; k < numRows; ++k) {
                if (k != i) {
                    temp = a[k][i];
                    for (j = 0; j < numCols; ++j) {
                        a[k][j] -= temp * a[i][j];
                    }
                    a[k][i] = -temp * a[i][i];
                }
            }
        }
        for (j = numCols - 1; j >= 0; --j) {
            k = index[j][0];
            m = index[j][1];
            if (k != m) {
                for (i = 0; i < numRows; ++i) {
                    temp = a[i][m];
                    a[i][m] = a[i][k];
                    a[i][k] = temp;
                }
            }
        }
    }
}