# Matrix.frink

``` /** 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." ```

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 18783 days, 12 hours, 2 minutes ago.