From 31aeaff69902b0b09a38019549c399043387fc39 Mon Sep 17 00:00:00 2001
From: David Thompson <dthompson2@worcester.edu>
Date: Fri, 7 Nov 2014 09:02:39 -0500
Subject: transform: Use GSL for much faster matrix multiplication.

* sly/transform (transform->pointer): New procedure.
  (transform*): Use GSL.
---
 sly/transform.scm | 23 ++++++++++++-----------
 1 file changed, 12 insertions(+), 11 deletions(-)

diff --git a/sly/transform.scm b/sly/transform.scm
index 96b83e1..80dd238 100644
--- a/sly/transform.scm
+++ b/sly/transform.scm
@@ -22,6 +22,7 @@
 ;;; Code:
 
 (define-module (sly transform)
+  #:use-module (system foreign)
   #:use-module (ice-9 match)
   #:use-module (srfi srfi-1)
   #:use-module (srfi srfi-9)
@@ -29,6 +30,7 @@
   #:use-module (sly math)
   #:use-module (sly math vector)
   #:use-module (sly quaternion)
+  #:use-module (sly wrappers gsl)
   #:export (make-transform null-transform identity-transform
             transform? transform-matrix
             transpose transform-vector2
@@ -124,21 +126,20 @@ null-transform if called without any arguments."
       (%make-transform m3)))
   (reduce add null-transform transforms))
 
+(define (transform->pointer t)
+  (bytevector->pointer (array-contents (transform-matrix t))))
+
 (define (transform* . transforms)
   "Return the product of all given transformation matrices.  Return
 identity-transform if called without any arguments."
   (define (mul a b)
-    (let ((m1 (transform-matrix a))
-          (m2 (transform-matrix b))
-          (m3 (make-4x4-matrix)))
-      (do-ec
-            (: r 4) (: c 4)
-            (let ((x (sum-ec
-                      (: k 4)
-                      (* (array-ref m1 r k)
-                         (array-ref m2 k c)))))
-              (array-set! m3 x r c)))
-      (%make-transform m3)))
+    (let ((result (%make-transform (make-4x4-matrix))))
+      (cblas-sgemm cblas-row-major cblas-no-trans cblas-no-trans
+                   4 4 4 1
+                   (transform->pointer a) 4
+                   (transform->pointer b) 4
+                   0 (transform->pointer result) 4)
+      result))
   (reduce mul identity-transform transforms))
 
 (define translate
-- 
cgit v1.2.3