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

   /** 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 18547 days, 19 hours, 14 minutes ago.