522 lines
18 KiB
Plaintext
522 lines
18 KiB
Plaintext
;; @module gsl.lsp
|
|
;; @description Selected functions from the GNU Scientific Library
|
|
;; @version 1.0 - initial release. Minimum newLISP version is 10.4.0.
|
|
;; @version 1.1 - added check for extended ffi enabled version.
|
|
;; @version 1.2 - changed ostype Win32 to Windows
|
|
;; @author Lutz Mueller 2012
|
|
|
|
;; <h2>Module GSL For The GNU Scientific Library</h2>
|
|
;; The @link http://www.gnu.org/software/gsl/ GSL GNU Scientific Library
|
|
;; implements over a 1000 functions from different subject areas. In this
|
|
;; release of the 'gsl.lsp' module only a few linear algebra functions
|
|
;; are implemented.
|
|
|
|
;; To use this module, the main GSL library 'libgsl' and a supporting
|
|
;; library 'libgslcblas' must be installed.
|
|
;;
|
|
;; This module requires newLISP version 10.4.0 as a minimum. For 64-bit
|
|
;; newLISP on Mac OSX, Linux and other Unix, 10.4.2 is the minimum.
|
|
;;
|
|
;; Precompiled 32-bit libraries for the binary distributions of newLISP
|
|
;; versions are available in the
|
|
;; @link http://www.nuevatec.com/GSL/ http://www.nuevatec.com/GSL/
|
|
;; directory. See the 'INSTALL.txt' file in that directory for instructions how
|
|
;; to install the library files.
|
|
;;
|
|
;; The module contains <tt>(test-gsl)</tt> to test all implemented functions.
|
|
|
|
;; @syntax (gsl:CholeskyD <matrix-A>)
|
|
;; @param <matrix-A> A square matrix of m row vectors with n = m column elements each.
|
|
;; @return The matrix L .
|
|
;; The function performs a Cholesky decomposition of <matrix-A>. The matrix A must
|
|
;; be symmetric and positive definite square.
|
|
;;
|
|
;; A = L * Lt
|
|
;;
|
|
;; Lt is the transposition of L.
|
|
;;
|
|
;; @example
|
|
;; (gsl:CholeskyD '((4 2 -2) (2 10 2) (-2 2 5)))
|
|
;;
|
|
;; gsl:Cholesky-L =>
|
|
;;
|
|
;; ( ( 2 0 0 )
|
|
;; ( 1 3 0 )
|
|
;; (-1 1 1.732050808) )
|
|
;;
|
|
|
|
;; @syntax (gsl:Cholesky-solve <matrix-A> <vector-b>)
|
|
;; @param <matrix-A> A square matrix of m row vectors with n = m column elements each.
|
|
;; @params <vector-b> A vector of n elements.
|
|
;; @return Vector x containing a solution for Ax = b .
|
|
;;
|
|
;; @example
|
|
;; (gsl:Cholesky-solve '((4 2 -2) (2 10 2) (-2 2 5)) '(1 2 3))
|
|
;;
|
|
;; => (0.8333333333 -0.1666666667 1)
|
|
;;
|
|
|
|
;; @syntax (gsl:QRD <matrix-A>)
|
|
;; @param <matrix-A> A list of m row vector lists with n column elements each.
|
|
;; @return Vector <em>tau</em> of k = min(n, m) elements.
|
|
;; The function does QR decomposition of a matrix A.
|
|
;;
|
|
;; A = Q * R
|
|
;;
|
|
;; The function works for both, square and rectangular matrices.
|
|
;; The number of rows m in A must be equal to or greater than the number of columns n.
|
|
;; The orthogonal m * m matrix Q can be found in the variable 'gsl:QR-Q'.
|
|
;; The m * n matrix R can be found in the variable 'gsl:QR-R'.
|
|
;;
|
|
;; @example
|
|
;; (gsl:QRD '((12 -51 4) (6 167 -68) (-4 24 -41)) ) => (1.857142857 1.993846154 0)
|
|
;;
|
|
;; gsl:QR-Q =>
|
|
;;
|
|
;; ( ( 0.8571428571 -0.3942857143 -0.3314285714 )
|
|
;; ( 0.4285714286 0.9028571429 0.03428571429)
|
|
;; (-0.2857142857 0.1714285714 -0.9428571429 ) )
|
|
;;
|
|
;; gsl:QR-R =>
|
|
;;
|
|
;; ( (14 21 -14 )
|
|
;; ( 0 175 -70 )
|
|
;; ( 0 0 35 ) )
|
|
|
|
;; @syntax (gsl:QR-solve <matrix-A> <vector-b>)
|
|
;; @param <matrix-A> A matrix of m row vectors with n column elements each.
|
|
;; @params <vector-b> A vector of n elements.
|
|
;; @return Vector x containing a solution for Ax = b .
|
|
;; Matrix A is either square or overdetermined with m > n, more rows than columns.
|
|
;; When m > n, then the variable 'gsl:QR-residual' contains a vector of residuals.
|
|
;; For a square matrix A this vector contains 0's.
|
|
;;
|
|
;; @example
|
|
;; (gsl:QR-solve '((12 -51 4) (6 167 -68) (-4 24 -41)) '(1 2 3))
|
|
;;
|
|
;; => (0.009387755102 -0.02432653061 -0.08832653061)
|
|
;;
|
|
;; (gsl:QR-solve '((1 2) (3 4) (5 6) (7 8)) '(1 2 3 4))
|
|
;; => (9.690821045e-16 0.5)
|
|
;;
|
|
;; gsl:QR-residual
|
|
;; => (2.512815377e-17 1.103917396e-16 -2.961679405e-16 1.606480472e-16)
|
|
|
|
;; @syntax (gsl:SVD <matrix-A>)
|
|
;; @param <matrix-A> A matrix of m row vectors with n column elements each.
|
|
;; @return A vector of diagonal elements from the S matrix.
|
|
;; The function does a SVD (Singular Value Decomposition) of the <matrix-A> into
|
|
;; its components U, S and V. Matrix U is orthogonal m*n. S is a diagonal
|
|
;; square matrix of n singular values. V is a n*n square orthogonal matrix.
|
|
;; The function works for both, square and rectangular matrices.
|
|
;;
|
|
;; A = U * S * Vt
|
|
;;
|
|
;; Vt is the transposition of V.
|
|
;;
|
|
;; The number of rows m in A must be equal to or greater than the number of columns n.
|
|
;; The m*n matrix U can be found in the variable 'gsl:SVD-U'.
|
|
;; The diagonal elements of matrix S can be found in the vector variable 'gsl:SVD-S'.
|
|
;; The n*n matrix V can be found in the variable 'gsl:SVD-V'.
|
|
;;
|
|
;; @example
|
|
;; (gsl:SVD '((1 2) (3 4) (5 6) (7 8))) => (14.2690955 0.6268282324)
|
|
;;
|
|
;; gsl:SVD-U =>
|
|
;; ( (-0.1524832333 -0.8226474722 )
|
|
;; (-0.3499183718 -0.4213752877 )
|
|
;; (-0.5473535103 -0.02010310314)
|
|
;; (-0.7447886488 0.3811690814 ) )
|
|
;;
|
|
;; gsl:SVD-S =>
|
|
;;
|
|
;; (14.2690955 0.6268282324)
|
|
;;
|
|
;; gsl:SVD-V =>
|
|
;;
|
|
;; ( (-0.641423028 0.7671873951)
|
|
;; (-0.7671873951 -0.641423028 ) )
|
|
|
|
;; @syntax (gsl:SVD-solve <matrix-A> <vector-b>)
|
|
;; @param <matrix-A> A matrix of m row vectors with n column elements each.
|
|
;; @params <vector-b> A vector of n elements
|
|
;; @return a vector x containing a solution for Ax = b .
|
|
;; The number of rows m in A must be equal to or greater than the number n of columns.
|
|
;;
|
|
;; @example
|
|
;; (gsl:SV-solve '((1 2) (3 4) (5 6) (7 8)) '(1 2 3 4))
|
|
;;
|
|
;; => (5.551115123e-16 0.5 0 0)
|
|
;;
|
|
|
|
; on Mac OSX when using 32-bit newLISP from the normal binary Mac OSX installer
|
|
; use the following to make the libraries:
|
|
;
|
|
; ./configure CFLAGS="-O2 -m32"
|
|
; make
|
|
; sudo make install
|
|
;
|
|
; install needs the root password
|
|
;
|
|
|
|
(when (!= 1024 (& 1024 (sys-info -1)))
|
|
(println "This module needs newLISP compiled for extended ffi.")
|
|
(println "Must exit")
|
|
(exit))
|
|
|
|
; the following assumes the libararies installed in the system library path
|
|
(set 'LIB
|
|
(if
|
|
(= ostype "Windows") "libgsl-0.dll" ; 32-bit or 64-bit
|
|
(= ostype "OSX") "libgsl.dylib" ; 32-bit or 64-bit
|
|
(= ostype "Linux") "/usr/local/lib/libgsl.so" ; 32-bit or 64-bit
|
|
))
|
|
|
|
; load libgslcblas which contans functions referenced by libgsl
|
|
; the symbol cblas_sdsdot is not needed but newLISP versions before
|
|
; 10.4.2 can not use the 'import' statement without a function name
|
|
; libgslcblas is mainly needed internally by libgsl.
|
|
; On windows the library is automatically loaded by libgsl-0.dll.
|
|
(if
|
|
(= ostype "OSX") (import "libgslcblas.dylib" "cblas_sdsdot")
|
|
(= ostype "Linux") (import "/usr/local/lib/libgslcblas.so" "cblas_sdsdot")
|
|
)
|
|
|
|
; structs are defined but only needed for debugging, instead use "void*"
|
|
(struct 'complex "double" "double") ; complex numbers
|
|
(struct 'block "long" "void*") ; vectors and matrices allocation block, size data-ptr
|
|
(struct 'vector "long" "long" "void*" "block" "int") ; size stride data block owner
|
|
(struct 'matrix "long" "long" "long" "void*" "block" "int") ; size1 size2 tda data block owner
|
|
|
|
(import LIB "gsl_set_error_handler_off" "void*") ; turn of error handler
|
|
(import LIB "gsl_strerror" "char*" "int") ; get error text
|
|
|
|
; pointers are only used passed around internally
|
|
; for debugging (unpack vector vptr) can always be used
|
|
(import LIB "gsl_vector_alloc" "void*" "long")
|
|
(import LIB "gsl_vector_calloc" "void*" "long")
|
|
(import LIB "gsl_vector_free" "void" "void*")
|
|
(import LIB "gsl_vector_get" "double" "void*" "long")
|
|
(import LIB "gsl_vector_set" "void" "void*" "long" "double")
|
|
|
|
; in matrix functions "void*" can be used instead of struct "matrix"
|
|
; because pointers are only used passed around internally
|
|
; for debugging (unpack matrix Xptr) can always be used
|
|
(import LIB "gsl_matrix_alloc" "void*" "long" "long")
|
|
(import LIB "gsl_matrix_calloc" "void*" "long" "long")
|
|
(import LIB "gsl_matrix_free" "void" "void*")
|
|
(import LIB "gsl_matrix_get" "double" "void*" "long" "long")
|
|
(import LIB "gsl_matrix_set" "void" "void*" "long" "long" "double")
|
|
(import LIB "gsl_matrix_scale" "int" "void*" "double")
|
|
|
|
; linear algebra functions
|
|
(import LIB "gsl_linalg_cholesky_decomp" "int" "void*")
|
|
(import LIB "gsl_linalg_cholesky_solve" "int" "void*" "void*" "void*")
|
|
(import LIB "gsl_linalg_QR_decomp" "int" "void*" "void*")
|
|
(import LIB "gsl_linalg_QR_unpack" "int" "void*" "void*" "void*" "void*")
|
|
(import LIB "gsl_linalg_QR_solve" "int" "void*" "void*" "void*" "void*")
|
|
(import LIB "gsl_linalg_QR_lssolve" "int" "void*" "void*" "void*" "void*" "void*")
|
|
(import LIB "gsl_linalg_SV_decomp" "int" "void*" "void*" "void*" "void*")
|
|
(import LIB "gsl_linalg_SV_solve" "int" "void*" "void*" "void*" "void*" "void*")
|
|
|
|
; turn of error handler
|
|
; instead check return values from gsl_xxx_xxx functions
|
|
; and do 'throw-error' retrieving error text with 'gsl_strerror'
|
|
(gsl_set_error_handler_off)
|
|
|
|
; helper functions for translating C-arrays and vectors to lists and back
|
|
|
|
(define (get-vector-from-list x)
|
|
(when (array? x) (set 'x (array-list x)))
|
|
(letn (m (length x) xptr (gsl_vector_calloc m))
|
|
(dotimes (i m) (gsl_vector_set xptr i (pop x)))
|
|
(list xptr m))
|
|
)
|
|
|
|
(define (get-matrix-from-list X)
|
|
(when (array? X) (set 'X (array-list X)))
|
|
(let (m (length X) n (length (X 0)) Xptr 0)
|
|
(set 'Xptr (gsl_matrix_calloc m n))
|
|
(dotimes (i m) (set 'row (pop X)) ; row
|
|
(dotimes (j n) (gsl_matrix_set Xptr i j (pop row)))) ; col
|
|
(list Xptr m n))
|
|
)
|
|
|
|
(define (get-list-from-vector xr)
|
|
(let (xptr (xr 0) n (xr 1) result nil)
|
|
(dotimes (i n)
|
|
(push (gsl_vector_get xptr i) result -1))
|
|
result
|
|
)
|
|
)
|
|
|
|
(define (get-list-from-matrix Xr)
|
|
(let (row nil col nil result nil Xptr (Xr 0) m (Xr 1) n (Xr 2))
|
|
(dotimes (i m)
|
|
(set 'row nil)
|
|
(dotimes (j n)
|
|
(push (gsl_matrix_get Xptr i j) row -1))
|
|
(push row result -1))
|
|
result
|
|
)
|
|
)
|
|
|
|
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; API functions ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
|
|
|
|
|
(define (gsl:CholeskyD A)
|
|
(letn (
|
|
Ar (get-matrix-from-list A)
|
|
Aptr (Ar 0)
|
|
m (Ar 1)
|
|
n (Ar 2)
|
|
result 0
|
|
)
|
|
(set 'result (gsl_linalg_cholesky_decomp Aptr))
|
|
(unless (zero? result)
|
|
(throw-error (gsl_strerror result)))
|
|
; zero the triangle above diagonal
|
|
(for (i 0 (- m 2))
|
|
(for (j (+ i 1) (- n 1))
|
|
(gsl_matrix_set Aptr i j 0.0)))
|
|
(set 'result (get-list-from-matrix Ar))
|
|
(gsl_matrix_free Aptr)
|
|
result) ; matix L of A = LLt
|
|
)
|
|
|
|
(define (gsl:Cholesky-solve A b)
|
|
(letn (
|
|
Ar (get-matrix-from-list A)
|
|
Aptr (Ar 0)
|
|
m (Ar 1)
|
|
n (Ar 2)
|
|
br (get-vector-from-list b)
|
|
bptr (br 0)
|
|
xptr (gsl_vector_calloc n)
|
|
result 0
|
|
)
|
|
(set 'result (gsl_linalg_cholesky_decomp Aptr))
|
|
(unless (zero? result)
|
|
(throw-error (gsl_strerror result)))
|
|
(set 'result (gsl_linalg_cholesky_solve Aptr bptr xptr))
|
|
(unless (zero? result)
|
|
(throw-error (gsl_strerror result)))
|
|
(set 'result (get-list-from-vector (list xptr n)))
|
|
(gsl_matrix_free Aptr)
|
|
(gsl_vector_free bptr)
|
|
(gsl_vector_free xptr)
|
|
result)
|
|
)
|
|
|
|
(define (gsl:QRD A flag)
|
|
(letn (
|
|
Ar (get-matrix-from-list A)
|
|
Aptr (Ar 0)
|
|
m (Ar 1)
|
|
n (Ar 2)
|
|
nt (min m n)
|
|
Tptr (gsl_vector_calloc nt)
|
|
Qptr (gsl_matrix_calloc m m)
|
|
Rptr (gsl_matrix_calloc m n)
|
|
result 0
|
|
)
|
|
(set 'result (gsl_linalg_QR_decomp Aptr Tptr))
|
|
(unless (zero? result)
|
|
(throw-error (gsl_strerror result)))
|
|
; Aptr contains QR
|
|
(set 'result (gsl_linalg_QR_unpack Aptr Tptr Qptr Rptr))
|
|
(unless (zero? result)
|
|
(throw-error (gsl_strerror result)))
|
|
(set 'gsl:QR-tau (get-list-from-vector (list Tptr nt)))
|
|
; Q and R have reverse signs compared with published
|
|
; examples but values are the same and verify A = QR.
|
|
; reverse signs for Q and R unless flag set true
|
|
(unless flag
|
|
(gsl_matrix_scale Qptr -1)
|
|
(gsl_matrix_scale Rptr -1)
|
|
; avoid -0 in lower triangular
|
|
(for (i 1 (- m 1)) (for (j 0 (- i 1))
|
|
(gsl_matrix_set Rptr i j 0))))
|
|
(set 'gsl:QR-Q (get-list-from-matrix (list Qptr m m)))
|
|
(set 'gsl:QR-R (get-list-from-matrix (list Rptr m n)))
|
|
(gsl_matrix_free Aptr)
|
|
(gsl_vector_free Tptr)
|
|
(gsl_matrix_free Qptr)
|
|
(gsl_matrix_free Rptr)
|
|
gsl:QR-tau); vector tau
|
|
)
|
|
|
|
(define (gsl:QR-solve A b)
|
|
(letn (
|
|
Ar (get-matrix-from-list A)
|
|
Aptr (Ar 0)
|
|
m (Ar 1)
|
|
n (Ar 2)
|
|
nt (min m n)
|
|
tptr (gsl_vector_calloc nt)
|
|
br (get-vector-from-list b)
|
|
bptr (br 0)
|
|
xptr (gsl_vector_calloc n)
|
|
resptr (gsl_vector_calloc m)
|
|
result 0
|
|
)
|
|
(set 'result (gsl_linalg_QR_decomp Aptr tptr))
|
|
(unless (zero? result)
|
|
(throw-error (gsl_strerror result)))
|
|
; Aptr contains QR
|
|
(if (> m n) ; more rows than cols
|
|
(set 'result (gsl_linalg_QR_lssolve Aptr tptr bptr xptr resptr))
|
|
(set 'result (gsl_linalg_QR_solve Aptr tptr bptr xptr)) )
|
|
(unless (zero? result)
|
|
(throw-error (gsl_strerror result)))
|
|
(set 'result (get-list-from-vector (list xptr n)))
|
|
(set 'gsl:QR-residual (get-list-from-vector (list resptr m)))
|
|
(gsl_matrix_free Aptr)
|
|
(gsl_vector_free tptr)
|
|
(gsl_vector_free bptr)
|
|
(gsl_vector_free xptr)
|
|
(gsl_vector_free resptr)
|
|
result)
|
|
)
|
|
|
|
(define (gsl:SVD A)
|
|
(letn (
|
|
Ar (get-matrix-from-list A)
|
|
Aptr (Ar 0)
|
|
m (Ar 1)
|
|
n (Ar 2)
|
|
Sptr (gsl_vector_calloc n)
|
|
Vptr (gsl_matrix_calloc n n)
|
|
work (gsl_vector_calloc n)
|
|
result 0)
|
|
(set 'result (gsl_linalg_SV_decomp Aptr Vptr Sptr work))
|
|
(unless (zero? result)
|
|
(throw-error (gsl_strerror result)))
|
|
(set 'gsl:SVD-U (get-list-from-matrix Ar))
|
|
(set 'gsl:SVD-S (get-list-from-vector (list Sptr n)))
|
|
(set 'gsl:SVD-V (get-list-from-matrix (list Vptr n n)))
|
|
(gsl_matrix_free Aptr)
|
|
(gsl_vector_free Sptr)
|
|
(gsl_matrix_free Vptr)
|
|
(gsl_vector_free work)
|
|
gsl:SVD-S) ; vector of diagonals of S of A = USVt
|
|
)
|
|
|
|
(define (gsl:SV-solve A b)
|
|
(letn (
|
|
Ar (get-matrix-from-list A)
|
|
Aptr (Ar 0)
|
|
m (Ar 1)
|
|
n (Ar 2)
|
|
Sptr (gsl_vector_calloc n)
|
|
Vptr (gsl_matrix_calloc n n)
|
|
work (gsl_vector_calloc n)
|
|
xptr (gsl_vector_calloc n)
|
|
br (get-vector-from-list b)
|
|
bptr (br 0)
|
|
result 0)
|
|
(set 'result (gsl_linalg_SV_decomp Aptr Vptr Sptr work))
|
|
(unless (zero? result)
|
|
(throw-error (gsl_strerror result)))
|
|
(set 'result (gsl_linalg_SV_solve Aptr Vptr Sptr bptr xptr))
|
|
(unless (zero? result)
|
|
(throw-error (gsl_strerror result)))
|
|
(set 'result (get-list-from-vector (list xptr n)))
|
|
(gsl_matrix_free Aptr)
|
|
(gsl_vector_free Sptr)
|
|
(gsl_matrix_free Vptr)
|
|
(gsl_vector_free work)
|
|
(gsl_vector_free bptr)
|
|
(gsl_vector_free xptr)
|
|
result)
|
|
)
|
|
|
|
; User helper functions
|
|
|
|
(define (gsl:diagonal x flag)
|
|
(when (array? x) (set 'x (array-list x)))
|
|
(letn (len (length x) A (array len len (dup 0 len)))
|
|
(dotimes (i len) (setf (A i i) (pop x)))
|
|
(unless flag A (array-list A)))
|
|
)
|
|
|
|
(context MAIN)
|
|
|
|
;;;;;;;;;;;;;;;;;;;; tests ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
|
|
|
(define (test-gsl)
|
|
; Cholesky decomp
|
|
(println)
|
|
(println "========== Cholesky example A = L * Lt")
|
|
(println "(gsl:CholeskyD '((4 2 -2) (2 10 2) (-2 2 5)))" )
|
|
(set 'L (gsl:CholeskyD '((4 2 -1) (2 10 2) (-2 2 5))))
|
|
(println)
|
|
(println "L -> ") (map println L)
|
|
(println)
|
|
(println "verify LLt ->") (map println (multiply L (transpose L)))
|
|
(println)
|
|
(println "(gsl:Cholesky-solve '((4 2 -2) (2 10 2) (-2 2 5)) '(1 2 3))")
|
|
(set 'x (gsl:Cholesky-solve '((4 2 -2) (2 10 2) (-2 2 5)) '(1 2 3)))
|
|
(println)
|
|
(println "x -> " x)
|
|
(println "Ax -> " (multiply '((4 2 -2) (2 10 2) (-2 2 5)) (transpose (list x))))
|
|
|
|
; QR decomp example
|
|
(println)
|
|
(println "========== QR example A = Q * R")
|
|
(println "(gsl:QRD '((12 -51 4) (6 167 -68) (-4 24 -41)) ) tau vector ===> "
|
|
(gsl:QRD '((12 -51 4) (6 167 -68) (-4 24 -41)) ))
|
|
|
|
(println)
|
|
(println "Q -> ") (map println gsl:QR-Q)
|
|
(println)
|
|
(println "R -> ") (map println gsl:QR-R)
|
|
(println)
|
|
(println "verify QR -> ") (map println (multiply gsl:QR-Q gsl:QR-R))
|
|
(println)
|
|
(println "(gsl:QR-solve '((12 -51 4) (6 167 -68) (-4 24 -41)) '(1 2 3))")
|
|
(set 'x (gsl:QR-solve '((12 -51 4) (6 167 -68) (-4 24 -41)) '(1 2 3)))
|
|
(println)
|
|
(println "x -> " x)
|
|
(println "Ax -> " (multiply '((12 -51 4) (6 167 -68) (-4 24 -41)) (transpose (list x))))
|
|
(println "residual -> " gsl:QR-residual)
|
|
(println)
|
|
(println "(gsl:QR-solve '((1 2) (3 4) (5 6) (7 8)) '(1 2 3 4))")
|
|
(set 'x (gsl:QR-solve '((1 2) (3 4) (5 6) (7 8)) '(1 2 3 4)))
|
|
(println)
|
|
(println "x -> " x)
|
|
(println "Ax -> " (multiply '((1 2) (3 4) (5 6) (7 8)) (transpose (list x))))
|
|
(println "residual -> " gsl:QR-residual)
|
|
|
|
; SV decomp example
|
|
(println)
|
|
(println "========== SVD example A = U * S * Vt")
|
|
(println "(gsl:SVD '((1 2) (3 4) (5 6) (7 8))) => " (gsl:SVD '((1 2) (3 4) (5 6) (7 8))))
|
|
(println)
|
|
(println "U -> ") (map println gsl:SVD-U)
|
|
(println)
|
|
(println "diagonal of S -> ") (println gsl:SVD-S)
|
|
(println)
|
|
(println "V -> ") (map println gsl:SVD-V)
|
|
(println)
|
|
(println "verify USVt -> ") (map println
|
|
(multiply (multiply gsl:SVD-U (gsl:diagonal gsl:SVD-S)) (transpose gsl:SVD-V)))
|
|
(println)
|
|
(println "(gsl:SV-solve '((1 2) (3 4) (5 6) (7 8)) '(1 2 3 4))")
|
|
(set 'x (gsl:SV-solve '((1 2) (3 4) (5 6) (7 8)) '(1 2 3 4)))
|
|
(println)
|
|
(println "x -> " x)
|
|
(println "Ax -> " (multiply '((1 2) (3 4) (5 6) (7 8)) (transpose (list x))))
|
|
(println)
|
|
(println "(gsl:SV-solve '((12 -51 4) (6 167 -68) (-4 24 -41)) '(1 2 3))")
|
|
(set 'x (gsl:SV-solve '((12 -51 4) (6 167 -68) (-4 24 -41)) '(1 2 3)))
|
|
(println)
|
|
(println "x -> " x)
|
|
(println "Ax -> " (multiply '((12 -51 4) (6 167 -68) (-4 24 -41)) (transpose (list x))))
|
|
true
|
|
)
|
|
|
|
;(test-gsl) (exit)
|
|
|
|
; eof ;
|