Hopefully most recent set of changes, including each sparse with indi…

…ces for yale.

ext/nmatrix/data/data.h

@@ -615,6 +615,23 @@ namespace nm {
 {fun<nm::EW_GEQ, uint32_t, uint8_t>,fun<nm::EW_GEQ, uint32_t, int8_t>,fun<nm::EW_GEQ, uint32_t, int16_t>,fun<nm::EW_GEQ, uint32_t, int32_t>,fun<nm::EW_GEQ, uint32_t, int64_t>,fun<nm::EW_GEQ, uint32_t, float32_t>,fun<nm::EW_GEQ, uint32_t, float64_t>,fun<nm::EW_GEQ, uint32_t, nm::Complex64>,fun<nm::EW_GEQ, uint32_t, nm::Complex128>,fun<nm::EW_GEQ, uint32_t, nm::Rational32>,fun<nm::EW_GEQ, uint32_t, nm::Rational64>,fun<nm::EW_GEQ, uint32_t, nm::Rational128>,fun<nm::EW_GEQ, uint32_t, nm::RubyObject>},\
 {fun<nm::EW_GEQ, uint64_t, uint8_t>,fun<nm::EW_GEQ, uint64_t, int8_t>,fun<nm::EW_GEQ, uint64_t, int16_t>,fun<nm::EW_GEQ, uint64_t, int32_t>,fun<nm::EW_GEQ, uint64_t, int64_t>,fun<nm::EW_GEQ, uint64_t, float32_t>,fun<nm::EW_GEQ, uint64_t, float64_t>,fun<nm::EW_GEQ, uint64_t, nm::Complex64>,fun<nm::EW_GEQ, uint64_t, nm::Complex128>,fun<nm::EW_GEQ, uint64_t, nm::Rational32>,fun<nm::EW_GEQ, uint64_t, nm::Rational64>,fun<nm::EW_GEQ, uint64_t, nm::Rational128>,fun<nm::EW_GEQ, uint64_t, nm::RubyObject>}}};
 
