;; Householder Transformations ;; Liam Healy, Wed May 10 2006 - 10:03 ;; Time-stamp: <2011-02-19 13:08:23EST householder.lisp> ;; ;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy ;; Distributed under the terms of the GNU General Public License ;; ;; This program is free software: you can redistribute it and/or modify ;; it under the terms of the GNU General Public License as published by ;; the Free Software Foundation, either version 3 of the License, or ;; (at your option) any later version. ;; ;; This program is distributed in the hope that it will be useful, ;; but WITHOUT ANY WARRANTY; without even the implied warranty of ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ;; GNU General Public License for more details. ;; ;; You should have received a copy of the GNU General Public License ;; along with this program. If not, see . ;;; For householder-transform, it would be nice to be able to pick the ;;; order of arguments returned (reverse). (in-package :gsl) ;;;;**************************************************************************** ;;;; Householder transformations ;;;;**************************************************************************** ;;; FDL ;;; A Householder transformation is a rank-1 modification of the identity ;;; matrix which can be used to zero out selected elements of a vector. A ;;; Householder matrix P takes the form, ;;; P = I - tau v v^T ;;; where v is a vector (called the Householder vector) and ;;; tau = 2/(v^T v). The functions described in this section use the ;;; rank-1 structure of the Householder matrix to create and apply ;;; Householder transformations efficiently. (defmfun householder-transform (v) "gsl_linalg_householder_transform" (((mpointer v) :pointer)) :inputs (v) :outputs (v) :c-return (ret :double) :return (v ret) :documentation ; FDL "Prepare a Householder transformation P = I - tau v v^T which can be used to zero all the elements of the input vector except the first. Returned values are the transformation, which is stored in the vector v, and the scalar tau.") (defmfun householder-HM (tau v A) "gsl_linalg_householder_hm" ((tau :double) ((mpointer v) :pointer) ((mpointer A) :pointer)) :inputs (v A) :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 matrix A. On output the result P A is stored in A.") (defmfun householder-MH (tau v A) "gsl_linalg_householder_mh" ((tau :double) ((mpointer v) :pointer) ((mpointer A) :pointer)) :inputs (v A) :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 matrix A. On output the result A P is stored in A.") (defmfun householder-Hv (tau v w) "gsl_linalg_householder_hv" ((tau :double) ((mpointer v) :pointer) ((mpointer w) :pointer)) :inputs (v w) :outputs (w) :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.") ;;;;**************************************************************************** ;;;; Householder solver for linear systems ;;;;**************************************************************************** (defmfun householder-solve (A b &optional x-spec &aux (x (grid:ensure-foreign-array x-spec (grid:dimensions b) t))) ("gsl_linalg_HH_svx" "gsl_linalg_HH_solve") ((((mpointer A) :pointer) ((mpointer b) :pointer)) (((mpointer A) :pointer) ((mpointer b) :pointer) ((mpointer x) :pointer))) :inputs (A b) :outputs (x) :return ((or x b)) :documentation ; FDL "Solve the system A x = b directly using Householder transformations. If x-spec is NIL (default), the solution will replace b. If x-spec is T, then an array will be created and the solution returned in it. If x-spec is a grid:foreign-array, the solution will be returned in it. If x-spec is non-NIL, on output the solution is stored in x and b is not modified. The matrix A is destroyed by the Householder transformations. The solution is returned from the function call.") ;;;;**************************************************************************** ;;;; Examples and unit test, from linalg/test.c ;;;;**************************************************************************** (defun test-hh-solve-dim (matrix) "Solve the linear equation using Householder with the supplied matrix and a right-hand side vector which is the reciprocal of one more than the index." (householder-solve (copy matrix) (create-rhs-vector (dim0 matrix)) T)) (save-test householder (test-hh-solve-dim *hilb2*) (test-hh-solve-dim *hilb3*) (test-hh-solve-dim *hilb4*) (test-hh-solve-dim *hilb12*) (test-hh-solve-dim *vander2*) (test-hh-solve-dim *vander3*) (test-hh-solve-dim *vander4*) (test-hh-solve-dim *vander12*))