View or download Matrix.frink in plain text format
/** 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 with rows separated with newlines. */
toString[] :=
{
result = ""
for r = 0 to rows-1
result = result + " " + array@r + "\n"
return result
}
/** Returns the matrix formatted by formatTable with rows separated with
newlines. */
toTable[] :=
{
return formatTable[array]
}
/** Formats a matrix with Unicode box drawing characters. The result is a
single string. */
formatMatrix[] :=
{
return joinln[formatMatrixLines[]]
}
/** Returns the matrix formatted as a matrix with Unicode box drawing
characters. The result is an array of lines that can be further
formatted.
*/
formatMatrixLines[] :=
{
m = formatTableLines[array]
rows = length[m]
width = length[m@0]
result = new array
result.push["\u250c" + repeat[" ", width] + "\u2510"]
for n=0 to rows-1
result.push["\u2502" + m@n + "\u2502"]
result.push["\u2514" + repeat[" ", width] + "\u2518"]
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
}
/** Create a square identity matrix with the specified number of rows and
columns. The elements on the diagonal will be set to 1, the rest to
zero. This requires Frink 2020-04-19 or later. */
class makeIdentity[dimension] :=
{
return new Matrix[makeArray[[dimension, dimension], {|a,b| a==b ? 1 : 0}]]
}
/** Makes a diagonal matrix. This is passed in an array of elements (e.g.
[1,2,3] that will make up the diagonal elements. This requires Frink
2020-04-22 or later.
*/
class makeDiagonal[array] :=
{
d = length[array]
return new Matrix[makeArray[[d,d], {|a,b,data| a==b ? data@a : 0}, array]]
}
/** Returns a matrix with the specified row and column removed. Row and
column numbers are 1-indexed. This is used in many matrix operations
including calculation of determinant or of the adjugate/adjoint matrix.
*/
removeRowColumn[rowToRemove, colToRemove] :=
{
a = makeArray[[rows-1, cols-1]]
newRow = 0
ROW:
for origRow = 0 to rows-1
{
if origRow == rowToRemove-1
next ROW
newCol = 0
COL:
for origCol = 0 to rows-1
{
if origCol == colToRemove-1
next COL
a@newRow@newCol = array@origRow@origCol
newCol = newCol + 1
}
newRow = newRow + 1
}
return new Matrix[a]
}
/** Returns the adjugate or adjoint matrix. See:
https://semath.info/src/inverse-cofactor-ex4.html
*/
adjugate[] :=
{
a = makeArray[[rows,cols]]
for i = 0 to rows-1
for j = 0 to cols-1
{
a@i@j = (-1)^((i+1)+(j+1)) * removeRowColumn[j+1, i+1].det[]
if a@i@j == undef
a@i@j = 0
}
return new Matrix[a]
}
}
"class Matrix loaded OK."
View or download Matrix.frink in plain text format
This is a program written in the programming language Frink.
For more information, view the Frink
Documentation or see More Sample Frink Programs.
Alan Eliasen was born 18731 days, 1 hours, 56 minutes ago.