+
+#define NAMED_BOOL_ITYPE_DTYPE_TEMPLATE_TABLE(name,  fun,  ret,  ...) \
+	static ret (*(name)[2][nm::NUM_ITYPES][nm::NUM_DTYPES])(__VA_ARGS__) = { \
+		{ \
+      {fun<false, uint8_t, uint8_t>,fun<false, uint8_t, int8_t>,fun<false, uint8_t, int16_t>,fun<false, uint8_t, int32_t>,fun<false, uint8_t, int64_t>,fun<false, uint8_t, float32_t>,fun<false, uint8_t, float64_t>,fun<false, uint8_t, nm::Complex64>,fun<false, uint8_t, nm::Complex128>,fun<false, uint8_t, nm::Rational32>,fun<false, uint8_t, nm::Rational64>,fun<false, uint8_t, nm::Rational128>,fun<false, uint8_t, nm::RubyObject>},\
+      {fun<false, uint16_t, uint8_t>,fun<false, uint16_t, int8_t>,fun<false, uint16_t, int16_t>,fun<false, uint16_t, int32_t>,fun<false, uint16_t, int64_t>,fun<false, uint16_t, float32_t>,fun<false, uint16_t, float64_t>,fun<false, uint16_t, nm::Complex64>,fun<false, uint16_t, nm::Complex128>,fun<false, uint16_t, nm::Rational32>,fun<false, uint16_t, nm::Rational64>,fun<false, uint16_t, nm::Rational128>,fun<false, uint16_t, nm::RubyObject>},\
+      {fun<false, uint32_t, uint8_t>,fun<false, uint32_t, int8_t>,fun<false, uint32_t, int16_t>,fun<false, uint32_t, int32_t>,fun<false, uint32_t, int64_t>,fun<false, uint32_t, float32_t>,fun<false, uint32_t, float64_t>,fun<false, uint32_t, nm::Complex64>,fun<false, uint32_t, nm::Complex128>,fun<false, uint32_t, nm::Rational32>,fun<false, uint32_t, nm::Rational64>,fun<false, uint32_t, nm::Rational128>,fun<false, uint32_t, nm::RubyObject>},\
+      {fun<false, uint64_t, uint8_t>,fun<false, uint64_t, int8_t>,fun<false, uint64_t, int16_t>,fun<false, uint64_t, int32_t>,fun<false, uint64_t, int64_t>,fun<false, uint64_t, float32_t>,fun<false, uint64_t, float64_t>,fun<false, uint64_t, nm::Complex64>,fun<false, uint64_t, nm::Complex128>,fun<false, uint64_t, nm::Rational32>,fun<false, uint64_t, nm::Rational64>,fun<false, uint64_t, nm::Rational128>,fun<false, uint64_t, nm::RubyObject>}  \
+    }, \
+    {  \
+      {fun<true, uint8_t, uint8_t>,fun<true, uint8_t, int8_t>,fun<true, uint8_t, int16_t>,fun<true, uint8_t, int32_t>,fun<true, uint8_t, int64_t>,fun<true, uint8_t, float32_t>,fun<true, uint8_t, float64_t>,fun<true, uint8_t, nm::Complex64>,fun<true, uint8_t, nm::Complex128>,fun<true, uint8_t, nm::Rational32>,fun<true, uint8_t, nm::Rational64>,fun<true, uint8_t, nm::Rational128>,fun<true, uint8_t, nm::RubyObject>},\
+      {fun<true, uint16_t, uint8_t>,fun<true, uint16_t, int8_t>,fun<true, uint16_t, int16_t>,fun<true, uint16_t, int32_t>,fun<true, uint16_t, int64_t>,fun<true, uint16_t, float32_t>,fun<true, uint16_t, float64_t>,fun<true, uint16_t, nm::Complex64>,fun<true, uint16_t, nm::Complex128>,fun<true, uint16_t, nm::Rational32>,fun<true, uint16_t, nm::Rational64>,fun<true, uint16_t, nm::Rational128>,fun<true, uint16_t, nm::RubyObject>},\
+      {fun<true, uint32_t, uint8_t>,fun<true, uint32_t, int8_t>,fun<true, uint32_t, int16_t>,fun<true, uint32_t, int32_t>,fun<true, uint32_t, int64_t>,fun<true, uint32_t, float32_t>,fun<true, uint32_t, float64_t>,fun<true, uint32_t, nm::Complex64>,fun<true, uint32_t, nm::Complex128>,fun<true, uint32_t, nm::Rational32>,fun<true, uint32_t, nm::Rational64>,fun<true, uint32_t, nm::Rational128>,fun<true, uint32_t, nm::RubyObject>},\
+      {fun<true, uint64_t, uint8_t>,fun<true, uint64_t, int8_t>,fun<true, uint64_t, int16_t>,fun<true, uint64_t, int32_t>,fun<true, uint64_t, int64_t>,fun<true, uint64_t, float32_t>,fun<true, uint64_t, float64_t>,fun<true, uint64_t, nm::Complex64>,fun<true, uint64_t, nm::Complex128>,fun<true, uint64_t, nm::Rational32>,fun<true, uint64_t, nm::Rational64>,fun<true, uint64_t, nm::Rational128>,fun<true, uint64_t, nm::RubyObject>} \
+    } \
+  };
+
 /*
  * Defines a static array that holds function pointers to left dtype, right
  * dtype, and itype templated versions of the specified function.

ext/nmatrix/nmatrix.cpp

@@ -428,7 +428,7 @@ static double get_time(void);
 
 void Init_nmatrix() {
 
-	#ifdef 0
+	#ifdef RDOC
 		rb_define_singleton_method(cNMatrix, "upcast", nm_upcast, 2);
 		rb_define_singleton_method(cNMatrix, "itype_by_shape", nm_itype_by_shape, 1);	
 		rb_define_method(cNMatrix, "initialize", nm_init, -1);
@@ -672,7 +672,7 @@ static void nm_delete_ref(NMATRIX* mat) {
  *
  * Get the data type (dtype) of a matrix, e.g., :byte, :int8, :int16, :int32,
  * :int64, :float32, :float64, :complex64, :complex128, :rational32,
- * :rational64, :rational128, or :object (a Ruby object).
+ * :rational64, :rational128, or :object (the last is a Ruby object).
  */
 static VALUE nm_dtype(VALUE self) {
   ID dtype = rb_intern(DTYPE_NAMES[NM_DTYPE(self)]);
@@ -723,58 +723,7 @@ static VALUE nm_upcast(VALUE self, VALUE t1, VALUE t2) {
   return ID2SYM(rb_intern( DTYPE_NAMES[ Upcast[d1][d2] ] ));
 }
 
-/*
- * Each: Yield objects directly (suitable only for a dense matrix of Ruby objects).
- */
-static VALUE nm_dense_each_direct(VALUE nm) {
-  DENSE_STORAGE* s = NM_STORAGE_DENSE(nm);
-
-  RETURN_ENUMERATOR(nm, 0, 0);
 
-  for (size_t i = 0; i < nm_storage_count_max_elements(s); ++i)
-    rb_yield( reinterpret_cast<VALUE*>(s->elements)[i] );
-
-  return nm;
-}
-
-/*
- * Each: Copy matrix elements into Ruby VALUEs before operating on them (suitable for a dense matrix).
- */
-static VALUE nm_dense_each_indirect(VALUE nm) {
-  DENSE_STORAGE* s = NM_STORAGE_DENSE(nm);
-
-  RETURN_ENUMERATOR(nm, 0, 0);
-
-  for (size_t i = 0; i < nm_storage_count_max_elements(s); ++i) {
-    VALUE v = rubyobj_from_cval((char*)(s->elements) + i*DTYPE_SIZES[NM_DTYPE(nm)], NM_DTYPE(nm)).rval;
-    rb_yield( v ); // yield to the copy we made
-  }
-
-  return nm;
-}
-
-/*
- * Borrowed this function from NArray. Handles 'each' iteration on a dense
- * matrix.
- *
- * Additionally, handles separately matrices containing VALUEs and matrices
- * containing other types of data.
- */
-static VALUE nm_dense_each(VALUE nmatrix) {
-  volatile VALUE nm = nmatrix; // Not sure this actually does anything.
-
-  if (NM_DTYPE(nm) == nm::RUBYOBJ) {
-
-    // matrix of Ruby objects -- yield those objects directly
-    return nm_dense_each_direct(nm);
-
-  } else {
-
-    // We're going to copy the matrix element into a Ruby VALUE and then operate on it. This way user can't accidentally
-    // modify it and cause a seg fault.
-    return nm_dense_each_indirect(nm);
-  }
-}
 
 /*
  * call-seq:
@@ -795,6 +744,24 @@ static VALUE nm_each(VALUE nmatrix) {
   }
 }
 
+/*
+ * Iterate over the sparse entries of any matrix. For dense and yale, this iterates over non-zero
+ * entries; for list, this iterates over non-default entries. Yields dim+1 values for each entry:
+ * i, j, ..., and the entry itself.
+ */
+static VALUE nm_each_sparse_with_indices(VALUE nmatrix) {
+  volatile VALUE nm = nmatrix;
+
+  switch(NM_STYPE(nm)) {
+  case nm::YALE_STORE:
+    return nm_yale_each_sparse_with_indices(nm);
+  default:
+    rb_raise(rb_eNotImpError, "only yale matrix's each_sparse_with_indices method works right now");
+  }
+}
+
+
+
 /*
  * Equality operator. Returns a single true or false value indicating whether
  * the matrices are equivalent.

ext/nmatrix/nmatrix.h

@@ -104,6 +104,21 @@
 
 #ifdef __cplusplus /* These are the C++ versions of the macros. */
 
+  /*
+   * If no block is given, return an enumerator. This copied straight out of ruby's include/ruby/intern.h.
+   *
+   * rb_enumeratorize is located in enumerator.c.
+   *
+   *    VALUE rb_enumeratorize(VALUE obj, VALUE meth, int argc, VALUE *argv) {
+   *      return enumerator_init(enumerator_allocate(rb_cEnumerator), obj, meth, argc, argv);
+   *    }
+   */
+  #define RETURN_ENUMERATOR(obj, argc, argv) do {				            \
+    if (!rb_block_given_p())					                              \
+      return rb_enumeratorize((obj), ID2SYM(rb_frame_this_func()),  \
+              (argc), (argv));			                                \
+    } while (0)
+
   #define NM_DECL_ENUM(enum_type, name)   nm::enum_type name
   #define NM_DECL_STRUCT(type, name)      type          name;
 

