From b009f3e6347a96c325bab2facecc8ae790defc05 Mon Sep 17 00:00:00 2001
From: Liam Healy <liam@thinkpad.local>
Date: Sun, 15 Mar 2009 18:34:44 -0400
Subject: [PATCH] Non-square matrix multiplication; marking cl-invalid

Fixed a problem in #'matrix-product-dimensions reported by
Norman Werner that prevented taking the product of non-square
matrices, and added a test to check for the multiplication for all
four element types supported.  This revealed a problem in non-native
marray initialization (found in CCL) which has been fixed by marking
foreign-arrays :cl-invalid when neither :initial-contents nor
:initial-element have been specified.
SBCL 64: TOTAL: 1449 assertions passed, 7 failed, 0 execution errors.
CCL64: TOTAL: 1450 assertions passed, 6 failed, 0 execution errors.
---
 data/foreign-array.lisp             |  7 +--
 gsll-tests.asd                      |  3 +-
 linear-algebra/blas2.lisp           |  4 +-
 linear-algebra/blas3.lisp           |  7 ++-
 tests/matrix-product-nonsquare.lisp | 78 +++++++++++++++++++++++++++++
 5 files changed, 92 insertions(+), 7 deletions(-)
 create mode 100644 tests/matrix-product-nonsquare.lisp

diff --git a/data/foreign-array.lisp b/data/foreign-array.lisp
index 9c391ac1..58b38ce9 100644
--- a/data/foreign-array.lisp
+++ b/data/foreign-array.lisp
@@ -1,6 +1,6 @@
 ;; Foreign arrays (usually in C)
 ;; Liam Healy 2008-12-28 10:44:22EST foreign-array.lisp
-;; Time-stamp: <2009-02-23 15:18:04EST foreign-array.lisp>
+;; Time-stamp: <2009-03-15 18:22:52EDT foreign-array.lisp>
 ;; $Id: $
 
 (in-package :gsl)
