Fix inverse_exact

ext/nmatrix/math/math.h

@@ -102,8 +102,14 @@ extern "C" {
   void nm_math_solve(VALUE lu, VALUE b, VALUE x, VALUE ipiv);
   void nm_math_inverse(const int M, void* A_elements, nm::dtype_t dtype);
   void nm_math_hessenberg(VALUE a);
-  void nm_math_det_exact(const int M, const void* elements, const int lda, nm::dtype_t dtype, void* result);
-  void nm_math_inverse_exact(const int M, const void* A_elements, const int lda, void* B_elements, const int ldb, nm::dtype_t dtype);
+  void nm_math_det_exact_from_dense(const int M, const void* elements, 
+      const int lda, nm::dtype_t dtype, void* result);
+  void nm_math_det_exact_from_yale(const int M, const YALE_STORAGE* storage, 
+      const int lda, nm::dtype_t dtype, void* result);
+  void nm_math_inverse_exact_from_dense(const int M, const void* A_elements, 
+      const int lda, void* B_elements, const int ldb, nm::dtype_t dtype);
+  void nm_math_inverse_exact_from_yale(const int M, const YALE_STORAGE* storage, 
+      const int lda, YALE_STORAGE* inverse, const int ldb, nm::dtype_t dtype);
 }
 
 

ext/nmatrix/ruby_nmatrix.c