ext/nmatrix/storage/dense.cpp

@@ -188,6 +188,62 @@ void nm_dense_storage_mark(void* storage_base) {
 // Accessors //
 ///////////////
 
+
+/*
+ * Each: Yield objects directly (suitable only for a dense matrix of Ruby objects).
+ */
+static VALUE nm_dense_each_direct(VALUE nm) {
+  DENSE_STORAGE* s = NM_STORAGE_DENSE(nm);
+
+  RETURN_ENUMERATOR(nm, 0, 0);
+
+  for (size_t i = 0; i < nm_storage_count_max_elements(s); ++i)
+    rb_yield( reinterpret_cast<VALUE*>(s->elements)[i] );
+
+  return nm;
+}
+
+/*
+ * Each: Copy matrix elements into Ruby VALUEs before operating on them (suitable for a dense matrix).
+ */
+static VALUE nm_dense_each_indirect(VALUE nm) {
+  DENSE_STORAGE* s = NM_STORAGE_DENSE(nm);
+
+  RETURN_ENUMERATOR(nm, 0, 0);
+
+  for (size_t i = 0; i < nm_storage_count_max_elements(s); ++i) {
+    VALUE v = rubyobj_from_cval((char*)(s->elements) + i*DTYPE_SIZES[NM_DTYPE(nm)], NM_DTYPE(nm)).rval;
+    rb_yield( v ); // yield to the copy we made
+  }
+
+  return nm;
+}
+
+
+/*
+ * Borrowed this function from NArray. Handles 'each' iteration on a dense
+ * matrix.
+ *
+ * Additionally, handles separately matrices containing VALUEs and matrices
+ * containing other types of data.
+ */
+VALUE nm_dense_each(VALUE nmatrix) {
+  volatile VALUE nm = nmatrix; // Not sure this actually does anything.
+
+  if (NM_DTYPE(nm) == nm::RUBYOBJ) {
+
+    // matrix of Ruby objects -- yield those objects directly
+    return nm_dense_each_direct(nm);
+
+  } else {
+
+    // We're going to copy the matrix element into a Ruby VALUE and then operate on it. This way user can't accidentally
+    // modify it and cause a seg fault.
+    return nm_dense_each_indirect(nm);
+  }
+}
+
+
 /*
  * Get a slice or one element, using copying.
  *

ext/nmatrix/storage/dense.h

@@ -78,6 +78,7 @@ void						nm_dense_storage_mark(void*);
 // Accessors //
 ///////////////
 
+VALUE nm_dense_each(VALUE nmatrix);
 void*	nm_dense_storage_get(STORAGE* s, SLICE* slice);
 void*	nm_dense_storage_ref(STORAGE* s, SLICE* slice);
 void	nm_dense_storage_set(STORAGE* s, SLICE* slice, void* val);

ext/nmatrix/storage/yale.cpp

@@ -293,7 +293,7 @@ size_t max_size(YALE_STORAGE* s) {
 ///////////////
 
 /*
- * Returns a slice of YALE_STORAGE object by coppy
+ * Returns a slice of YALE_STORAGE object by copy
  *
  * Slicing-related.
  */
@@ -1138,7 +1138,57 @@ static STORAGE* matrix_multiply(const STORAGE_PAIR& casted_storage, size_t* resu
   return reinterpret_cast<STORAGE*>(result);
 }
 
-}} // end of namespace nm::yale_storage.
+
+} // end of namespace nm::yale_storage
+
+
+template <bool Direct, typename IType, typename DType>
+static VALUE yale_each_sparse_with_indices(VALUE nm) {
+  YALE_STORAGE* s = NM_STORAGE_YALE(nm);
+  DType* a    = reinterpret_cast<DType*>(s->a);
+  IType* ija  = reinterpret_cast<IType*>(s->ija);
+
+  // If we don't have a block, return an enumerator.
+  RETURN_ENUMERATOR(nm, 0, 0);
+
+  // Iterate along diagonal
+  for (size_t k = 0; k < s->shape[0]; ++k) {
+    VALUE ii = LONG2NUM(k),
+          jj = LONG2NUM(k);
+    if (Direct) {
+      rb_yield_values(3, &(a[k]), ii, jj );  // yield element, i, j
+    } else {
+      VALUE v = rubyobj_from_cval(&(a[k]), NM_DTYPE(nm)).rval;
+      rb_yield_values(3, v, ii, jj );
+    }
+  }
+
+  // Iterate through non-diagonal elements, row by row
+  for (long i = 0; i < s->shape[0]; ++i) {
+    long p      = static_cast<long>( ija[i]   ),
+         next_p = static_cast<long>( ija[i+1] );
+
+    if (next_p == p) continue; // empty row
+
+    for (; p < next_p; ++p) {
+      long j = static_cast<long>(ija[p]);
+      VALUE ii = LONG2NUM(i),
+            jj = LONG2NUM(j);
+
+      if (Direct) {
+        rb_yield_values(3, &(ija[p]), ii, jj );
+      } else {
+        VALUE v = rubyobj_from_cval(&(a[p]), NM_DTYPE(nm)).rval;
+        rb_yield_values(3, v, ii, jj );
+      }
+    }
+  }
+
+  return nm;
+}
+
+
+} // end of namespace nm.
 
 ///////////////////
 // Ruby Bindings //
