Permalink
Browse files

Added LAPACK functions: potrf, potri, potrs, posv. Added test for potrf.

  • Loading branch information...
1 parent 6429ac7 commit 4c4a064a8367edcc961397cf512a0c2d462a9afb @mohawkjohn committed Oct 30, 2012
Showing with 355 additions and 10 deletions.
  1. +133 −0 ext/nmatrix/util/math.cpp
  2. +185 −10 ext/nmatrix/util/math.h
  3. +23 −0 lib/nmatrix/lapack.rb
  4. +14 −0 spec/lapack_spec.rb
View
@@ -133,9 +133,13 @@ extern "C" {
VALUE x, VALUE incx, VALUE vBeta, VALUE y, VALUE incy);
static VALUE nm_cblas_trsm(VALUE self, VALUE order, VALUE side, VALUE uplo, VALUE trans_a, VALUE diag, VALUE m, VALUE n,
VALUE vAlpha, VALUE a, VALUE lda, VALUE b, VALUE ldb);
+
static VALUE nm_clapack_getrf(VALUE self, VALUE order, VALUE m, VALUE n, VALUE a, VALUE lda);
+ static VALUE nm_clapack_potrf(VALUE self, VALUE order, VALUE uplo, VALUE n, VALUE a, VALUE lda);
static VALUE nm_clapack_getrs(VALUE self, VALUE order, VALUE trans, VALUE n, VALUE nrhs, VALUE a, VALUE lda, VALUE ipiv, VALUE b, VALUE ldb);
+ static VALUE nm_clapack_potrs(VALUE self, VALUE order, VALUE uplo, VALUE n, VALUE nrhs, VALUE a, VALUE lda, VALUE b, VALUE ldb);
static VALUE nm_clapack_getri(VALUE self, VALUE order, VALUE n, VALUE a, VALUE lda, VALUE ipiv);
+ static VALUE nm_clapack_potri(VALUE self, VALUE order, VALUE uplo, VALUE n, VALUE a, VALUE lda);
static VALUE nm_clapack_laswp(VALUE self, VALUE n, VALUE a, VALUE lda, VALUE k1, VALUE k2, VALUE ipiv, VALUE incx);
static VALUE nm_clapack_scal(VALUE self, VALUE n, VALUE scale, VALUE vector, VALUE incx);
@@ -251,8 +255,11 @@ void nm_math_init_blas() {
cNMatrix_LAPACK = rb_define_module_under(cNMatrix, "LAPACK");
rb_define_singleton_method(cNMatrix_LAPACK, "clapack_getrf", (METHOD)nm_clapack_getrf, 5);
+ rb_define_singleton_method(cNMatrix_LAPACK, "clapack_potrf", (METHOD)nm_clapack_potrf, 5);
rb_define_singleton_method(cNMatrix_LAPACK, "clapack_getrs", (METHOD)nm_clapack_getrs, 9);
+ rb_define_singleton_method(cNMatrix_LAPACK, "clapack_potrs", (METHOD)nm_clapack_potrs, 8);
rb_define_singleton_method(cNMatrix_LAPACK, "clapack_getri", (METHOD)nm_clapack_getri, 5);
+ rb_define_singleton_method(cNMatrix_LAPACK, "clapack_potri", (METHOD)nm_clapack_potri, 5);
rb_define_singleton_method(cNMatrix_LAPACK, "clapack_laswp", (METHOD)nm_clapack_laswp, 7);
rb_define_singleton_method(cNMatrix_LAPACK, "clapack_scal", (METHOD)nm_clapack_scal, 4);
@@ -538,6 +545,51 @@ static VALUE nm_clapack_getrf(VALUE self, VALUE order, VALUE m, VALUE n, VALUE a
}
+/* Call any of the clapack_xpotrf functions as directly as possible.
+ *
+ * You probably don't want to call this function. Instead, why don't you try clapack_potrf, which is more flexible
+ * with its arguments?
+ *
+ * This function does almost no type checking. Seriously, be really careful when you call it! There's no exception
+ * handling, so you can easily crash Ruby!
+ *
+ * Returns an array giving the pivot indices (normally these are argument #5).
+ */
+static VALUE nm_clapack_potrf(VALUE self, VALUE order, VALUE uplo, VALUE n, VALUE a, VALUE lda) {
+#ifndef HAVE_CLAPACK_H
+ rb_raise(rb_eNotImpError, "potrf currently requires LAPACK");
+#endif
+
+ static int (*ttable[nm::NUM_DTYPES])(const enum CBLAS_ORDER, const enum CBLAS_UPLO, const int n, void* a, const int lda) = {
+ NULL, NULL, NULL, NULL, NULL, // integers not allowed due to division
+ nm::math::clapack_potrf<float>,
+ nm::math::clapack_potrf<double>,
+#ifdef HAVE_CLAPACK_H
+ clapack_cpotrf, clapack_zpotrf, // call directly, same function signature!
+#else // Especially important for Mac OS, which doesn't seem to include the ATLAS clapack interface.
+ nm::math::clapack_potrf<nm::Complex64>,
+ nm::math::clapack_potrf<nm::Complex128>,
+#endif
+ NULL, NULL, NULL, NULL /*
+ nm::math::clapack_potrf<nm::Rational32>,
+ nm::math::clapack_potrf<nm::Rational64>,
+ nm::math::clapack_potrf<nm::Rational128>,
+ nm::math::clapack_potrf<nm::RubyObject> */
+ };
+
+ if (!ttable[NM_DTYPE(a)]) {
+ rb_raise(rb_eNotImpError, "this operation not yet implemented for non-BLAS dtypes");
+ // FIXME: Once BLAS dtypes are implemented, replace error above with the error below.
+ //rb_raise(nm_eDataTypeError, "this matrix operation undefined for integer matrices");
+ } else {
+ // Call either our version of potrf or the LAPACK version.
+ ttable[NM_DTYPE(a)](blas_order_sym(order), blas_uplo_sym(uplo), FIX2INT(n), NM_STORAGE_DENSE(a)->elements, FIX2INT(lda));
+ }
+
+ return a;
+}
+
+
/*
* Call any of the clapack_xgetrs functions as directly as possible.
*/
@@ -586,6 +638,42 @@ static VALUE nm_clapack_getrs(VALUE self, VALUE order, VALUE trans, VALUE n, VAL
}
+/*
+ * Call any of the clapack_xpotrs functions as directly as possible.
+ */
+static VALUE nm_clapack_potrs(VALUE self, VALUE order, VALUE uplo, VALUE n, VALUE nrhs, VALUE a, VALUE lda, VALUE b, VALUE ldb) {
+ static int (*ttable[nm::NUM_DTYPES])(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo, const int N,
+ const int NRHS, const void* A, const int lda, void* B, const int ldb) = {
+ NULL, NULL, NULL, NULL, NULL, // integers not allowed due to division
+ nm::math::clapack_potrs<float,false>,
+ nm::math::clapack_potrs<double,false>,
+#ifdef HAVE_CLAPACK_H
+ clapack_cpotrs, clapack_zpotrs, // call directly, same function signature!
+#else // Especially important for Mac OS, which doesn't seem to include the ATLAS clapack interface.
+ nm::math::clapack_potrs<nm::Complex64,true>,
+ nm::math::clapack_potrs<nm::Complex128,true>,
+#endif
+ nm::math::clapack_potrs<nm::Rational32,false>,
+ nm::math::clapack_potrs<nm::Rational64,false>,
+ nm::math::clapack_potrs<nm::Rational128,false>,
+ nm::math::clapack_potrs<nm::RubyObject,false>
+ };
+
+
+ if (!ttable[NM_DTYPE(a)]) {
+ rb_raise(nm_eDataTypeError, "this matrix operation undefined for integer matrices");
+ } else {
+
+ // Call either our version of potrs or the LAPACK version.
+ ttable[NM_DTYPE(a)](blas_order_sym(order), blas_uplo_sym(uplo), FIX2INT(n), FIX2INT(nrhs), NM_STORAGE_DENSE(a)->elements, FIX2INT(lda),
+ NM_STORAGE_DENSE(b)->elements, FIX2INT(ldb));
+ }
+
+ // b is both returned and modified directly in the argument list.
+ return b;
+}
+
+
/* Call any of the clapack_xgetri functions as directly as possible.
*
* You probably don't want to call this function. Instead, why don't you try clapack_getri, which is more flexible
@@ -643,6 +731,51 @@ static VALUE nm_clapack_getri(VALUE self, VALUE order, VALUE n, VALUE a, VALUE l
}
+/* Call any of the clapack_xpotri functions as directly as possible.
+ *
+ * You probably don't want to call this function. Instead, why don't you try clapack_potri, which is more flexible
+ * with its arguments?
+ *
+ * This function does almost no type checking. Seriously, be really careful when you call it! There's no exception
+ * handling, so you can easily crash Ruby!
+ *
+ * Returns an array giving the pivot indices (normally these are argument #5).
+ */
+static VALUE nm_clapack_potri(VALUE self, VALUE order, VALUE uplo, VALUE n, VALUE a, VALUE lda) {
+#ifndef HAVE_CLAPACK_H
+ rb_raise(rb_eNotImpError, "getri currently requires LAPACK");
+#endif
+
+ static int (*ttable[nm::NUM_DTYPES])(const enum CBLAS_ORDER, const enum CBLAS_UPLO, const int n, void* a, const int lda) = {
+ NULL, NULL, NULL, NULL, NULL, // integers not allowed due to division
+ nm::math::clapack_potri<float>,
+ nm::math::clapack_potri<double>,
+#ifdef HAVE_CLAPACK_H
+ clapack_cpotri, clapack_zpotri, // call directly, same function signature!
+#else // Especially important for Mac OS, which doesn't seem to include the ATLAS clapack interface.
+ nm::math::clapack_potri<nm::Complex64>,
+ nm::math::clapack_potri<nm::Complex128>,
+#endif
+ NULL, NULL, NULL, NULL /*
+ nm::math::clapack_getri<nm::Rational32>,
+ nm::math::clapack_getri<nm::Rational64>,
+ nm::math::clapack_getri<nm::Rational128>,
+ nm::math::clapack_getri<nm::RubyObject> */
+ };
+
+ if (!ttable[NM_DTYPE(a)]) {
+ rb_raise(rb_eNotImpError, "this operation not yet implemented for non-BLAS dtypes");
+ // FIXME: Once BLAS dtypes are implemented, replace error above with the error below.
+ //rb_raise(nm_eDataTypeError, "this matrix operation undefined for integer matrices");
+ } else {
+ // Call either our version of getri or the LAPACK version.
+ ttable[NM_DTYPE(a)](blas_order_sym(order), blas_uplo_sym(uplo), FIX2INT(n), NM_STORAGE_DENSE(a)->elements, FIX2INT(lda));
+ }
+
+ return a;
+}
+
+
/*
* Call any of the clapack_xlaswp functions as directly as possible.
*
Oops, something went wrong.

0 comments on commit 4c4a064

Please sign in to comment.