Commit 10eda520 authored by Liam Healy's avatar Liam Healy

Port a bunch of linear algebra to ffa

Port qrpt, householder, cholesky, svd.  Modify lu and qr to match the
naming patterns adopted.
parent 3ee9947c
;; Definition of GSLL system
;; Liam Healy
;; Time-stamp: <2008-08-10 22:46:33EDT gsll.asd>
;; Time-stamp: <2008-08-11 23:09:38EDT gsll.asd>
;; $Id$
(asdf:defsystem "gsll"
......@@ -85,12 +85,11 @@
(:file "exponential")
(:file "lu")
(:file "qr")
;(:file "qrpt")
;(:file "svd")
; (:file "cholesky")
; (:file "diagonal")
;(:file "householder")
))
(:file "qrpt")
(:file "svd")
(:file "cholesky")
;(:file "diagonal")
(:file "householder")))
#+no
(:file "eigensystems" :depends-on (init data))
;; Skip fft for now, I'm not sure how it works in C
......
;; Cholesky Decomposition
;; Liam Healy, Wed May 3 2006 - 16:38
;; Time-stamp: <2008-02-17 10:43:34EST cholesky.lisp>
;; Time-stamp: <2008-08-11 22:49:09EDT cholesky.lisp>
;; $Id$
(in-package :gsl)
......@@ -17,9 +17,11 @@
;;; (L y = b, L^T x = y), which can be solved by forward and
;;; back-substitution.
(defmfun cholesky-decomp (A)
"gsl_linalg_cholesky_decomp" (((pointer A) gsl-matrix-c))
:invalidate (A)
(defmfun cholesky-decomposition (A)
"gsl_linalg_cholesky_decomp" (((mpointer A) :pointer))
:inputs (A)
:outputs (A)
:return (A)
:documentation ; FDL
"Factorize the positive-definite square matrix A
into the Cholesky decomposition A = L L^T. On output the diagonal
......@@ -27,25 +29,29 @@
L. The upper triangular part of the input matrix contains
L^T, the diagonal terms being identical for both L and
L^T. If the matrix is not positive-definite then the
decomposition will fail, returning the error code :EDOM.")
decomposition will fail, returning the error EDOM.")
(defmfun cholesky-solve (cholesky b x)
(defmfun cholesky-solve (A x b)
"gsl_linalg_cholesky_solve"
(((pointer cholesky) gsl-matrix-c) ((pointer b) gsl-vector-c)
((pointer x) gsl-vector-c))
:invalidate (x)
(((mpointer A) :pointer) ((mpointer b) :pointer)
((mpointer x) :pointer))
:inputs (A b)
:outputs (x)
:return (x)
:documentation ; FDL
"Solve the system A x = b using the Cholesky
decomposition of A into the matrix cholesky given by
#'cholesky-decomp.")
#'cholesky-decomposition.")
(defmfun cholesky-svx (cholesky x)
(defmfun cholesky-solvex (cholesky x)
"gsl_linalg_cholesky_svx"
(((pointer cholesky) gsl-matrix-c) ((pointer x) gsl-vector-c))
:invalidate (x)
(((mpointer cholesky) :pointer) ((mpointer x) :pointer))
:inputs (A x)
:outputs (x)
:return (x)
:documentation ; FDL
"Solve the system A x = b in-place using the
Cholesky decomposition of A into the matrix cholesky given
by #'cholesky-decomp. On input x should contain
by #'cholesky-decomposition. On input x should contain
the right-hand side B, which is replaced by the solution on
output.")
;; Householder Transformations
;; Liam Healy, Wed May 10 2006 - 10:03
;; Time-stamp: <2008-02-17 11:00:45EST householder.lisp>
;; Time-stamp: <2008-08-11 22:31:45EDT householder.lisp>
;; $Id$
;;; For householder-transform, it would be nice to be able to pick the
......@@ -24,8 +24,8 @@
(defmfun householder-transform (v)
"gsl_linalg_householder_transform"
(((pointer v) gsl-vector-c))
:invalidate (v)
(((mpointer v) :pointer))
:outputs (v)
:c-return (ret :double)
:return (v ret)
:documentation ; FDL
......@@ -36,8 +36,8 @@
(defmfun householder-HM (tau v A)
"gsl_linalg_householder_hm"
((tau :double) ((pointer v) gsl-vector-c) ((pointer A) gsl-matrix-c))
:invalidate (A)
((tau :double) ((mpointer v) :pointer) ((mpointer A) :pointer))
:outputs (A)
:documentation ; FDL
"Apply the Householder matrix P defined by the
scalar tau and the vector v to the left-hand side of the
......@@ -45,8 +45,8 @@
(defmfun householder-MH (tau v A)
"gsl_linalg_householder_mh"
((tau :double) ((pointer v) gsl-vector-c) ((pointer A) gsl-matrix-c))
:invalidate (A)
((tau :double) ((mpointer v) :pointer) ((mpointer A) :pointer))
:outputs (A)
:documentation ; FDL
"Apply the Householder matrix P defined by the
scalar tau and the vector v to the right-hand side of the
......@@ -54,12 +54,12 @@
(defmfun householder-Hv (tau v w)
"gsl_linalg_householder_hv"
((tau :double) ((pointer v) gsl-vector-c) ((pointer w) gsl-vector-c))
((tau :double) ((mpointer v) :pointer) ((mpointer w) :pointer))
:documentation ; FDL
"Apply the Householder transformation P defined by
the scalar tau and the vector v to the vector w. On
output the result P w is stored in w."
:invalidate (w))
:outputs (w))
;;;;****************************************************************************
;;;; Householder solver for linear systems
......@@ -67,9 +67,9 @@
(defmfun householder-solve (A b x)
"gsl_linalg_HH_solve"
(((pointer A) gsl-matrix-c) ((pointer b) gsl-vector-c)
((pointer x) gsl-vector-c))
:invalidate (x A)
(((mpointer A) :pointer) ((mpointer b) :pointer)
((mpointer x) :pointer))
:outputs (x A)
:return (x)
:documentation ; FDL
"Solve the system A x = b directly using
......@@ -79,8 +79,8 @@
(defmfun householder-svx (A x)
"gsl_linalg_HH_svx"
(((pointer A) gsl-matrix-c) ((pointer x) gsl-vector-c))
:invalidate (x A)
(((mpointer A) :pointer) ((mpointer x) :pointer))
:outputs (x A)
:return (x)
:documentation ; FDL
"Solve the system A x = b in-place using
......
;; LU decomposition
;; Liam Healy, Thu Apr 27 2006 - 12:42
;; Time-stamp: <2008-08-10 22:20:37EDT lu.lisp>
;; Time-stamp: <2008-08-11 23:10:56EDT lu.lisp>
;; $Id$
(in-package :gsl)
......@@ -32,7 +32,7 @@
partial pivoting (Golub & Van Loan, Matrix Computations,
Algorithm 3.4.1).")
(defmfun solve-LU ((LU matrix) p (b vector) (x vector))
(defmfun LU-solve ((LU matrix) p (b vector) (x vector))
("gsl_linalg" :complex "_LU_solve")
(((mpointer LU) :pointer) ((mpointer p) :pointer)
((mpointer b) :pointer) ((mpointer x) :pointer))
......@@ -45,7 +45,7 @@
"Solve the square system A x = b using the LU
decomposition of A into (LU, p) given by LU-decomp.")
(defmfun svx-LU ((LU matrix) p (x vector))
(defmfun LU-solvex ((LU matrix) p (x vector))
("gsl_linalg" :complex "_LU_svx")
(((mpointer LU) :pointer) ((mpointer p) :pointer)
((mpointer x) :pointer))
......@@ -60,7 +60,7 @@
(LU, p). On input x should contain the right-hand
side b, which is replaced by the solution on output.")
(defmfun refine-LU ((A matrix) LU p (b vector) (x vector) residual)
(defmfun LU-refine ((A matrix) LU p (b vector) (x vector) residual)
("gsl_linalg" :complex "_LU_refine")
(((mpointer A) :pointer) ((mpointer LU) :pointer)
((mpointer p) :pointer)
......@@ -76,7 +76,7 @@
A x = b, using the LU decomposition of A into (LU,p). The initial
residual r = A x - b is also computed and stored in residual. ")
(defmfun invert-LU ((LU matrix) p inverse)
(defmfun LU-invert ((LU matrix) p inverse)
("gsl_linalg" :complex "_LU_invert")
(((mpointer LU) :pointer) ((mpointer p) :pointer)
((mpointer inverse) :pointer))
......@@ -94,7 +94,7 @@
the same result more efficiently and reliably (consult any
introductory textbook on numerical linear algebra for details).")
(defmfun determinant-LU ((LU matrix) signum)
(defmfun LU-determinant ((LU matrix) signum)
("gsl_linalg" :complex "_LU_det")
(((mpointer LU) gsl-matrix-c) (signum :int))
:c-return :double
......@@ -107,7 +107,7 @@
diagonal elements of U and the sign of the row permutation signum.")
(defmfun log-determinant-LU ((LU matrix))
(defmfun LU-log-determinant ((LU matrix))
("gsl_linalg" :complex "_LU_lndet")
(((mpointer LU) gsl-matrix-c))
:c-return :double
......@@ -140,7 +140,7 @@
(per (permutation dim))
(inv (matrix-double-float (list dim dim))))
(LU-decomposition mmat per)
(invert-lu mmat per inv)
(lu-invert mmat per inv)
(cl-array inv)))
#|
......
;; QR decomposition
;; Liam Healy 2008-02-17 11:05:20EST qr.lisp
;; Time-stamp: <2008-08-10 22:46:16EDT qr.lisp>
;; Time-stamp: <2008-08-11 22:22:09EDT qr.lisp>
;; $Id$
(in-package :gsl)
......@@ -53,7 +53,7 @@
QR-decomp. The least-squares solution for rectangular systems can
be found using QR-lssolve.")
(defmfun QR-svx (QR tau x)
(defmfun QR-solvex (QR tau x)
"gsl_linalg_QR_svx"
(((mpointer QR) :pointer) ((mpointer tau) :pointer)
((mpointer x) :pointer))
......@@ -121,7 +121,7 @@
"Solve the triangular system R x = b for x. It may be useful if the
product b' = Q^T b has already been computed using QR-QTvec}.")
(defmfun QR-Rsvx (QR x)
(defmfun QR-Rsolvex (QR x)
"gsl_linalg_QR_Rsvx"
(((mpointer QR) :pointer) ((mpointer x) :pointer))
:inputs (QR x)
......@@ -181,7 +181,7 @@
:documentation ; FDL
"Solves the triangular system R x = b for the N-by-N matrix R.")
(defmfun R-svx (R x)
(defmfun R-solvex (R x)
"gsl_linalg_R_svx"
(((mpointer R) :pointer) ((mpointer x) :pointer))
:inputs (R x)
......
;; QR with column pivoting
;; Liam Healy, Fri Apr 28 2006 - 16:53
;; Time-stamp: <2008-02-17 11:24:34EST qrpt.lisp>
;; Time-stamp: <2008-08-11 22:21:22EDT qrpt.lisp>
;; $Id$
(in-package :gsl)
(defmfun QRPT-decomp (A tau p signum norm)
(defmfun QRPT-decomposition (A tau p signum norm)
"gsl_linalg_QRPT_decomp"
(((pointer A) gsl-matrix-c) ((pointer tau) gsl-vector-c)
((pointer p) gsl-permutation-c)
(signum :pointer) ((pointer norm) gsl-vector-c))
:invalidate (A tau p norm)
(((mpointer A) :pointer) ((mpointer tau) :pointer)
((mpointer p) :pointer)
(signum :pointer) ((mpointer norm) :pointer))
:inputs (A)
:outputs (A tau p norm)
:return (A tau signum p)
:documentation ; FDL
"Factorizes the M-by-N matrix A into
......@@ -33,13 +34,14 @@
column pivoting (Golub & Van Loan, Matrix Computations, Algorithm
5.4.1).")
(defmfun QRPT-decomp2 (A q r tau p signum norm)
(defmfun QRPT-decomposition* (A q r tau p signum norm)
"gsl_linalg_QRPT_decomp2"
(((pointer A) gsl-matrix-c) ((pointer q) gsl-matrix-c)
((pointer r) gsl-matrix-c) ((pointer tau) gsl-vector-c)
((pointer p) gsl-permutation-c) (signum :pointer)
((pointer norm) gsl-vector-c))
:invalidate (q r norm)
(((mpointer A) :pointer) ((mpointer q) :pointer)
((mpointer r) :pointer) ((mpointer tau) :pointer)
((mpointer p) :pointer) (signum :pointer)
((mpointer norm) :pointer))
:inputs (A)
:outputs (q r norm)
:return (q r p signum)
:documentation ; FDL
"Factorize the matrix A into the decomposition
......@@ -48,20 +50,22 @@
(defmfun QRPT-solve (QR tau p b x)
"gsl_linalg_QRPT_solve"
(((pointer QR) gsl-matrix-c) ((pointer tau) gsl-vector-c)
((pointer p) gsl-permutation-c)
((pointer b) gsl-vector-c) ((pointer x) gsl-vector-c))
:invalidate (x)
(((mpointer QR) :pointer) ((mpointer tau) :pointer)
((mpointer p) :pointer)
((mpointer b) :pointer) ((mpointer x) :pointer))
:inputs (QR tau p b)
:outputs (x)
:return (x)
:documentation ; FDL
"Solve the square system A x = b using the QRP^T
decomposition of A into (QR, tau, p) given by #'QRPT-decomp.")
decomposition of A into (QR, tau, p) given by #'QRPT-decomposition.")
(defmfun QRPT-svx (QR tau p x)
(defmfun QRPT-solvex (QR tau p x)
"gsl_linalg_QRPT_svx"
(((pointer QR) gsl-matrix-c) ((pointer tau) gsl-vector-c)
((pointer p) gsl-permutation-c) ((pointer x) gsl-vector-c))
:invalidate (x)
(((mpointer QR) :pointer) ((mpointer tau) :pointer)
((mpointer p) :pointer) ((mpointer x) :pointer))
:inputs (QR tau p x)
:outputs (x)
:return (x)
:documentation ; FDL
"Solve the square system A x = b in-place using the
......@@ -70,9 +74,10 @@
(defmfun QRPT-QRsolve (QR p b x)
"gsl_linalg_QRPT_QRsolve"
(((pointer QR) gsl-matrix-c) ((pointer p) gsl-permutation-c)
((pointer b) gsl-vector-c) ((pointer x) gsl-vector-c))
:invalidate (x)
(((mpointer QR) :pointer) ((mpointer p) :pointer)
((mpointer b) :pointer) ((mpointer x) :pointer))
:inputs (QR p b)
:outputs (x)
:return (x)
:documentation ; FDL
"Solve the square system R P^T x = Q^T b for
......@@ -81,10 +86,11 @@
(defmfun QRPT-update (Q R p w v)
"gsl_linalg_QRPT_update"
(((pointer Q) gsl-matrix-c) ((pointer R) gsl-matrix-c)
((pointer p) gsl-permutation-c)
((pointer w) gsl-vector-c) ((pointer v) gsl-vector-c))
:invalidate (w Q R)
(((mpointer Q) :pointer) ((mpointer R) :pointer)
((mpointer p) :pointer)
((mpointer w) :pointer) ((mpointer v) :pointer))
:inputs (Q R p w v)
:outputs (w Q R)
:return (Q R)
:documentation ; FDL
"Perform a rank-1 update w v^T of the QRP^T
......@@ -95,19 +101,21 @@
(defmfun QRPT-Rsolve (QR p b x)
"gsl_linalg_QRPT_Rsolve"
(((pointer QR) gsl-matrix-c) ((pointer p) gsl-permutation-c)
((pointer b) gsl-vector-c) ((pointer x) gsl-vector-c))
:invalidate (x)
(((mpointer QR) :pointer) ((mpointer p) :pointer)
((mpointer b) :pointer) ((mpointer x) :pointer))
:inputs (QR p b)
:outputs (x)
:return (x)
:documentation ; FDL
"Solve the triangular system R P^T x = b for the
N-by-N matrix R contained in QR.")
(defmfun QRPT-Rsvx (QR p x)
(defmfun QRPT-Rsolvex (QR p x)
"gsl_linalg_QRPT_Rsvx"
(((pointer QR) gsl-matrix-c) ((pointer p) gsl-permutation-c)
((pointer x) gsl-vector-c))
:invalidate (x)
(((mpointer QR) :pointer) ((mpointer p) :pointer)
((mpointer x) :pointer))
:inputs (QR p x)
:outputs (x)
:return (x)
:documentation ; FDL
"Solve the triangular system R P^T x = b in-place
......
;; Singular Value Decomposition
;; Liam Healy, Tue May 2 2006 - 12:15
;; Time-stamp: <2008-02-17 11:29:50EST svd.lisp>
;; Time-stamp: <2008-08-11 22:46:49EDT svd.lisp>
;; $Id$
(in-package :gsl)
......@@ -25,11 +25,13 @@
;;; precision. Small singular values should be edited by choosing a suitable
;;; tolerance.
(defmfun SV-decomp (A V S work)
(defmfun SV-decomposition (A S V work)
"gsl_linalg_SV_decomp"
(((pointer A) gsl-matrix-c) ((pointer V) gsl-matrix-c)
((pointer S) gsl-vector-c) ((pointer work) gsl-vector-c))
:invalidate (A)
(((mpointer A) :pointer) ((mpointer V) :pointer)
((mpointer S) :pointer) ((mpointer work) :pointer))
:inputs (A)
:outputs (A S V)
:return (A S V)
:documentation ; FDL
"Factorize the M-by-N matrix A into
the singular value decomposition A = U S V^T for M >= N.
......@@ -42,34 +44,40 @@
transpose of V. A workspace of length N is required in work.
This routine uses the Golub-Reinsch SVD algorithm.")
(defmfun SV-decomp-mod (A X V S work)
(defmfun SV-modified-decomposition (A S V X work)
"gsl_linalg_SV_decomp_mod"
(((pointer A) gsl-matrix-c) ((pointer X) gsl-matrix-c)
((pointer V) gsl-matrix-c)
((pointer S) gsl-vector-c) ((pointer work) gsl-vector-c))
:invalidate (A)
(((mpointer A) :pointer) ((mpointer X) :pointer)
((mpointer V) :pointer)
((mpointer S) :pointer) ((mpointer work) :pointer))
:inputs (A)
:outputs (A S V)
:return (A S V)
:documentation ; FDL
"The SVD using the modified Golub-Reinsch algorithm, which is
faster for M >> N. It requires the vector work of length N and the
N-by-N matrix X as additional working space.")
(defmfun SV-decomp-jacobi (A V S)
(defmfun SV-jacobi-decomposition (A S V)
"gsl_linalg_SV_decomp_jacobi"
(((pointer A) gsl-matrix-c) ((pointer V) gsl-matrix-c)
((pointer S) gsl-vector-c))
:invalidate (A)
(((mpointer A) :pointer) ((mpointer V) :pointer)
((mpointer S) :pointer))
:inputs (A)
:outputs (A S V)
:return (A S V)
:documentation ; FDL
"The SVD of the M-by-N matrix A using one-sided Jacobi
orthogonalization for M >= N. The Jacobi method can compute singular
values to higher relative accuracy than Golub-Reinsch algorithms (see
references for details).")
(defmfun SV-solve (U V S b x)
(defmfun SV-solve (U S V b x)
"gsl_linalg_SV_solve"
(((pointer U) gsl-matrix-c) ((pointer V) gsl-matrix-c)
((pointer S) gsl-vector-c)
((pointer b) gsl-vector-c) ((pointer x) gsl-vector-c))
:invalidate (x)
(((mpointer U) :pointer) ((mpointer V) :pointer)
((mpointer S) :pointer)
((mpointer b) :pointer) ((mpointer x) :pointer))
:inputs (U S V b)
:outputs (x)
:return (x)
:documentation ; FDL
"Solve the system A x = b using the singular value
decomposition (U, S, V) of A given by #'SV-decomp.
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment