Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Browse files

Improved cblas_rotg behavior, partially fixing issue #49. Still needs…

… rational-safe algorithm.

Updated History and README with new features
  • Loading branch information...
commit 53895c8d7e8de667f5d493d35eecd558412e9b54 1 parent ad8a325
@mohawkjohn mohawkjohn authored
View
6 History.txt
@@ -48,11 +48,11 @@
* Shortcuts for matrix creation (by @agarie)
- * Access to most ATLAS-implemented LAPACK functions for those with ATLAS' clapack interface: xGETRF, xGETRI, xGETRS, xGESV, xPOTRF, xPOTRI, xPOTRS, xPOSV, xLASWP, xSCAL, xLAUUM
+ * Access to most ATLAS-implemented LAPACK functions for those with ATLAS' CLAPACK interface: xGETRF, xGETRI, xGETRS, xGESV, xPOTRF, xPOTRI, xPOTRS, xPOSV, xLASWP, xSCAL, xLAUUM
- * Access to additional ATLAS-implemented BLAS functions: xTRMM, xSYRK, xHERK
+ * Access to additional ATLAS-implemented BLAS functions: xTRMM, xSYRK, xHERK, xROT, xROTG
- * Non-ATLAS versions of CLAPACK functions: xLASWP, xSCAL, xLAUUM
+ * Non-ATLAS versions of CLAPACK functions: xLASWP, xSCAL, xLAUUM, xROT
* Matrix inversion (LU and Cholesky; requires CLAPACK)
View
12 README.rdoc
@@ -33,15 +33,16 @@ The following features exist in the current version of NMatrix (0.0.2):
* Data types: uint8, int8, int16, int32, int64, float32, float64, complex64, complex128, rational64, rational128
(incomplete)
* Conversion between storage and data types (except from-complex, and from-float-to-rational)
-* Element-wise operations and comparisons
+* Element-wise operations and comparisons for dense and yale
* Matrix-matrix multiplication for dense (using ATLAS) and yale
* Matrix-vector multiplication for dense (using ATLAS)
* Dense and list matrix slicing and referencing
* Matlab .MAT v5 file input
* C and C++ API
-* Level 3 BLAS internal implementations (no library) and ATLAS (with library) access:
- * xGEMM
- * xTRSM
+* BLAS internal implementations (no library) and ATLAS (with library) access:
+ * Level 1: xROT, xROTG (BLAS dtypes only)
+ * Level 2: xGEMV
+ * Level 3: xGEMM, xTRSM
* LAPACK ATLAS access:
* xGETRF, xGETRI, xGETRS, xGESV (Gaussian elimination)
* xPOTRF, xPOTRI, xPOTRS, xPOSV (Cholesky factorization)
@@ -52,12 +53,13 @@ The following features exist in the current version of NMatrix (0.0.2):
* xLAUUM (no LAPACK needed, but BLAS dtypes only)
* LU decomposition
* Matrix inversions (requires LAPACK; BLAS dtypes only)
+* Determinant calculation for BLAS dtypes
=== Planned Features (Short-to-Medium Term)
These are features planned for NMatrix 0.1.0, our first beta.
-* calculation of determinant, trace, and eigenvalues (characteristic polynomial) (0.1.0)
+* calculation of determinant (LAPACK-free), trace, and eigenvalues (characteristic polynomial) (0.1.0)
* exponentials and square roots
* matrix inversions (LAPACK-free)
* matrix decomposition/factorization
View
36 ext/nmatrix/util/math.cpp
@@ -128,7 +128,7 @@ extern "C" {
#endif
static VALUE nm_cblas_rot(VALUE self, VALUE n, VALUE x, VALUE incx, VALUE y, VALUE incy, VALUE c, VALUE s);
- static VALUE nm_cblas_rotg(VALUE self, VALUE a, VALUE b); //, VALUE c, VALUE s);
+ static VALUE nm_cblas_rotg(VALUE self, VALUE ab);
static VALUE nm_cblas_gemm(VALUE self, VALUE order, VALUE trans_a, VALUE trans_b, VALUE m, VALUE n, VALUE k, VALUE vAlpha,
VALUE a, VALUE lda, VALUE b, VALUE ldb, VALUE vBeta, VALUE c, VALUE ldc);
@@ -308,7 +308,7 @@ void nm_math_init_blas() {
cNMatrix_BLAS = rb_define_module_under(cNMatrix, "BLAS");
rb_define_singleton_method(cNMatrix_BLAS, "cblas_rot", (METHOD)nm_cblas_rot, 7);
- rb_define_singleton_method(cNMatrix_BLAS, "cblas_rotg", (METHOD)nm_cblas_rotg, 2);
+ rb_define_singleton_method(cNMatrix_BLAS, "cblas_rotg", (METHOD)nm_cblas_rotg, 1);
rb_define_singleton_method(cNMatrix_BLAS, "cblas_gemm", (METHOD)nm_cblas_gemm, 14);
rb_define_singleton_method(cNMatrix_BLAS, "cblas_gemv", (METHOD)nm_cblas_gemv, 11);
@@ -395,20 +395,17 @@ static inline enum CBLAS_ORDER blas_order_sym(VALUE op) {
* The Givens plane rotation can be used to introduce zero elements into a matrix selectively.
*
* This function differs from most of the other raw BLAS accessors. Instead of providing a, b, c, s as arguments, you
- * should only provide a and b (the inputs). The outputs [a,b,c,s] will be returned in a Ruby Array at the end.
+ * should only provide a and b (the inputs), and you should provide them as a single NVector (or the first two elements
+ * of any dense NMatrix or NVector type, specifically).
*
- * The type for b is inferred from a's type, so make sure these are actually the same Ruby object types.
+ * The outputs [c,s] will be returned in a Ruby Array at the end; the input NVector will also be modified in-place.
*
* If you provide rationals, be aware that there's a high probability of an error, since rotg includes a square root --
* and most rationals' square roots are irrational. You're better off converting to Float first.
*
- * You probably don't want to call this function. Instead, why don't you try rotg, which is more flexible
- * with its arguments? Then you don't have to even worry about the above paragraph.
- *
- * 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!
+ * This function, like the other cblas_ functions, does minimal type-checking.
*/
-static VALUE nm_cblas_rotg(VALUE self, VALUE a, VALUE b) {
+static VALUE nm_cblas_rotg(VALUE self, VALUE ab) {
static void (*ttable[nm::NUM_DTYPES])(void* a, void* b, void* c, void* s) = {
NULL, NULL, NULL, NULL, NULL, // can't represent c and s as integers, so no point in having integer operations.
nm::math::cblas_rotg<float>,
@@ -421,29 +418,26 @@ static VALUE nm_cblas_rotg(VALUE self, VALUE a, VALUE b) {
nm::math::cblas_rotg<nm::RubyObject>
};
- dtype_t dtype = nm_dtype_guess(a);
+ dtype_t dtype = NM_DTYPE(ab);
if (!ttable[dtype]) {
rb_raise(nm_eDataTypeError, "this matrix operation undefined for integer matrices");
return Qnil;
} else {
- void *pA = ALLOCA_N(char, DTYPE_SIZES[dtype]),
- *pB = ALLOCA_N(char, DTYPE_SIZES[dtype]),
- *pC = ALLOCA_N(char, DTYPE_SIZES[dtype]),
+ void *pC = ALLOCA_N(char, DTYPE_SIZES[dtype]),
*pS = ALLOCA_N(char, DTYPE_SIZES[dtype]);
- rubyval_to_cval(a, dtype, pA);
- rubyval_to_cval(b, dtype, pB);
+ // extract A and B from the NVector (first two elements)
+ void* pA = NM_STORAGE_DENSE(ab)->elements;
+ void* pB = (char*)(NM_STORAGE_DENSE(ab)->elements) + DTYPE_SIZES[dtype];
// c and s are output
ttable[dtype](pA, pB, pC, pS);
- VALUE result = rb_ary_new2(4);
- rb_ary_store(result, 0, rubyobj_from_cval(pA, dtype).rval);
- rb_ary_store(result, 1, rubyobj_from_cval(pB, dtype).rval);
- rb_ary_store(result, 2, rubyobj_from_cval(pC, dtype).rval);
- rb_ary_store(result, 3, rubyobj_from_cval(pS, dtype).rval);
+ VALUE result = rb_ary_new2(2);
+ rb_ary_store(result, 0, rubyobj_from_cval(pC, dtype).rval);
+ rb_ary_store(result, 1, rubyobj_from_cval(pS, dtype).rval);
return result;
}
View
20 spec/blas_spec.rb
@@ -82,16 +82,16 @@
y[4].should be_within(1e-4).of(-1.3660254037844386)
end
- end
- end
+ # FIXME: Need to write new Rational algorithm, which doesn't choke quite so often on irrational square roots (which often eventually cancel).
+ it "exposes cblas rotg" do
+ ab = NVector.new(2, [6,-8], dtype)
+ c,s = NMatrix::BLAS::cblas_rotg(ab)
+ ab[0].should be_within(1e-6).of(-10)
+ ab[1].should be_within(1e-6).of(-5.quo(3))
+ c.should be_within(1e-6).of(-3.quo(5))
+ s.should be_within(1e-6).of(4.quo(5))
+ end
- # FIXME: Need a way to force rotg to call :float32 instead of :float64, and :complex64 instead of :complex128
- # FIXME: Need to write new Rational algorithm, which doesn't choke quite so often on irrational square roots (which often eventually cancel).
- it "exposes cblas rotg" do
- a,b,c,s = NMatrix::BLAS::cblas_rotg(6.0, -8.0)
- a.should be_within(1e-6).of(-10)
- b.should be_within(1e-6).of(-5.quo(3))
- c.should be_within(1e-6).of(-3.quo(5))
- s.should be_within(1e-6).of(4.quo(5))
+ end
end
end
Please sign in to comment.
Something went wrong with that request. Please try again.