@@ -1456,7 +1456,8 @@ static VALUE nm_init(int argc, VALUE* argv, VALUE nm) {
 
 
 /*
- * Helper for nm_cast which uses the C types instead of the Ruby objects. Called by nm_cast.
+ * Helper for nm_cast_with_types which uses the C types instead of the Ruby objects. 
+ * Called by nm_cast_with_types.
  */
 NMATRIX* nm_cast_with_ctype_args(NMATRIX* self, nm::stype_t new_stype, nm::dtype_t new_dtype, void* init_ptr) {
 
@@ -1474,6 +1475,23 @@ NMATRIX* nm_cast_with_ctype_args(NMATRIX* self, nm::stype_t new_stype, nm::dtype
   return lhs;
 }
 
+/*
+ * Cast NMatrix with given new_stype and new_dtype. Called by nm_cast.
+ */
+VALUE nm_cast_with_types(VALUE self, nm::stype_t new_stype, nm::dtype_t new_dtype, 
+        void* init_ptr) {
+  NMATRIX *rhs;
+
+  UnwrapNMatrix( self, rhs );
+
+  NMATRIX* m = nm_cast_with_ctype_args(rhs, new_stype, new_dtype, init_ptr);
+  nm_register_nmatrix(m);
+
+  VALUE to_return = Data_Wrap_Struct(CLASS_OF(self), nm_mark, nm_delete, m);
+
+  nm_unregister_nmatrix(m);
+  return to_return;
+}
 
 /*
  * call-seq:
@@ -1490,19 +1508,11 @@ VALUE nm_cast(VALUE self, VALUE new_stype_symbol, VALUE new_dtype_symbol, VALUE
   nm::stype_t new_stype = nm_stype_from_rbsymbol(new_stype_symbol);
 
   CheckNMatrixType(self);
-  NMATRIX *rhs;
-
-  UnwrapNMatrix( self, rhs );
-
   void* init_ptr = NM_ALLOCA_N(char, DTYPE_SIZES[new_dtype]);
   rubyval_to_cval(init, new_dtype, init_ptr);
 
-  NMATRIX* m = nm_cast_with_ctype_args(rhs, new_stype, new_dtype, init_ptr);
-  nm_register_nmatrix(m);
-
-  VALUE to_return = Data_Wrap_Struct(CLASS_OF(self), nm_mark, nm_delete, m);
+  VALUE to_return = nm_cast_with_types(self, new_stype, new_dtype, init_ptr);
 
-  nm_unregister_nmatrix(m);
   NM_CONSERVATIVE(nm_unregister_value(&self));
   NM_CONSERVATIVE(nm_unregister_value(&init));
   return to_return;
@@ -2946,21 +2956,38 @@ static VALUE nm_inverse(VALUE self, VALUE inverse, VALUE bang) {
  * Does not test for invertibility!
  */
 static VALUE nm_inverse_exact(VALUE self, VALUE inverse, VALUE lda, VALUE ldb) {
-
-  if (NM_STYPE(self) != nm::DENSE_STORE) {
-    rb_raise(rb_eNotImpError, "needs exact determinant implementation for this matrix stype");
-    return Qnil;
-  }
-
   if (NM_DIM(self) != 2 || NM_SHAPE0(self) != NM_SHAPE1(self)) {
     rb_raise(nm_eShapeError, "matrices must be square to have an inverse defined");
     return Qnil;
   }
 
-  nm_math_inverse_exact(NM_SHAPE0(self), 
-    NM_STORAGE_DENSE(self)->elements, FIX2INT(lda), 
-    NM_STORAGE_DENSE(inverse)->elements, FIX2INT(ldb), NM_DTYPE(self));
+  nm::dtype_t dtype = NM_DTYPE(self);
+  void* result = NM_ALLOCA_N(char, DTYPE_SIZES[dtype]);
+  if (dtype == nm::RUBYOBJ) {
+    nm_register_values(reinterpret_cast<VALUE*>(result), 1);
+  }
+  nm::stype_t old_stype = NM_STYPE(self);
+  if (old_stype == nm::LIST_STORE) {
+    self = nm_cast_with_types(self, nm::YALE_STORE, dtype, result);
+    inverse = nm_cast_with_types(inverse, nm::YALE_STORE, dtype, result);
+  }
+
+  if (NM_STYPE(self) == nm::DENSE_STORE) {
+    nm_math_inverse_exact_from_dense(NM_SHAPE0(self), 
+      NM_STORAGE_DENSE(self)->elements, FIX2INT(lda), 
+      NM_STORAGE_DENSE(inverse)->elements, FIX2INT(ldb), dtype);
+  } else {
+    nm_math_inverse_exact_from_yale(NM_SHAPE0(self), 
+      NM_STORAGE_YALE(self), FIX2INT(lda), 
+      NM_STORAGE_YALE(inverse), FIX2INT(ldb), dtype);
+  }
 
+  if (old_stype == nm::LIST_STORE) {
+    inverse = nm_cast_with_types(inverse, nm::LIST_STORE, dtype, result);
+  }
+  if (dtype == nm::RUBYOBJ) {
+    nm_unregister_values(reinterpret_cast<VALUE*>(result), 1);
+  }
   return inverse;
 }
 
@@ -2973,21 +3000,27 @@ static VALUE nm_inverse_exact(VALUE self, VALUE inverse, VALUE lda, VALUE ldb) {
  */
 static VALUE nm_det_exact(VALUE self) {
 
-  if (NM_STYPE(self) != nm::DENSE_STORE) {
-    rb_raise(rb_eNotImpError, "can only calculate exact determinant for dense matrices");
-    return Qnil;
-  }
   if (NM_DIM(self) != 2 || NM_SHAPE0(self) != NM_SHAPE1(self)) {
     rb_raise(nm_eShapeError, "matrices must be square to have a determinant defined");
     return Qnil;
   }
 
+  nm::dtype_t dtype = NM_DTYPE(self);
+  void* result = NM_ALLOCA_N(char, DTYPE_SIZES[dtype]);
+  if (NM_STYPE(self) == nm::LIST_STORE) {
+    self = nm_cast_with_types(self, nm::YALE_STORE, dtype, result);
+  }
+
   NM_CONSERVATIVE(nm_register_value(&self));
 
   // Calculate the determinant and then assign it to the return value
-  void* result = NM_ALLOCA_N(char, DTYPE_SIZES[NM_DTYPE(self)]);
-  nm::dtype_t dtype = NM_DTYPE(self);
-  nm_math_det_exact(NM_SHAPE0(self), NM_STORAGE_DENSE(self)->elements, NM_SHAPE0(self), NM_DTYPE(self), result);
+  if (NM_STYPE(self) == nm::DENSE_STORE) {
+    nm_math_det_exact_from_dense(NM_SHAPE0(self), NM_STORAGE_DENSE(self)->elements, 
+          NM_SHAPE0(self), NM_DTYPE(self), result);
+  } else {
+    nm_math_det_exact_from_yale(NM_SHAPE0(self), NM_STORAGE_YALE(self), 
+          NM_SHAPE0(self), NM_DTYPE(self), result);
+  }
 
   VALUE to_return;
   if (dtype == nm::RUBYOBJ) {

lib/nmatrix/math.rb

@@ -112,6 +112,62 @@ def invert
   end
   alias :inverse :invert
 
+  #
+  # call-seq:
+  #     invert_exact! -> NMatrix
+  #
+  # Calulates inverse_exact of a matrix of size 2 or 3.
+  # Only works on dense matrices.
+  #
+  # * *Raises* :
+  #   - +StorageTypeError+ -> only implemented on dense matrices.
+  #   - +ShapeError+ -> matrix must be square.
+  #   - +DataTypeError+ -> cannot invert an integer matrix in-place.
+  #   - +NotImplementedError+ -> cannot find exact inverse of matrix with size greater than 3
+  #
+  def invert_exact!
+    raise(StorageTypeError, "invert only works on dense matrices currently") unless self.dense?
+    raise(ShapeError, "Cannot invert non-square matrix") unless self.dim == 2 && self.shape[0] == self.shape[1]
+    raise(DataTypeError, "Cannot invert an integer matrix in-place") if self.integer_dtype?
+    #No internal implementation of getri, so use this other function
+    n = self.shape[0]
+    if n>3
+      raise(NotImplementedError, "Cannot find exact inverse of matrix of size greater than 3")
+    else
+      clond=self.clone
+      __inverse_exact__(clond, n, n)
+    end
+  end
+
+  #
+  # call-seq:
+  #     invert_exact -> NMatrix
+  #
+  # Make a copy of the matrix, then invert using exact_inverse
+  #
+  # * *Returns* :
+  #   - A dense NMatrix. Will be the same type as the input NMatrix,
+  #   except if the input is an integral dtype, in which case it will be a
+  #   :float64 NMatrix.
+  #
+  # * *Raises* :
+  #   - +StorageTypeError+ -> only implemented on dense matrices.
+  #   - +ShapeError+ -> matrix must be square.
+  #   - +NotImplementedError+ -> cannot find exact inverse of matrix with size greater than 3
+  #
+  def invert_exact
+    #write this in terms of invert_exact! so plugins will only have to overwrite
+    #invert_exact! and not invert_exact
+    if self.integer_dtype?
+      cloned = self.cast(dtype: :float64)
+      cloned.invert_exact!
+    else
+      cloned = self.clone
+      cloned.invert_exact!
+    end
+  end
+  alias :inverse_exact :invert_exact
+
   #
   # call-seq:
   #     adjugate! -> NMatrix
@@ -1393,4 +1449,4 @@ def dtype_for_floor_or_ceil
       self.__dense_map__ { |l| l.send(op,rhs) }
     end
   end
-end
+end

lib/nmatrix/shortcuts.rb

@@ -536,7 +536,12 @@ def block_diagonal(*params)
     def random(shape, opts={})
       scale = opts.delete(:scale) || 1.0
 
-      rng = Random.new
+      if opts[:seed].nil?
+        rng = Random.new
+      else
+        rng = Random.new(opts[:seed])
+      end
+
 
       random_values = []
 

spec/00_nmatrix_spec.rb

@@ -41,8 +41,71 @@
     expect { n[0] }.to raise_error(ArgumentError)
   end
 
-  it "calculates exact determinants on small square matrices" do
+  it "calculates exact determinants on small dense matrices" do
     expect(NMatrix.new(2, [1,2,3,4], stype: :dense, dtype: :int64).det_exact).to eq(-2)
+    expect(NMatrix.new(3, [1,2,3,0,5,6,7,8,0], stype: :dense, dtype: :int64)
+           .det_exact).to eq(-69)
+  end
+
+  it "calculates exact determinants on small yale square matrices" do
+    expect(NMatrix.new(2, [1,2,3,4], stype: :yale, dtype: :int64).det_exact).to eq(-2)
+    expect(NMatrix.new(3, [1,2,3,0,5,6,7,8,0], stype: :yale, dtype: :int64)
+           .det_exact).to eq(-69)
+  end
+
+  it "calculates exact determinants on small list square matrices" do
+    expect(NMatrix.new(2, [1,2,3,4], stype: :list, dtype: :int64).det_exact).to eq(-2)
+    expect(NMatrix.new(3, [1,2,3,0,5,6,7,8,0], stype: :list, dtype: :int64)
+           .det_exact).to eq(-69)
+  end
+
+  it "calculates inverse exact determinants on small dense matrices" do
+    pending("not yet implemented for NMatrix-JRuby") if jruby?
+    a = NMatrix.new(3, [1,2,3,0,1,4,5,6,0], stype: :dense, dtype: :int64)
+    inversed = a.method(:__inverse_exact__).call(a.clone, 3, 3)
+    b = NMatrix.new(3, [-24,18,5,20,-15,-4,-5,4,1], stype: :dense, dtype: :int64)
+    expect(inversed).to eq(b)
+
+    c = NMatrix.new(3, [1,0,3,0,0,1,0,6,0], stype: :dense, dtype: :int64)
+    inversed = c.method(:__inverse_exact__).call(c.clone, 3, 3)
+    d = NMatrix.new(3, [1,-3,0,0,0,0,0,1,0], stype: :dense, dtype: :int64)
+    expect(inversed).to eq(d)
+
+    e = NMatrix.new(2, [3,1,2,1], stype: :dense, dtype: :int64)
+    inversed = e.method(:__inverse_exact__).call(e.clone, 2, 2)
+    f = NMatrix.new(2, [1,-1,-2,-3], stype: :dense, dtype: :int64)
+    expect(inversed).to eq(f)
+  end
+
+  it "calculates inverse exact determinants on small yale matrices" do
+    pending("not yet implemented for NMatrix-JRuby") if jruby?
+    a = NMatrix.new(3, [1,2,3,0,1,4,5,6,0], stype: :yale, dtype: :int64)
+    inversed = a.method(:__inverse_exact__).call(a.clone, 3, 3)
+    b = NMatrix.new(3, [-24,18,5,20,-15,-4,-5,4,1], stype: :yale, dtype: :int64)
+    expect(inversed).to eq(b)
+
+    c = NMatrix.new(3, [1,0,3,0,0,1,0,6,0], stype: :yale, dtype: :int64)
+    inversed = c.method(:__inverse_exact__).call(c.clone, 3, 3)
+    d = NMatrix.new(3, [1,-3,0,0,0,0,0,1,0], stype: :yale, dtype: :int64)
+    expect(inversed).to eq(d)
+
+    e = NMatrix.new(2, [3,1,2,1], stype: :yale, dtype: :int64)
+    inversed = e.method(:__inverse_exact__).call(e.clone, 2, 2)
+    f = NMatrix.new(2, [1,-1,-2,-3], stype: :yale, dtype: :int64)
+    expect(inversed).to eq(f)
+  end
+
+  it "calculates inverse exact determinants on small list matrices" do
+    pending("not yet implemented for NMatrix-JRuby") if jruby?
+    a = NMatrix.new(3, [1,2,3,0,1,4,5,6,0], stype: :list, dtype: :int64)
+    inversed = a.method(:__inverse_exact__).call(a.clone, 3, 3)
+    b = NMatrix.new(3, [-24,18,5,20,-15,-4,-5,4,1], stype: :list, dtype: :int64)
+    expect(inversed).to eq(b)
+
+    c = NMatrix.new(2, [3,1,2,1], stype: :list, dtype: :int64)
+    inversed = c.method(:__inverse_exact__).call(c.clone, 2, 2)
+    d = NMatrix.new(2, [1,-1,-2,-3], stype: :list, dtype: :int64)
+    expect(inversed).to eq(d)
   end
 
   it "calculates determinants" do

spec/math_spec.rb

@@ -477,7 +477,7 @@
         end
       end
 
-      it "should correctly invert a matrix out-of-place" do
+      it "should correctly invert a dense matrix out-of-place" do
         a = NMatrix.new(:dense, 3, [1,2,3,0,1,4,5,6,0], dtype)
 
         if a.integer_dtype?
@@ -488,6 +488,23 @@
 
         expect(a.invert).to be_within(err).of(b)
       end
+
+      it "should correctly find exact inverse" do
+        pending("not yet implemented for NMatrix-JRuby") if jruby?
+        a = NMatrix.new(:dense, 3, [1,2,3,0,1,4,5,6,0], dtype)
+        b = NMatrix.new(:dense, 3, [-24,18,5,20,-15,-4,-5,4,1], dtype)
+
+        expect(a.invert_exact).to be_within(err).of(b)
+      end
+
+      it "should correctly find exact inverse" do
+        pending("not yet implemented for NMatrix-JRuby") if jruby?
+        a = NMatrix.new(:dense, 2, [1,3,3,8,], dtype)
+        b = NMatrix.new(:dense, 2, [-8,3,3,-1], dtype)
+
+        expect(a.invert_exact).to be_within(err).of(b)
+      end
+
     end
   end
 

spec/shortcuts_spec.rb

@@ -137,6 +137,17 @@
       expect(m.dtype).to eq(:float64)
     end
 
+    it "creates a matrix of random numbers with defined seed value" do
+      m1 = NMatrix.random(2,:seed => 62)
+      m2 = NMatrix.random(2,:seed => 62)
+      m3 = NMatrix.random(2,:seed => 65)
+
+
+      expect(m1).to eq(m2)
+      expect(m1).not_to eq(m3)
+
+    end
+
     it "creates a complex matrix of random numbers" do
       m = NMatrix.random(2, :dtype => :complex128)
     end