Matrix.frink

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]]
   }
}

"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 18620 days, 6 hours, 21 minutes ago.