Merge pull request #30 from Uditgulati/linalg

Implement NumRuby::Linalg

.gitignore

@@ -3,4 +3,6 @@
 Gemfile.lock
 /.yardoc
 .rake_tasks~
-*.so
+*.so
+.vscode
+

lib/nmatrix/lapack.rb

@@ -1,8 +1,20 @@
 module NumRuby::Linalg
-  def self.inv(obj)
-    if obj.is_a?(NMatrix)
-      return obj.invert
+  def self.inv(matrix)
+    if not matrix.is_a?(NMatrix)
+      raise("Invalid matrix. Not of type NMatrix.")
+    end
+    if matrix.dim != 2
+      raise("Invalid shape of matrix. Should be 2.")
+    end
+    if matrix.shape[0] != matrix.shape[1]
+      raise("Invalid shape. Expected square matrix.")
     end
+    m, n = matrix.shape
+
+    lu, ipiv = NumRuby::Lapack.getrf(matrix)
+    inv_a = NumRuby::Lapack.getri(lu, ipiv)
+
+    return inv_a
   end
 
   def self.dot(lha, rha)
@@ -13,12 +25,22 @@ def self.norm
 
   end
 
-  def self.solve
-
+  def self.solve(a, b, sym_pos: False, lower: False, assume_a: "gen", transposed: False)
+    # TODO: implement this and remove NMatrix.solve
   end
 
-  def self.det
+  def self.det(matrix)
+    if not matrix.is_a?(NMatrix)
+      raise("Invalid matrix. Not of type NMatrix.")
+    end
+    if matrix.dim != 2
+      raise("Invalid shape of matrix. Should be 2.")
+    end
+    if matrix.shape[0] != matrix.shape[1]
+      raise("Invalid shape. Expected square matrix.")
+    end
 
+    return matrix.det
   end
 
   def self.least_square
@@ -48,23 +70,47 @@ def self.eigvalsh
   # Matrix Decomposition
 
 
-  def self.lu(matrix)
+  def self.lu(matrix, permute_l: False)
+    if not matrix.is_a?(NMatrix)
+      raise("Invalid matrix. Not of type NMatrix.")
+    end
+    if matrix.dim != 2
+      raise("Invalid shape of matrix. Should be 2.")
+    end
+
+    lu, ipiv = NumRuby::Linalg.getrf(matrix)
 
+    # TODO: calulate p, l, u
   end
 
   def self.lu_factor(matrix)
+    if not matrix.is_a?(NMatrix)
+      raise("Invalid matrix. Not of type NMatrix.")
+    end
+    if matrix.dim != 2
+      raise("Invalid shape of matrix. Should be 2.")
+    end
+    if matrix.shape[0] != matrix.shape[1]
+      raise("Invalid shape. Expected square matrix.")
+    end
 
+    lu, ipiv = NumRuby::Linalg.getrf(matrix)
+
+    return [lu, ipiv]
   end
 
-  def self.lu_solve(matrix, rhs_val)
+  def self.lu_solve(lu, ipiv, b, trans: 0)
+    if lu.shape[0] != b.shape[0]
+      raise("Incompatibel dimensions.")
+    end
 
+    x = NumRuby::Lapack.getrs(lu, ipiv, b, trans)
+    return x
   end
 
-  # Computes the QR decomposition of matrix.
+  # Computes the SVD decomposition of matrix.
   # Args:
   # - input matrix, type: NMatrix
-  # - mode, type: String
-  # - pivoting, type: Boolean
   def self.svd(matrix)
 
   end
@@ -89,11 +135,24 @@ def self.cholesky_solve(matrix)
 
   end
 
-  # Computes the QR decomposition of matrix.
+  # Computes QR decomposition of a matrix.
+  #
+  # Calculates the decomposition A = Q*R where Q is unitary/orthogonal and R is upper triangular.
+  #
   # Args:
-  # - input matrix, type: NMatrix
+  # - matrix, type: NMatrix
+  #     Matrix to be decomposed
   # - mode, type: String
+  #     Determines what information is to be returned: either both Q and R
+  #     ('full', default), only R ('r') or both Q and R but computed in
+  #     economy-size ('economic', see Notes). The final option 'raw'
+  #     (added in Scipy 0.11) makes the function return two matrices
+  #     (Q, TAU) in the internal format used by LAPACK.
   # - pivoting, type: Boolean
+  #     Whether or not factorization should include pivoting for rank-revealing
+  #     qr decomposition. If pivoting, compute the decomposition
+  #     A*P = Q*R as above, but where P is chosen such that the diagonal
+  #     of R is non-increasing.
   def self.qr(matrix, mode: "full", pivoting: false)
     if not ['full', 'r', 'economic', 'raw'].include?(mode.downcase)
       raise("Invalid mode. Should be one of ['full', 'r', 'economic', 'raw']")
@@ -106,9 +165,42 @@ def self.qr(matrix, mode: "full", pivoting: false)
     end
     m, n = matrix.shape
 
-    if pivoting == false
+    if pivoting == true
+      qr, tau, jpvt = NumRuby::Lapack.geqp3(matrix)
+      jpvt -= 1
+    else
       qr, tau = NumRuby::Lapack.geqrf(matrix)
+    end
+
+    # calculate R here for both pivot true & false
+
+    if ['economic', 'raw'].include?(mode.downcase) or m < n
+      r = NumRuby.triu(matrix)
+    else
+      r = NumRuby.triu(matrix[0...n, 0...n])
+    end
+
+    if pivoting == true
+      rj = r, jpvt
+    else
+      rj = r
+    end
+
+    if mode == 'r'
+      return rj
+    elsif mode == 'raw'
       return [qr, tau]
     end
+
+    if m < n
+      q = NumRuby::Lapack.orgqr(qr[0...m, 0...m], tau)
+    elsif mode == 'economic'
+      q = NumRuby::Lapack.orgqr(qr, tau)
+    else
+      # TODO: Implement slice view and set slice
+      q = NumRuby::Lapack.orgqr(qr, tau)
+    end
+
+    return q, rj
   end
 end