lu.lisp 7.47 KB
Newer Older
liam's avatar
liam committed
1
2
;; LU decomposition
;; Liam Healy, Thu Apr 27 2006 - 12:42
3
;; Time-stamp: <2010-02-05 09:15:11EST lu.lisp>
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
;;
;; Copyright 2006, 2007, 2008, 2009 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 <http://www.gnu.org/licenses/>.
20
21
22

(in-package :gsl)

23
24
;;; /usr/include/gsl/gsl_linalg.h

25
26
(defmfun LU-decomposition
    ((A matrix) &optional (permutation (make-permutation (dim0 A))))
Liam Healy's avatar
Liam Healy committed
27
  ("gsl_linalg" :complex "_LU_decomp")
28
29
  (((mpointer A) :pointer) ((mpointer permutation) :pointer)
   (signum (:pointer :int)))
Liam Healy's avatar
Liam Healy committed
30
  :definition :generic
31
  :inputs (A)
32
33
  :outputs (A permutation)
  :return (A permutation signum)
Liam Healy's avatar
Liam Healy committed
34
  :element-types :doubles
liam's avatar
liam committed
35
  :documentation			; FDL
Liam Healy's avatar
Liam Healy committed
36
37
38
39
40
41
  "Factorize the square matrix A into the LU decomposition PA = LU,
  and return the sign of the permutation.  On output the diagonal and
  upper triangular part of the input matrix A contain the matrix U.
  The lower triangular part of the input matrix (excluding the
  diagonal) contains L.  The diagonal elements of L are unity, and are
  not stored.
42

43
44
  The permutation matrix P is encoded in the permutation supplied as
  the second argument and returned as the second value.  The j-th
Liam Healy's avatar
Liam Healy committed
45
46
  column of the matrix P is given by the k-th column of the identity
  matrix, where k = p_j the j-th element of the permutation
47
48
49
  vector. The sign of the permutation is returned as the second value;
  it is the value (-1)^n, where n is the number of interchanges in the
  permutation.
50
51

  The algorithm used in the decomposition is Gaussian Elimination with
liam's avatar
liam committed
52
53
  partial pivoting (Golub & Van Loan, Matrix Computations,
  Algorithm 3.4.1).")
54

55
56
57
58
59
60
61
62
63
64
65
66
67
(defmfun LU-solve
    ((A matrix) (b vector) permutation &optional x-spec
     &aux
     (x (if (eq x-spec t)
	    (make-marray element-type :dimensions (dimensions b))
	    x-spec)))
  (("gsl_linalg" :complex "_LU_svx")
   ("gsl_linalg" :complex "_LU_solve"))
  ((((mpointer A) :pointer)
    ((mpointer permutation) :pointer) ((mpointer b) :pointer))
   (((mpointer A) :pointer)
    ((mpointer permutation) :pointer) ((mpointer b) :pointer)
    ((mpointer x) :pointer)))
Liam Healy's avatar
Liam Healy committed
68
69
  :definition :generic
  :element-types :doubles
70
71
72
  :inputs (A b permutation)
  :outputs (x b)		  ; depends on switch; both to be sure
  :return ((or x b))
liam's avatar
liam committed
73
74
  :documentation			; FDL
  "Solve the square system A x = b using the LU
75
76
77
78
  decomposition of A into (LU, p) given by LU-decomp.
  If x-spec is nil, the solution will be computed in-place replacing b,
  if it is T, an appropriate vector will be created and the solution
  will be computed there.  Otherwise it should be a supplied vector.")
79

Liam Healy's avatar
Liam Healy committed
80
81
82
(defmfun LU-refine
    ((A matrix) LU p (b vector) (x vector)
     &optional (residual (make-marray element-type :dimensions (dim0 A))))
Liam Healy's avatar
Liam Healy committed
83
84
85
86
87
88
89
90
91
92
  ("gsl_linalg" :complex "_LU_refine")
  (((mpointer A) :pointer) ((mpointer LU) :pointer)
   ((mpointer p) :pointer)
   ((mpointer b) :pointer) ((mpointer x) :pointer)
   ((mpointer residual) :pointer))
  :definition :generic
  :inputs (LU p)
  :outputs (x residual)
  :return (x)
  :element-types :doubles
liam's avatar
liam committed
93
94
  :documentation			; FDL
  "Apply an iterative improvement to x, the solution of
95
  A x = b, using the LU decomposition of A into (LU,p). The initial
liam's avatar
liam committed
96
  residual r = A x - b is also computed and stored in residual. ")
97

98
(defmfun LU-invert ((LU matrix) p inverse)
Liam Healy's avatar
Liam Healy committed
99
100
101
102
103
104
105
  ("gsl_linalg" :complex "_LU_invert")
  (((mpointer LU) :pointer) ((mpointer p) :pointer)
   ((mpointer inverse) :pointer))
  :definition :generic
  :inputs (LU p)
  :outputs (inverse)
  :element-types :doubles
liam's avatar
liam committed
106
107
  :documentation			; FDL
  "Compute the inverse of a matrix A from its LU
108
109
110
111
112
   decomposition (LU,p), storing the result in the matrix inverse. The
   inverse is computed by solving the system A x = b for each column of
   the identity matrix. It is preferable to avoid direct use of the
   inverse whenever possible, as the linear solver functions can obtain
   the same result more efficiently and reliably (consult any
liam's avatar
liam committed
113
   introductory textbook on numerical linear algebra for details).")
114

115
 (defmfun LU-determinant ((LU matrix) signum)
Liam Healy's avatar
Liam Healy committed
116
  ("gsl_linalg" :complex "_LU_det")
117
  (((mpointer LU) :pointer) (signum :int))
Liam Healy's avatar
Liam Healy committed
118
119
120
  :definition :generic
  :inputs (LU)
  :element-types :doubles
Liam Healy's avatar
Liam Healy committed
121
  :c-return :double
Liam Healy's avatar
Liam Healy committed
122
123
  :documentation			; FDL
  "Compute the determinant of a matrix from its LU
124
  decomposition, LU. The determinant is computed as the product of the
Liam Healy's avatar
Liam Healy committed
125
126
  diagonal elements of U and the sign of the row permutation signum.")

127
(defmfun LU-log-determinant ((LU matrix))
Liam Healy's avatar
Liam Healy committed
128
  ("gsl_linalg" :complex "_LU_lndet")
129
  (((mpointer LU) :pointer))
Liam Healy's avatar
Liam Healy committed
130
131
132
  :definition :generic
  :inputs (LU)
  :element-types :doubles
Liam Healy's avatar
Liam Healy committed
133
  :c-return :double
liam's avatar
liam committed
134
135
136
  :documentation			; FDL
  "The logarithm of the absolute value of the
   determinant of a matrix A, ln|det(A)|, from its LU decomposition,
137
   LU. This function may be useful if the direct computation of the
liam's avatar
liam committed
138
   determinant would overflow or underflow.")
139

Liam Healy's avatar
Liam Healy committed
140
141
(defmfun LU-sgndet ((LU matrix) signum)
  ("gsl_linalg" :complex "_LU_sgndet")
142
  (((mpointer LU) :pointer) (signum (:pointer :int)))
Liam Healy's avatar
Liam Healy committed
143
144
145
  :definition :generic
  :inputs (LU)
  :element-types :doubles
Liam Healy's avatar
Liam Healy committed
146
  :c-return :int
liam's avatar
liam committed
147
  :documentation 			; FDL
148
  "Compute the sign or phase factor of the determinant of a matrix A,
liam's avatar
liam committed
149
  det(A)/|det(A)|, from its LU decomposition, LU.")
150

Liam Healy's avatar
Liam Healy committed
151
;;; Invert a matrix using LU
152
(export 'invert-matrix)
153
154
(defun invert-matrix (mat)
  "Invert the matrix."
155
  (let* ((dim (first (dimensions mat)))
156
	 (per (make-permutation dim))
157
	 (inv (make-marray 'double-float :dimensions (list dim dim))))
158
159
    (LU-decomposition mat per)
    (lu-invert mat per inv)))
160

Liam Healy's avatar
Liam Healy committed
161
162
163
164
;;; Examples and unit test

;;; These are direct tests of matrix inversion

165
166
167
168
169
170
171
(save-test
 lu
 (cl-array
  (invert-matrix
   (make-marray 'double-float
		:dimensions  '(2 2)
		:initial-contents '(1.0d0 2.0d0 3.0d0 4.0d0)))))
172
173

(generate-all-array-tests lu :doubles
174
175
176
177
178
179
180
181
 (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
182
183
	(matrix-product-triangular
	 matrix
184
185
186
	 (matrix-product-triangular matrix x 1 :upper :notrans :nonunit)
	 1 :lower :notrans :unit)))))))

Liam Healy's avatar
Liam Healy committed
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
;;; From linalg/test.c

(defun test-lu-solve-dim (matrix &optional vector)
  "Solve the linear equation using LU with the supplied matrix and
   a right-hand side vector which is the reciprocal of one more than
   the index."
  (let* ((dim (dim0 matrix))
	 (rhs (or vector
		  (create-rhs-vector dim (element-type matrix)))))
    (multiple-value-bind (upper permutation signum)
	(LU-decomposition (copy matrix))
      (declare (ignore signum))
      (let ((initial-solution (LU-solve upper rhs permutation T)))
	(LU-refine matrix upper permutation rhs initial-solution)))))

(save-test lu
203
204
205
206
207
208
209
210
 (test-lu-solve-dim *hilb2*)
 (test-lu-solve-dim *hilb3*)
 (test-lu-solve-dim *hilb4*)
 (test-lu-solve-dim *hilb12*)
 (test-lu-solve-dim *vander2*)
 (test-lu-solve-dim *vander3*)
 (test-lu-solve-dim *vander4*)
 (test-lu-solve-dim *vander12*)
Liam Healy's avatar
Liam Healy committed
211
 (test-lu-solve-dim (create-complex-matrix 7)))