/** This is an incomplete class that implements matrix operations. Please feel free to implement missing functions! */ class Matrix { /** The array that will actually store the matrix data. Note that this is a raw 2-dimensional Frink array, and the indices are zero-based in contrast with most matrix notation which is one-based. The class should try to enforce that this is a rectangular 2-dimensional array, but the code doesn't currently enforce that with some constructors. */ var array /** The number of rows in the matrix. */ var rows /** The number of columns in the matrix. */ var cols /** Construct a new Matrix with the specified number of rows and columns and the specified default value for each element. */ new[rowCount, colCount, default=0] := { rows = rowCount cols = colCount array = makeArray[[rows, cols], default] } /** Construct a new Matrix from a rectangular 2-dimensional array. This currently does not verify that the array is rectangular. */ new[a is array] := { rows = length[a] cols = length[a@0] array = a } /** Sets the value of the specified element, using 1-based coordinates, as is common in matrix notation. Stupid mathematicians. */ set[row, col, val] := array@(row-1)@(col-1) = val /** Gets the value of the specified element, using 1-based coordinates, as is common in matrix notation. */ get[row, col] := array@(row-1)@(col-1) /** Multiply two matrices. The number of columns in column a should equal the number of rows in column b. The resulting matrix will have the same number of rows as matrix a and the same number of columns as matrix b. TODO: Use a faster algorithm for large arrays? This is O(n^3). */ class multiply[a, b] := { if a.cols != b.rows return undef else items = a.cols resRows = a.rows resCols = b.cols resultArray = makeArray[[resRows, resCols], 0] aa = a.array bb = b.array for row = 0 to resRows-1 for col = 0 to resCols-1 { sum = 0 for k = 0 to items-1 sum = sum + aa@row@k * bb@k@col resultArray@row@col = sum } return new Matrix[resultArray] } /** Multiply this matrix by the specified matrix and return the result. */ multiply[b] := multiply[this, b] /** Multiplies all elements of a matrix by a scalar. */ multiplyByScalar[s] := { a = makeArray[[rows, cols]] for rowNum = 0 to rows-1 { row = array@rowNum for colNum = 0 to cols-1 a@rowNum@colNum = row@colNum * s } return new Matrix[a] } /** Transposes the Matrix and returns the new transposed Matrix. This uses the built-in array.transpose[] method. */ transpose[] := new Matrix[array.transpose[]] /** Returns true if this Matrix is square, false otherwise. */ isSquare[] := rows == cols /** Returns a matrix as a string. */ toString[] := { result = "" for r = 0 to rows-1 result = result + " " + array@r + "\n" return result } /** Creates the LU (lower-upper) decomposition of the matrix. This creates two matrices, L and U, where L has the value 1 on the diagonal from top left to bottom right, like: 1 0 0 L = L21 1 0 L31 L32 1 U11 U12 U13 U = 0 U22 U23 0 0 U33 This is basically like solving equations by Gaussian elimination. */ LUDecompose[] := { } /** This uses Crout's method to decompose a square matrix into an lower and upper triangular matrix. This creates two matrices, L and U, where U has the value 1 on the diagonal from top left to bottom right, like: L11 0 0 L = L21 L22 0 L31 L32 L33 1 U12 U13 U = 0 1 U23 0 0 1 https://en.wikipedia.org/wiki/Crout_matrix_decomposition */ LUDecomposeCrout[] := { if ! isSquare[] { println["Matrix.LUDecomposeCrout only works on square arrays."] return undef } n = rows L = new array[[rows, cols], 0] U = new array[[rows, cols], 0] sum = 0 for i=0 to n-1 U@i@i = 1 for j=0 to n-1 { for i=j to n-1 { sum = 0 for k=0 to j-1 sum = sum + L@i@k * U@k@j L@i@j = array@i@j - sum } for i=j to n-1 { sum = 0 for k=0 to j-1 sum = sum + L@j@k * U@k@i if L@j@j == 0 { // println["Matrix.LUDecomposeCrout: det(L) is zero. Can't divide by zero."] return undef } U@j@i = (array@j@i - sum) / L@j@j } } return [new Matrix[L], new Matrix[U]] } /** Returns the determinant of a square matrix. */ det[] := { product = 1 [L, U] = LUDecomposeCrout[] // This will happen if the matrix is singular. if L == undef or U == undef return undef // Multiply diagonals of the lower triangular matrix. The // upper triangular matrix has all ones on the diagonal. // Should there be a permutation matrix here to get the signs // right? for i=0 to rows-1 product = product * (L.array)@i@i return product } } "class Matrix loaded OK."