Implemented basic linear solvers.

Author: adam-rumpf

GameMakerScripts.yyp

@@ -25,8 +25,10 @@ Algorithms from computational mathematics (such as regression and interpolation)
 
 Common linear algebra algorithms for dealing with matrices and vectors. In all functions, _vectors_ are considered to be 1D arrays while _matrices_ are 2D arrays.
 
+* [`_linear_solve`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_linear_solve/_linear_solve.gml): Solves a linear system given a coefficient matrix and righthand vector.
 * [`_matrix_multiply`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_matrix_multiply/_matrix_multiply.gml): Performs matrix-matrix, matrix-vector, and vector-vector multiplication.
 * [`_matrix_transpose`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_matrix_transpose/_matrix_transpose.gml): Transposes a 2D matrix.
+* [`_triangular_solve`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_triangular_solve/_triangular_solve.gml): Solves an upper or lower triangular system.
 * [`_unit_vector`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_unit_vector/_unit_vector.gml): Defines a unit direction vector between two coordinates.
 * [`_vector_angle`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_vector_angle/_vector_angle.gml): Calculates the acute angle between two vectors.
 * [`_vector_distance`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_vector_distance/_vector_distance.gml): Calculates the Lp-distance between two vectors.

README.md

@@ -0,0 +1,42 @@
+/// @func _linear_solve(a, b)
+/// @desc Solves a linear system of the form Ax = b (via Gaussian elimination) for a given coefficient matrix and constant vector.
+/// @requires _triangular_solve
+/// @param {real[]} a Square coefficient matrix (as a 2D array).
+/// @param {real[]} b Constraint vector (as a 1D array whose length matches a).
+/// @return {real[]} Solution vector x of the linear system Ax = b (undefined in case of inconsistent dimensions or singular matrix).
+
+function _linear_solve(a, b)
+{
+	// Verify that the inputs have compatible dimensions
+	if ((array_length(a) != array_length(a[0])) || (array_length(a) != array_length(b)))
+		return undefined;
+	
+	// Initialize system
+	var n, l, u;
+	n = array_length(b); // size of system
+	l = array_create(n, array_create(n, 0)); // lower triangular matrix
+	for (var i = 0; i < n; i++)
+		l[i][i] = 1; // initialize as identity matrix
+	u = a; // upper triangular matrix
+	
+	// Perform Gaussian elimination to obtain LU decomposition
+	for (var i = 0; i < n-1; i++)
+	{
+		for (var j = i+1; j < n; j++)
+		{
+			// Verify that we are not dividing by zero
+			if (u[i][i] == 0)
+				return undefined;
+			
+			l[j][i] = u[j][i]/u[i][i];
+			for (var k = i; k < n; k++)
+				u[j][k] = u[j][k] - l[j][i]*u[i][k];
+		}
+	}
+	
+	// Solve Ly = b
+	var yy = _triangular_solve(l, b, false);
+	
+	// Solve Ux = y
+	return _triangular_solve(u, yy, true);
+}

scripts/_linear_solve/_linear_solve.gml

@@ -0,0 +1,53 @@
+/// @func _triangular_solve(a, b[, upper])
+/// @desc Solves a triangular system of the form Ax = b for a given triangular coefficient matrix and constant vector.
+/// @param {real[]} a Triangular, square coefficient matrix (as a 2D array).
+/// @param {real[]} b Constraint vector (as a 1D array whose length matches a).
+/// @param {bool} [upper=true] True if a is upper triangular, false if a is lower triangular.
+/// @return {real[]} Solution vector x of the triangular system Ax = b (undefined in case of inconsistent dimensions or singular matrix).
+
+function _triangular_solve(a, b)
+{
+	// Check for optional upper argument
+	var upper = (argument_count > 2 ? argument[2] : true);
+	
+	// Verify that the inputs have compatible dimensions
+	if ((array_length(a) != array_length(a[0])) || (array_length(a) != array_length(b)))
+		return undefined;
+	
+	// Initialize system
+	var n, xx;
+	n = array_length(b); // size of system
+	xx = array_create(n); // solution vector
+	
+	// Verify that no diagonal element is zero
+	for (var i = 0; i < n; i++)
+	{
+		if (a[i][i] == 0)
+			return undefined;
+	}
+	
+	// Separate procedure for upper or lower triangular systems
+	if (upper == true)
+	{
+		// Use back substitution for upper triangular
+		for (var i = n-1; i >= 0; i--)
+		{
+			xx[i] = b[i];
+			for (var j = i+1; j < n; j++)
+				xx[i] -= a[i][j]*xx[j];
+			xx[i] /= a[i][i];
+		}
+	}
+	else
+	{
+		// Use forward substitution for lower triangular
+		for (var i = 0; i < n; i++)
+		{
+			xx[i] = b[i];
+			for (var j = 0; j < i; j++)
+				xx[i] -= a[i][j]*xx[j];
+		}
+	}
+	
+	return xx;
+}