householder.lisp 4.9 KB
Newer Older
liam's avatar
liam committed
1 2
;; Householder Transformations
;; Liam Healy, Wed May 10 2006 - 10:03
3
;; Time-stamp: <2011-02-19 13:08:23EST householder.lisp>
4
;;
5
;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy
6 7 8 9 10 11 12 13 14 15 16 17 18 19
;; 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 <http://www.gnu.org/licenses/>.
20 21 22 23 24 25 26 27 28 29

;;; For householder-transform, it would be nice to be able to pick the
;;; order of arguments returned (reverse).

(in-package :gsl)

;;;;****************************************************************************
;;;; Householder transformations
;;;;****************************************************************************

liam's avatar
liam committed
30
;;; FDL
31 32
;;; 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
liam's avatar
liam committed
33 34 35 36
;;; 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
37 38 39
;;; rank-1 structure of the Householder matrix to create and apply
;;; Householder transformations efficiently.

liam's avatar
liam committed
40
(defmfun householder-transform (v)
41
  "gsl_linalg_householder_transform"
42
  (((mpointer v) :pointer))
43
  :inputs (v)
44
  :outputs (v)
45
  :c-return (ret :double)
liam's avatar
liam committed
46 47 48 49 50 51
  :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.")
52

liam's avatar
liam committed
53
(defmfun householder-HM (tau v A)
54
  "gsl_linalg_householder_hm"
55
  ((tau :double) ((mpointer v) :pointer) ((mpointer A) :pointer))
56
  :inputs (v A)
57
  :outputs (A)
liam's avatar
liam committed
58 59 60 61
  :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.")
62

liam's avatar
liam committed
63
(defmfun householder-MH (tau v A)
64
  "gsl_linalg_householder_mh"
65
  ((tau :double) ((mpointer v) :pointer) ((mpointer A) :pointer))
66
  :inputs (v A)
67
  :outputs (A)
liam's avatar
liam committed
68 69 70 71
  :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.")
72

liam's avatar
liam committed
73
(defmfun householder-Hv (tau v w)
74
  "gsl_linalg_householder_hv"
75
  ((tau :double) ((mpointer v) :pointer) ((mpointer w) :pointer))
76 77
  :inputs (v w)
  :outputs (w)
liam's avatar
liam committed
78 79 80
  :documentation			; FDL
  "Apply the Householder transformation P defined by
  the scalar tau and the vector v to the vector w.  On
81
  output the result P w is stored in w.")
82 83 84 85 86

;;;;****************************************************************************
;;;; Householder solver for linear systems
;;;;****************************************************************************

87 88 89
(defmfun householder-solve
    (A b &optional x-spec
       &aux
90
       (x (grid:ensure-foreign-array x-spec (grid:dimensions b) t)))
91 92 93 94 95 96 97
  ("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))
liam's avatar
liam committed
98
  :documentation			; FDL
99 100 101
  "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
102
   solution returned in it.  If x-spec is a grid:foreign-array, the solution will
103 104 105 106
   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.")
Liam Healy's avatar
Liam Healy committed
107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126

;;;;****************************************************************************
;;;; 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*))