diff --git a/data/permutation.lisp b/data/permutation.lisp index b1c78118db828be64702dad6061097e62475a9ea..a77d851616d68a52afe7a0b19bfa866395bad2f6 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 258e4d00986dc6f38b2d0c952fc1a81d87bd1e12..637e8b7ecc7d318d845ddf2cb36c7d5207731430 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 cda6c41e758265cbdf6dfdd1bc86d88091e619e5..ae9b69df00f80558fe1545446ee696eddf9f1ece 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 2f1b886a69358869afcf13252daec043bf7917c9..6e27e59e8321bcbe75973c5734865347529ebd32 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))))