@@ -1149,7 +1199,7 @@ static STORAGE* matrix_multiply(const STORAGE_PAIR& casted_storage, size_t* resu
 extern "C" {
 
 void nm_init_yale_functions() {
-	#ifdef 0
+	#ifdef RDOC
 	  rb_define_method(cNMatrix_YaleFunctions, "yale_ija", nm_ija, 0);
 	  rb_define_method(cNMatrix_YaleFunctions, "yale_a", nm_a, 0);
 	  rb_define_method(cNMatrix_YaleFunctions, "yale_size", nm_size, 0);
@@ -1180,6 +1230,26 @@ void nm_init_yale_functions() {
 // C ACCESSORS //
 /////////////////
 
+
+
+VALUE nm_yale_each_sparse_with_indices(VALUE nmatrix) {
+  nm::dtype_t d = NM_DTYPE(nmatrix);
+  nm::itype_t i = NM_ITYPE(nmatrix);
+
+  NAMED_BOOL_ITYPE_DTYPE_TEMPLATE_TABLE(ttable, nm::yale_each_sparse_with_indices, VALUE, VALUE)
+
+  if (NM_DTYPE(nmatrix) == nm::RUBYOBJ) {
+
+    // matrix of Ruby objects -- yield those objects directly
+    return ttable[true][i][d](nmatrix);
+
+  } // else:
+  // Copy the matrix element into a Ruby VALUE and then operate on it.
+  return ttable[false][i][d](nmatrix);
+
+}
+
+
 /*
  * C accessor for inserting some value in a matrix (or replacing an existing cell).
  */

ext/nmatrix/storage/yale.h

@@ -97,6 +97,7 @@ extern "C" {
   // Accessors //
   ///////////////
 
+  VALUE nm_yale_each_sparse_with_indices(VALUE nmatrix);
   void* nm_yale_storage_get(STORAGE* s, SLICE* slice);
   void*	nm_yale_storage_ref(STORAGE* s, SLICE* slice);
   char  nm_yale_storage_set(STORAGE* storage, SLICE* slice, void* v);

lib/nmatrix/version.rb

@@ -23,13 +23,6 @@
 #++
 
 class NMatrix
-  module VERSION #:nodoc:
-    # Based on Rails. https://github.com/rails/rails/blob/master/version.rb
-    MAJOR = '0'
-    MINOR = '0'
-    TINY = '2'
-
-    STRING = [MAJOR, MINOR, TINY].compact.join('.')
-  end
+  VERSION = "0.0.3"
 end
 

nmatrix.gemspec

@@ -5,7 +5,7 @@ require 'nmatrix/version'
 
 Gem::Specification.new do |gem|
   gem.name = "nmatrix"
-  gem.version = NMatrix::VERSION::STRING
+  gem.version = NMatrix::VERSION
   gem.summary = "NMatrix is an experimental linear algebra library for Ruby, written mostly in C." 
   gem.description = "NMatrix is an experimental linear algebra library for Ruby, written mostly in C." 
   gem.homepage = 'http://sciruby.com'