Permalink
Browse files

Determinant can now be calculated for basically any matrix.

  • Loading branch information...
1 parent 5727559 commit 770b464bb36a59887f3082d8e3e3889ac26c4a7b @mohawkjohn committed Nov 1, 2012
Showing with 50 additions and 5 deletions.
  1. +0 −1 History.txt
  2. +45 −4 lib/nmatrix/nmatrix.rb
  3. +5 −0 spec/nmatrix_spec.rb
View
@@ -61,4 +61,3 @@
* Yale-to-list casting
* Now requires packable-1.3.5 or higher, fixing a problem with MATLAB .mat v5 file I/O (specific to doubles)
-
View
@@ -84,22 +84,63 @@ def pretty_print(q = nil)
alias :pp :pretty_print
- # Use LAPACK to calculate the inverse of the matrix (in-place).
+ # Use LAPACK to calculate the inverse of the matrix (in-place). Only works on dense matrices.
#
# Note: If you don't have LAPACK, e.g., on a Mac, this may not work yet.
def invert!
# Get the pivot array; factor the matrix
- pivot = NMatrix::LAPACK::clapack_getrf(:row, self.shape[0], self.shape[1], self, self.shape[0])
+ pivot = self.getrf!
# Now calculate the inverse using the pivot array
NMatrix::LAPACK::clapack_getri(:row, self.shape[0], self, self.shape[0], pivot)
self
end
- # Make a copy of the matrix, then invert it (requires LAPACK)
+ # Make a copy of the matrix, then invert it (requires LAPACK). Returns a dense matrix.
def invert
- self.clone.invert!
+ self.cast(:dense, self.dtype).invert!
+ end
+
+ alias :inverse :invert
+
+ # Calls clapack_getrf and returns the pivot array (dense only).
+ def getrf!
+ raise(StorageTypeError, "ATLAS functions only work on dense matrices") unless self.stype == :dense
+ NMatrix::LAPACK::clapack_getrf(:row, self.shape[0], self.shape[1], self, self.shape[0])
+ end
+
+ # Calculate the determinant by way of LU decomposition. This is accomplished using
+ # clapack_getrf, and then by summing the diagonal elements. There is a risk of
+ # underflow/overflow.
+ #
+ # There are probably also more efficient ways to calculate the determinant. This method
+ # requires making a copy of the matrix, since clapack_getrf modifies its input.
+ #
+ # For smaller matrices, you may be able to use det_exact.
+ #
+ # This function is guaranteed to return the same type of data in the matrix upon which it is called.
+ # In other words, if you call it on a rational matrix, you'll get a rational number back.
+ #
+ # Integer matrices are converted to rational matrices for the purposes of performing the calculation,
+ # as xGETRF can't work on integer matrices.
+ def det
+ raise(NotImplementedError, "determinant can be calculated only for 2D matrices") unless self.dim == 2
+
+ # Cast to a dtype for which getrf is implemented
+ new_dtype = [:byte,:int8,:int16,:int32,:int64].include?(self.dtype) ? :rational128 : self.dtype
+ copy = self.cast(:dense, new_dtype)
+
+ # Need to know the number of permutations. We'll add up the diagonals of the factorized matrix.
+ pivot = copy.getrf!
+
+ prod = pivot.size % 2 == 1 ? -1 : 1 # odd permutations => negative
+ [shape[0],shape[1]].min.times do |i|
+ prod *= copy[i,i]
+ end
+
+ # Convert back to an integer if necessary
+ new_dtype != self.dtype ? prod.to_i : prod
end
View
@@ -35,6 +35,11 @@
x = a.det_exact
end
+ it "calculates determinants" do
+ m = NMatrix.new(:dense, 3, [-2,2,3,-1,1,3,2,0,-1])
+ m.det.should == 6
+ end
+
it "allows stype casting of a dim 2 matrix between dense, sparse, and list (different dtypes)" do
m = NMatrix.new(:dense, [3,3], [0,0,1,0,2,0,3,4,5], :int64).
cast(:yale, :int32).

0 comments on commit 770b464

Please sign in to comment.