From b7074c4335a9a9378637a6fffb5e8ff0894489b1 Mon Sep 17 00:00:00 2001 From: Liam Healy <liam@thinkpad.local> Date: Tue, 17 Mar 2009 20:48:51 -0400 Subject: [PATCH] LU remultiplication with permutation; unify permutation functions The LU remultiplication tests now permute the answer at the end so that it ends up with the original vector with the elements in the original order. The permutation functions #'permute and #'permute-inverse now include all arguments, including marrays and pointers to C arrays. Formerely, there were separate functions #'permute-vector and #'permute-vector-inverse. SBCL64: TOTAL: 1451 assertions passed, 5 failed, 0 execution errors. CCL64: TOTAL: 1452 assertions passed, 4 failed, 0 execution errors. --- data/permutation.lisp | 56 ++++++++++++++------------- linear-algebra/lu.lisp | 22 ++++++----- tests/lu.lisp | 86 +++++++++++++++++++++++++----------------- tests/permutation.lisp | 2 +- 4 files changed, 96 insertions(+), 70 deletions(-) diff --git a/data/permutation.lisp b/data/permutation.lisp index b1c78118..a77d8516 100644 --- a/data/permutation.lisp +++ b/data/permutation.lisp @@ -1,6 +1,6 @@ ;; Permutations ;; Liam Healy, Sun Mar 26 2006 - 11:51 -;; Time-stamp: <2009-01-11 21:18:15EST permutation.lisp> +;; Time-stamp: <2009-03-17 20:27:01EDT permutation.lisp> ;; $Id$ (in-package :gsl) @@ -177,27 +177,7 @@ ;;;; Applying Permutations ;;;;**************************************************************************** -(defmfun permute (p data stride n) - "gsl_permute" - (((mpointer p) :pointer) (data :pointer) (stride sizet) (n sizet)) - :inputs (p data) - :outputs (data) - :return (data) - :documentation ; FDL - "Apply the permutation p to the array data of - size n with stride stride.") - -(defmfun permute-inverse (p data stride n) - "gsl_permute_inverse" - (((mpointer p) :pointer) (data :pointer) (stride sizet) (n sizet)) - :inputs (p data) - :outputs (data) - :return (data) - :documentation ; FDL - "Apply the inverse of the permutation p to the array data of - size n with stride.") - -(defmfun permute-vector ((p permutation) (v vector)) +(defmfun permute ((p permutation) (v vector) &optional size stride) ("gsl_permute_vector" :type) (((mpointer p) :pointer) ((mpointer v) :pointer)) :definition :generic @@ -212,7 +192,20 @@ identity matrix. The permutation p and the vector v must have the same length.") -(defmfun permute-vector-inverse ((p permutation) (v vector)) +(defmfun permute + (p (data #.+foreign-pointer-class+) &optional (size 1) (stride 1)) + "gsl_permute" + (((mpointer p) :pointer) (data :pointer) (stride sizet) (size sizet)) + :definition :method + :inputs (p data) + :outputs (data) + :return (data) + :documentation ; FDL + "Apply the permutation p to the array data of + size n with stride stride.") + +(defmfun permute-inverse + ((p permutation) (v vector) &optional size stride) ("gsl_permute_vector" :type "_inverse") (((mpointer p) :pointer) ((mpointer v) :pointer)) :definition :generic @@ -227,6 +220,18 @@ identity matrix. The permutation p and the vector v must have the same length.") +(defmfun permute-inverse + (p (data #.+foreign-pointer-class+) &optional (size 1) (stride 1)) + "gsl_permute_inverse" + (((mpointer p) :pointer) (data :pointer) (stride sizet) (n sizet)) + :definition :method + :inputs (p data) + :outputs (data) + :return (data) + :documentation ; FDL + "Apply the inverse of the permutation p to the array data of + size n with stride.") + (defmfun permutation* (p pa pb) "gsl_permutation_mul" (((mpointer p) :pointer) @@ -307,8 +312,7 @@ (loop collect (copy-seq (cl-array perm)) while (permutation-previous perm)))) -(save-test - permutation +(save-test permutation (let ((perm-1 (make-permutation 4 t))) ;maref (maref perm-1 2)) (let ((perm-1 (make-permutation 4 t))) ;cl-array @@ -333,7 +337,7 @@ (set-identity perm-1) (swap-elements perm-1 1 3) (swap-elements perm-1 0 2) - (permute-vector perm-1 intvec) + (permute perm-1 intvec) (cl-array intvec)) (let ((perm-1 (make-permutation 4 t))) ;inversions (set-identity perm-1) diff --git a/linear-algebra/lu.lisp b/linear-algebra/lu.lisp index 258e4d00..637e8b7e 100644 --- a/linear-algebra/lu.lisp +++ b/linear-algebra/lu.lisp @@ -1,6 +1,6 @@ ;; LU decomposition ;; Liam Healy, Thu Apr 27 2006 - 12:42 -;; Time-stamp: <2009-02-23 13:50:29EST lu.lisp> +;; Time-stamp: <2009-03-17 20:42:37EDT lu.lisp> ;; $Id$ (in-package :gsl) @@ -148,13 +148,17 @@ :initial-contents '(1.0d0 2.0d0 3.0d0 4.0d0))))) (generate-all-array-tests lu :doubles - (let ((mat (array-default '(3 3))) - (vec (array-default '3))) - (multiple-value-bind (matrix perm) (lu-decomposition mat) - (let ((x (lu-solve matrix vec perm))) - (cl-array + (let ((matrix (array-default '(4 4))) + (vec (array-default '4 ))) + (multiple-value-bind (matrix perm) + (lu-decomposition matrix) + (let ((x (lu-solve matrix vec perm))) + (cl-array + (permute-inverse + perm (matrix-product-triangular matrix - (matrix-product-triangular - matrix x 1 :Upper :NoTrans :NonUnit) - 1 :Lower :NoTrans :Unit)))))) + (matrix-product-triangular matrix x 1 :upper :notrans :nonunit) + 1 :lower :notrans :unit))))))) + + diff --git a/tests/lu.lisp b/tests/lu.lisp index cda6c41e..ae9b69df 100644 --- a/tests/lu.lisp +++ b/tests/lu.lisp @@ -3,61 +3,79 @@ (in-package :gsl) (LISP-UNIT:DEFINE-TEST LU - (LISP-UNIT::ASSERT-NUMERICAL-EQUAL + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL (LIST - #( 42.730000000000004d0 -17.24d0 43.309999999999995d0)) + #(-39.65999999999999d0 -49.46000000000001d0 + 19.679999999999993d0 -5.549999999999997d0)) (MULTIPLE-VALUE-LIST - (LET ((MAT + (LET ((MATRIX (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS - '((-34.5d0 8.24d0 3.29d0) - (-8.93d0 34.12d0 -6.15d0) - (49.27d0 -13.49d0 32.5d0)))) + '((-34.5d0 8.24d0 3.29d0 -8.93d0) + (34.12d0 -6.15d0 49.27d0 + -13.49d0) + (32.5d0 42.73d0 -17.24d0 + 43.31d0) + (-16.12d0 -8.25d0 21.44d0 + -49.08d0)))) (VEC (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS - '(42.73d0 -17.24d0 43.31d0)))) + '(-39.66d0 -49.46d0 19.68d0 + -5.55d0)))) (MULTIPLE-VALUE-BIND (MATRIX PERM) - (LU-DECOMPOSITION MAT) + (LU-DECOMPOSITION MATRIX) (LET ((X (LU-SOLVE MATRIX VEC PERM))) (CL-ARRAY - (MATRIX-PRODUCT-TRIANGULAR MATRIX - (MATRIX-PRODUCT-TRIANGULAR - MATRIX X 1 :UPPER - :NOTRANS :NONUNIT) - 1 :LOWER :NOTRANS - :UNIT))))))) - (LISP-UNIT::ASSERT-NUMERICAL-EQUAL + (PERMUTE-INVERSE PERM + (MATRIX-PRODUCT-TRIANGULAR + MATRIX + (MATRIX-PRODUCT-TRIANGULAR + MATRIX X 1 :UPPER :NOTRANS + :NONUNIT) + 1 :LOWER :NOTRANS + :UNIT)))))))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL (LIST - #(#C(42.730000000000004d0 -17.240000000000006d0) - #C(43.310000000000024d0 -16.11999999999999d0) - #C(-8.249999999999986d0 21.440000000000005d0))) + #(#C(-39.65999999999999d0 -49.46000000000001d0) + #C(19.679999999999996d0 -5.549999999999995d0) + #C(-8.820000000000006d0 25.370000000000005d0) + #C(-30.580000000000002d0 31.67d0))) (MULTIPLE-VALUE-LIST - (LET ((MAT + (LET ((MATRIX (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) - (49.27d0 -13.49d0 32.5d0 42.73d0 - -17.24d0 43.31d0)))) + 34.12d0 -6.15d0 49.27d0 + -13.49d0) + (34.12d0 -6.15d0 49.27d0 + -13.49d0 32.5d0 42.73d0 + -17.24d0 43.31d0) + (32.5d0 42.73d0 -17.24d0 43.31d0 + -16.12d0 -8.25d0 21.44d0 + -49.08d0) + (-16.12d0 -8.25d0 21.44d0 + -49.08d0 -39.66d0 -49.46d0 + 19.68d0 -5.55d0)))) (VEC (MAKE-MARRAY '(COMPLEX DOUBLE-FLOAT) :INITIAL-CONTENTS - '(42.73d0 -17.24d0 43.31d0 - -16.12d0 -8.25d0 21.44d0)))) + '(-39.66d0 -49.46d0 19.68d0 + -5.55d0 -8.82d0 25.37d0 -30.58d0 + 31.67d0)))) (MULTIPLE-VALUE-BIND (MATRIX PERM) - (LU-DECOMPOSITION MAT) + (LU-DECOMPOSITION MATRIX) (LET ((X (LU-SOLVE MATRIX VEC PERM))) (CL-ARRAY - (MATRIX-PRODUCT-TRIANGULAR MATRIX - (MATRIX-PRODUCT-TRIANGULAR - MATRIX X 1 :UPPER - :NOTRANS :NONUNIT) - 1 :LOWER :NOTRANS - :UNIT))))))) - (LISP-UNIT::ASSERT-NUMERICAL-EQUAL + (PERMUTE-INVERSE PERM + (MATRIX-PRODUCT-TRIANGULAR + MATRIX + (MATRIX-PRODUCT-TRIANGULAR + MATRIX X 1 :UPPER :NOTRANS + :NONUNIT) + 1 :LOWER :NOTRANS + :UNIT)))))))) + (LISP-UNIT:ASSERT-NUMERICAL-EQUAL (LIST #2A((-1.9999999999999998d0 1.0d0) (1.4999999999999998d0 -0.49999999999999994d0))) diff --git a/tests/permutation.lisp b/tests/permutation.lisp index 2f1b886a..6e27e59e 100644 --- a/tests/permutation.lisp +++ b/tests/permutation.lisp @@ -76,7 +76,7 @@ PERM-1 1 3) (SWAP-ELEMENTS PERM-1 0 2) - (PERMUTE-VECTOR + (PERMUTE PERM-1 INTVEC) (CL-ARRAY INTVEC)))) -- GitLab