Fixes #det method for object dtypes

ext/nmatrix/data/complex.h

@@ -110,11 +110,11 @@ class Complex {
   }
 
   /*
-   * Complex inverse function -- creates a copy, but inverted.
+   * Complex inverse (reciprocal) function -- computes 1 / n.
    *
    * FIXME: Check that this doesn't duplicate functionality of NativeType / Complex<Type>
    */
-  inline Complex<Type> inverse() const {
+  inline Complex<Type> reciprocal() const {
     Complex<Type> conj = conjugate();
     Type denom = this->r * this->r + this->i * this->i;
     return Complex<Type>(conj.r / denom, conj.i / denom);

ext/nmatrix/data/ruby_object.h

@@ -54,6 +54,12 @@
  * Classes and Functions
  */
 
+extern "C" {
+  inline VALUE quo_reciprocal(const VALUE rval) {
+    return rb_funcall(INT2FIX(1), nm_rb_quo, 1, rval);
+  }
+}
+
 namespace nm {
 template<typename T, typename U>
 struct made_from_same_template : std::false_type {}; 
@@ -125,10 +131,32 @@ class RubyObject {
   inline RubyObject(const RubyObject& other) : rval(other.rval) {}
 
   /*
-   * Inverse operator.
+   * Quotient operator (typically for rational divisions)
+   */
+  inline RubyObject quo(const RubyObject& other) const {
+    return RubyObject(rb_funcall(this->rval, nm_rb_quo, 1, other.rval));
+  }
+
+  /*
+   * Numeric inverse (reciprocal) operator, usually used for matrix
+   * operations like getrf. For this to work, Fixnum.quo(this) must not
+   * return a TypeError.
    */
-  inline RubyObject inverse() const {
-    rb_raise(rb_eNotImpError, "RubyObject#inverse needs to be implemented");
+  inline RubyObject reciprocal() const {
+    int exception;
+
+    // Attempt to call 1.quo(this).
+    VALUE result = rb_protect(quo_reciprocal, this->rval, &exception);
+    if (exception) {
+      ID rb_reciprocal = rb_intern("reciprocal");
+      // quo failed, so let's see if the object has a reciprocal method.
+      if (rb_respond_to(this->rval, rb_reciprocal)) {
+	return RubyObject(rb_funcall(this->rval, rb_reciprocal, 0));
+      } else {
+	rb_raise(rb_eNoMethodError, "expected reciprocal method, since 1.quo(object) raises an error");
+      }
+    }
+    return result;
   }
 
   /*

ext/nmatrix/math.cpp

@@ -9,8 +9,8 @@
 //
 // == Copyright Information
 //
-// SciRuby is Copyright (c) 2010 - 2014, Ruby Science Foundation
-// NMatrix is Copyright (c) 2012 - 2014, John Woods and the Ruby Science Foundation
+// SciRuby is Copyright (c) 2010 - present, Ruby Science Foundation
+// NMatrix is Copyright (c) 2012 - present, John Woods and the Ruby Science Foundation
 //
 // Please see LICENSE.txt for additional copyright notices.
 //
@@ -134,6 +134,7 @@
 #include "math/cblas_enums.h"
 
 #include "data/data.h"
+#include "math/magnitude.h"
 #include "math/imax.h"
 #include "math/scal.h"
 #include "math/laswp.h"
@@ -206,7 +207,7 @@ namespace nm {
       } else if (M == 3) {
         x = A[lda+1] * A[2*lda+2] - A[lda+2] * A[2*lda+1]; // ei - fh
         y = A[lda] * A[2*lda+2] -   A[lda+2] * A[2*lda];   // fg - di
-        x = A[0]*x - A[1]*y ; // a*(ei-fh) - b*(fg-di)
+        x = A[0]*x - A[1]*y; // a*(ei-fh) - b*(fg-di)
 
         y = A[lda] * A[2*lda+1] - A[lda+1] * A[2*lda];    // dh - eg
         *result = A[2]*y + x; // c*(dh-eg) + _
@@ -237,11 +238,14 @@ namespace nm {
       int col_index[M];
 
       for (int k = 0;k < M; ++k) {
-        DType akk = std::abs( matrix[k * (M + 1)] ) ; // diagonal element
+        typename MagnitudeDType<DType>::type akk;
+        akk = magnitude( matrix[k * (M + 1)] ); // diagonal element
+
         int interchange = k;
 
         for (int row = k + 1; row < M; ++row) {
-          DType big = std::abs( matrix[M*row + k] ); // element below the temp pivot
+          typename MagnitudeDType<DType>::type big;
+          big = magnitude( matrix[M*row + k] ); // element below the temp pivot
 
           if ( big > akk ) {
             interchange = row;
@@ -748,16 +752,16 @@ static VALUE nm_cblas_nrm2(VALUE self, VALUE n, VALUE x, VALUE incx) {
 static VALUE nm_cblas_asum(VALUE self, VALUE n, VALUE x, VALUE incx) {
 
   static void (*ttable[nm::NUM_DTYPES])(const int N, const void* X, const int incX, void* sum) = {
-      nm::math::cblas_asum<uint8_t,uint8_t>,
-      nm::math::cblas_asum<int8_t,int8_t>,
-      nm::math::cblas_asum<int16_t,int16_t>,
-      nm::math::cblas_asum<int32_t,int32_t>,
-      nm::math::cblas_asum<int64_t,int64_t>,
-      nm::math::cblas_asum<float32_t,float32_t>,
-      nm::math::cblas_asum<float64_t,float64_t>,
-      nm::math::cblas_asum<float32_t,nm::Complex64>,
-      nm::math::cblas_asum<float64_t,nm::Complex128>,
-      nm::math::cblas_asum<nm::RubyObject,nm::RubyObject>
+      nm::math::cblas_asum<uint8_t>,
+      nm::math::cblas_asum<int8_t>,
+      nm::math::cblas_asum<int16_t>,
+      nm::math::cblas_asum<int32_t>,
+      nm::math::cblas_asum<int64_t>,
+      nm::math::cblas_asum<float32_t>,
+      nm::math::cblas_asum<float64_t>,
+      nm::math::cblas_asum<nm::Complex64>,
+      nm::math::cblas_asum<nm::Complex128>,
+      nm::math::cblas_asum<nm::RubyObject>
   };
 
   nm::dtype_t dtype  = NM_DTYPE(x);

ext/nmatrix/math/asum.h

@@ -9,8 +9,8 @@
 //
 // == Copyright Information
 //
-// SciRuby is Copyright (c) 2010 - 2014, Ruby Science Foundation
-// NMatrix is Copyright (c) 2012 - 2014, John Woods and the Ruby Science Foundation
+// SciRuby is Copyright (c) 2010 - present, Ruby Science Foundation
+// NMatrix is Copyright (c) 2012 - present, John Woods and the Ruby Science Foundation
 //
 // Please see LICENSE.txt for additional copyright notices.
 //
@@ -60,6 +60,8 @@
 #define ASUM_H
 
 
+#include "math/magnitude.h"
+
 namespace nm { namespace math {
 
 /*
@@ -73,44 +75,21 @@ namespace nm { namespace math {
  *    complex64 -> float or double
  *    complex128 -> double
  */
-template <typename ReturnDType, typename DType>
-inline ReturnDType asum(const int N, const DType* X, const int incX) {
-  ReturnDType sum = 0;
-  if ((N > 0) && (incX > 0)) {
-    for (int i = 0; i < N; ++i) {
-      sum += std::abs(X[i*incX]);
-    }
-  }
-  return sum;
-}
-
-
-template <>
-inline float asum(const int N, const Complex64* X, const int incX) {
-  float sum = 0;
-  if ((N > 0) && (incX > 0)) {
-    for (int i = 0; i < N; ++i) {
-      sum += std::abs(X[i*incX].r) + std::abs(X[i*incX].i);
-    }
-  }
-  return sum;
-}
-
-template <>
-inline double asum(const int N, const Complex128* X, const int incX) {
-  double sum = 0;
+template <typename DType, typename MDType = typename MagnitudeDType<DType>::type>
+inline MDType asum(const int N, const DType* X, const int incX) {
+  MDType sum = 0;
   if ((N > 0) && (incX > 0)) {
     for (int i = 0; i < N; ++i) {
-      sum += std::abs(X[i*incX].r) + std::abs(X[i*incX].i);
+      sum += magnitude(X[i*incX]);
     }
   }
   return sum;
 }
 
 
-template <typename ReturnDType, typename DType>
+template <typename DType, typename MDType = typename MagnitudeDType<DType>::type>
 inline void cblas_asum(const int N, const void* X, const int incX, void* sum) {
-  *reinterpret_cast<ReturnDType*>( sum ) = asum<ReturnDType, DType>( N, reinterpret_cast<const DType*>(X), incX );
+  *reinterpret_cast<MDType*>( sum ) = asum<DType,MDType>( N, reinterpret_cast<const DType*>(X), incX );
 }
 
 

ext/nmatrix/math/cblas_templates_core.h

@@ -107,9 +107,9 @@ inline void cblas_rot(const int N, void* X, const int incX, void* Y, const int i
  *    complex64 -> float or double
  *    complex128 -> double
  */
-template <typename ReturnDType, typename DType>
-inline ReturnDType asum(const int N, const DType* X, const int incX) {
-  return nm::math::asum<ReturnDType,DType>(N,X,incX);
+template <typename DType, typename MDType = typename MagnitudeDType<DType>::type>
+inline MDType asum(const int N, const DType* X, const int incX) {
+  return nm::math::asum<DType,MDType>(N,X,incX);
 }
 
 
@@ -134,9 +134,9 @@ inline double asum(const int N, const Complex128* X, const int incX) {
 }
 
 
-template <typename ReturnDType, typename DType>
+template <typename DType, typename MDType = typename MagnitudeDType<DType>::type>
 inline void cblas_asum(const int N, const void* X, const int incX, void* sum) {
-  *static_cast<ReturnDType*>( sum ) = asum<ReturnDType, DType>( N, static_cast<const DType*>(X), incX );
+  *static_cast<MDType*>( sum ) = asum<DType, MDType>( N, static_cast<const DType*>(X), incX );
 }
 
 /*
@@ -149,9 +149,9 @@ inline void cblas_asum(const int N, const void* X, const int incX, void* sum) {
  *    complex64 -> float or double
  *    complex128 -> double
  */
-template <typename ReturnDType, typename DType>
+template <typename DType, typename MDType = typename MagnitudeDType<DType>::type>
 inline ReturnDType nrm2(const int N, const DType* X, const int incX) {
-  return nm::math::nrm2<ReturnDType,DType>(N, X, incX);
+  return nm::math::nrm2<DType,MDType>(N, X, incX);
 }
 
 
@@ -175,9 +175,9 @@ inline double nrm2(const int N, const Complex128* X, const int incX) {
   return cblas_dznrm2(N, X, incX);
 }
 
-template <typename ReturnDType, typename DType>
+template <typename DType, typename MDType = typename MagnitudeDType<DType>::type>
 inline void cblas_nrm2(const int N, const void* X, const int incX, void* result) {
-  *static_cast<ReturnDType*>( result ) = nrm2<ReturnDType, DType>( N, static_cast<const DType*>(X), incX );
+  *static_cast<MDType*>( result ) = nrm2<DType, MDType>( N, static_cast<const DType*>(X), incX );
 }
 
 //imax

ext/nmatrix/math/getrf.h

@@ -9,8 +9,8 @@
 //
 // == Copyright Information
 //
-// SciRuby is Copyright (c) 2010 - 2014, Ruby Science Foundation
-// NMatrix is Copyright (c) 2012 - 2014, John Woods and the Ruby Science Foundation
+// SciRuby is Copyright (c) 2010 - present, Ruby Science Foundation
+// NMatrix is Copyright (c) 2012 - present, John Woods and the Ruby Science Foundation
 //
 // Please see LICENSE.txt for additional copyright notices.
 //
@@ -59,6 +59,7 @@
 #ifndef GETRF_H
 #define GETRF_H
 
+#include "math/reciprocal.h"
 #include "math/laswp.h"
 #include "math/math.h"
 #include "math/trsm.h"
@@ -68,14 +69,6 @@
 
 namespace nm { namespace math {
 
-/* Numeric inverse -- usually just 1 / f, but a little more complicated for complex. */
-template <typename DType>
-inline DType numeric_inverse(const DType& n) {
-  return n.inverse();
-}
-template <> inline float numeric_inverse(const float& n) { return 1 / n; }
-template <> inline double numeric_inverse(const double& n) { return 1 / n; }
-
 /*
  * Templated version of row-order and column-order getrf, derived from ATL_getrfR.c (from ATLAS 3.8.0).
  *
@@ -132,21 +125,17 @@ inline int getrf_nothrow(const int M, const int N, DType* A, const int lda, int*
       An = &(Ar[N_ul]);
 
       nm::math::laswp<DType>(N_dr, Ar, lda, 0, N_ul, ipiv, 1);
-
       nm::math::trsm<DType>(CblasRowMajor, CblasRight, CblasUpper, CblasNoTrans, CblasUnit, N_dr, N_ul, one, A, lda, Ar, lda);
       nm::math::gemm<DType>(CblasRowMajor, CblasNoTrans, CblasNoTrans, N_dr, N-N_ul, N_ul, &neg_one, Ar, lda, Ac, lda, &one, An, lda);
-
       i = getrf_nothrow<true,DType>(N_dr, N-N_ul, An, lda, ipiv+N_ul);
     } else {
       Ar = NULL;
       Ac = &(A[N_ul * lda]);
       An = &(Ac[N_ul]);
 
       nm::math::laswp<DType>(N_dr, Ac, lda, 0, N_ul, ipiv, 1);
-
       nm::math::trsm<DType>(CblasColMajor, CblasLeft, CblasLower, CblasNoTrans, CblasUnit, N_ul, N_dr, one, A, lda, Ac, lda);
       nm::math::gemm<DType>(CblasColMajor, CblasNoTrans, CblasNoTrans, M-N_ul, N_dr, N_ul, &neg_one, &(A[N_ul]), lda, Ac, lda, &one, An, lda);
-
       i = getrf_nothrow<false,DType>(M-N_ul, N_dr, An, lda, ipiv+N_ul);
     }
 
@@ -157,7 +146,6 @@ inline int getrf_nothrow(const int M, const int N, DType* A, const int lda, int*
     }
 
     nm::math::laswp<DType>(N_ul, A, lda, N_ul, MN, ipiv, 1);  /* apply pivots */
-
   } else if (MN == 1) { // there's another case for the colmajor version, but it doesn't seem to be necessary.
 
     int i;
@@ -170,13 +158,14 @@ inline int getrf_nothrow(const int M, const int N, DType* A, const int lda, int*
     DType tmp = A[i];
     if (tmp != 0) {
 
-      nm::math::scal<DType>((RowMajor ? N : M), nm::math::numeric_inverse(tmp), A, 1);
+      nm::math::scal<DType>((RowMajor ? N : M), nm::math::reciprocal(tmp), A, 1);
       A[i] = *A;
       *A   = tmp;
 
     } else ierr = 1;
 
   }
+
   return(ierr);
 }
 

ext/nmatrix/math/imax.h

@@ -9,8 +9,8 @@
 //
 // == Copyright Information
 //
-// SciRuby is Copyright (c) 2010 - 2014, Ruby Science Foundation
-// NMatrix is Copyright (c) 2012 - 2014, John Woods and the Ruby Science Foundation
+// SciRuby is Copyright (c) 2010 - present, Ruby Science Foundation
+// NMatrix is Copyright (c) 2012 - present, John Woods and the Ruby Science Foundation
 //
 // Please see LICENSE.txt for additional copyright notices.
 //
@@ -29,8 +29,11 @@
 #ifndef IMAX_H
 #define IMAX_H
 
+#include "math/magnitude.h"
+
 namespace nm { namespace math {
 
+
 template<typename DType>
 inline int imax(const int n, const DType *x, const int incx) {
 
@@ -41,28 +44,28 @@ inline int imax(const int n, const DType *x, const int incx) {
     return 0;
   }
 
-  DType dmax;
+  typename MagnitudeDType<DType>::type dmax;
   int imax = 0;
 
   if (incx == 1) { // if incrementing by 1
 
-    dmax = abs(x[0]);
+    dmax = magnitude(x[0]);
 
     for (int i = 1; i < n; ++i) {
-      if (std::abs(x[i]) > dmax) {
+      if (magnitude(x[i]) > dmax) {
         imax = i;
-        dmax = std::abs(x[i]);
+        dmax = magnitude(x[i]);
       }
     }
 
   } else { // if incrementing by more than 1
 
-    dmax = std::abs(x[0]);
+    dmax = magnitude(x[0]);
 
     for (int i = 1, ix = incx; i < n; ++i, ix += incx) {
-      if (std::abs(x[ix]) > dmax) {
+      if (magnitude(x[ix]) > dmax) {
         imax = i;
-        dmax = std::abs(x[ix]);
+        dmax = magnitude(x[ix]);
       }
     }
   }