@@ -52,7 +52,7 @@
 (defmethod initialize-instance :after
     ((object foreign-array) &rest initargs
      &key dimensions initial-contents initial-element)
-  (declare (ignore dimensions initial-contents initial-element))
+  (declare (ignore dimensions #+native initial-contents #+native initial-element))
   (with-slots (cl-array dimensions original-array offset total-size) object
     (unless (and (slot-boundp object 'cl-array) cl-array)
       (setf cl-array (apply #'make-ffa (element-type object) initargs)))
@@ -64,7 +64,8 @@
 		 :count (total-size object))))
       (setf (c-pointer object) cptr)
       (tg:finalize object (lambda () (cffi:foreign-free cptr))))
-    #-native (setf (cl-invalid object) nil)
+    #-native
+    (setf (cl-invalid object) (not (or initial-contents initial-element)))
     (multiple-value-bind  (oa index-offset)
 	(find-original-array (cl-array object))
       (setf original-array oa
diff --git a/gsll-tests.asd b/gsll-tests.asd
index f60c2073..4037481d 100644
--- a/gsll-tests.asd
+++ b/gsll-tests.asd
@@ -1,6 +1,6 @@
 ;; Definition of GSLL system 
 ;; Liam Healy
-;; Time-stamp: <2009-02-24 15:06:55EST gsll-tests.asd>
+;; Time-stamp: <2009-03-15 17:21:18EDT gsll-tests.asd>
 ;; $Id$
 
 (asdf:defsystem "gsll-tests"
@@ -108,6 +108,7 @@
 	     (:file "matrix-mult")
 	     (:file "matrix-product-hermitian")
 	     (:file "matrix-product")
+	     (:file "matrix-product-nonsquare")
 	     (:file "matrix-product-symmetric")
 	     (:file "matrix-product-triangular")
 	     (:file "matrix-set-all-add")
diff --git a/linear-algebra/blas2.lisp b/linear-algebra/blas2.lisp
index 04e9af97..7f43287c 100644
--- a/linear-algebra/blas2.lisp
+++ b/linear-algebra/blas2.lisp
@@ -1,6 +1,6 @@
 ;; BLAS level 2, Matrix-vector operations
 ;; Liam Healy, Wed Apr 26 2006 - 21:08
-;; Time-stamp: <2009-01-15 21:45:11EST blas2.lisp>
+;; Time-stamp: <2009-03-15 17:09:00EDT blas2.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -30,7 +30,7 @@
 (defun matrix-product-dimensions (a b)
   (if (typep b 'matrix)
       (list (first (dimensions a))
-	    (second (dimensions a)))
+	    (second (dimensions b)))
       (first (dimensions a))))
 
 ;;; To do: the y should be an optional argument, default to a zero vector.
diff --git a/linear-algebra/blas3.lisp b/linear-algebra/blas3.lisp
index 3430d0e0..186792bb 100644
--- a/linear-algebra/blas3.lisp
+++ b/linear-algebra/blas3.lisp
@@ -1,6 +1,6 @@
 ;; BLAS level 3, Matrix-matrix operations
 ;; Liam Healy, Wed Apr 26 2006 - 21:08
-;; Time-stamp: <2009-01-15 21:49:10EST blas3.lisp>
+;; Time-stamp: <2009-03-15 17:15:40EDT blas3.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -153,6 +153,11 @@
        (s2 (scalar-default)))
    (cl-array (matrix-product m1 m2 m3 s1 s2))))
 
+(generate-all-array-tests matrix-product-nonsquare :float-complex
+ (let ((m1 (array-default '(2 3)))
+       (m2 (array-default '(3 2))))
+   (cl-array (matrix-product m1 m2))))
+
 (generate-all-array-tests matrix-product-triangular :float-complex
  (let ((m1 (array-default '(3 3)))
        (m2 (array-default '(3 3)))
diff --git a/tests/matrix-product-nonsquare.lisp b/tests/matrix-product-nonsquare.lisp
new file mode 100644
index 00000000..14ee6e71
--- /dev/null
+++ b/tests/matrix-product-nonsquare.lisp
@@ -0,0 +1,78 @@
+;; Regression test MATRIX-PRODUCT-NONSQUARE for GSLL, automatically generated
+
+(in-package :gsl)
+
+(LISP-UNIT:DEFINE-TEST MATRIX-PRODUCT-NONSQUARE
+                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+                        (LIST #2A((-1488.7347 959.9901) (774.94495 1312.0566)))
+                        (MULTIPLE-VALUE-LIST
+                         (LET ((M1
+                                (MAKE-MARRAY 'SINGLE-FLOAT :INITIAL-CONTENTS
+                                             '((-34.5 8.24 3.29)
+                                               (-8.93 34.12 -6.15))))
+                               (M2
+                                (MAKE-MARRAY 'SINGLE-FLOAT :INITIAL-CONTENTS
+                                             '((49.27 -13.49) (32.5 42.73)
+                                               (-17.24 43.31)))))
+                           (CL-ARRAY (MATRIX-PRODUCT M1 M2)))))
+                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+                        (LIST
+                         #2A((-1488.7346d0 959.9901d0)
+                             (774.9448999999998d0 1312.0567999999996d0)))
+                        (MULTIPLE-VALUE-LIST
+                         (LET ((M1
+                                (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
+                                             '((-34.5d0 8.24d0 3.29d0)
+                                               (-8.93d0 34.12d0 -6.15d0))))
+                               (M2
+                                (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
+                                             '((49.27d0 -13.49d0)
+                                               (32.5d0 42.73d0)
+                                               (-17.24d0 43.31d0)))))
+                           (CL-ARRAY (MATRIX-PRODUCT M1 M2)))))
+                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+                        (LIST
+                         #2A((#C(-1422.0259 2305.5098)
+                              #C(-1744.0583 -1092.294))
+                             (#C(-3459.892 1995.4918)
+                              #C(-3290.4465 -801.0577))))
+                        (MULTIPLE-VALUE-LIST
+                         (LET ((M1
+                                (MAKE-MARRAY '(COMPLEX SINGLE-FLOAT)
+                                             :INITIAL-CONTENTS
+                                             '((-34.5 8.24 3.29 -8.93 34.12
+                                                -6.15)
+                                               (-8.93 34.12 -6.15 49.27 -13.49
+                                                32.5))))
+                               (M2
+                                (MAKE-MARRAY '(COMPLEX SINGLE-FLOAT)
+                                             :INITIAL-CONTENTS
+                                             '((49.27 -13.49 32.5 42.73)
+                                               (32.5 42.73 -17.24 43.31)
+                                               (-17.24 43.31 -16.12 -8.25)))))
+                           (CL-ARRAY (MATRIX-PRODUCT M1 M2)))))
+                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+                        (LIST
+                         #2A((#C(-1422.0258d0 2305.5097000000005d0)
+                              #C(-1744.0584d0 -1092.2939d0))
+                             (#C(-3459.8918d0 1995.4917d0)
+                              #C(-3290.4465d0 -801.0577000000002d0))))
+                        (MULTIPLE-VALUE-LIST
+                         (LET ((M1
+                                (MAKE-MARRAY '(COMPLEX DOUBLE-FLOAT)
+                                             :INITIAL-CONTENTS
+                                             '((-34.5d0 8.24d0 3.29d0 -8.93d0
+                                                34.12d0 -6.15d0)
+                                               (-8.93d0 34.12d0 -6.15d0 49.27d0
+                                                -13.49d0 32.5d0))))
+                               (M2
+                                (MAKE-MARRAY '(COMPLEX DOUBLE-FLOAT)
+                                             :INITIAL-CONTENTS
+                                             '((49.27d0 -13.49d0 32.5d0
+                                                42.73d0)
+                                               (32.5d0 42.73d0 -17.24d0
+                                                43.31d0)
+                                               (-17.24d0 43.31d0 -16.12d0
+                                                -8.25d0)))))
+                           (CL-ARRAY (MATRIX-PRODUCT M1 M2))))))
+
-- 
GitLab