lu_base.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /* linalg/lu.c
24  *
25  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
26  *
27  * This program is free software; you can redistribute it and/or modify
28  * it under the terms of the GNU General Public License as published by
29  * the Free Software Foundation; either version 3 of the License, or (at
30  * your option) any later version.
31  *
32  * This program is distributed in the hope that it will be useful, but
33  * WITHOUT ANY WARRANTY; without even the implied warranty of
34  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
35  * General Public License for more details.
36  *
37  * You should have received a copy of the GNU General Public License
38  * along with this program; if not, write to the Free Software
39  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
40  * 02110-1301, USA.
41  */
42 /** \file lu_base.h
43  \brief Functions related to LU decomposition
44 */
45 
46 #ifdef DOXYGEN
47 namespace o2scl_linalg {
48 #endif
49 
50  /** \brief Return 1 if at least one of the elements in the
51  diagonal is zero
52  */
53  template<class mat_t>
54  int diagonal_has_zero(const size_t N, mat_t &A) {
55 
56  for(size_t i=0;i<N;i++) {
57  if (O2SCL_IX2(A,i,i)==0.0) return 1;
58  }
59  return 0;
60  }
61 
62  /** \brief Compute the LU decomposition of the matrix \c A
63 
64  On output the diagonal and upper triangular part of the input
65  matrix A contain the matrix U. The lower triangular part of the
66  input matrix (excluding the diagonal) contains L. The diagonal
67  elements of L are unity, and are not stored.
68 
69  The permutation matrix P is encoded in the permutation p. The
70  j-th column of the matrix P is given by the k-th column of the
71  identity matrix, where k = p_j the j-th element of the
72  permutation vector. The sign of the permutation is given by
73  signum. It has the value (-1)^n, where n is the number of
74  interchanges in the permutation.
75 
76  The algorithm used in the decomposition is Gaussian Elimination
77  with partial pivoting (Golub & Van Loan, Matrix Computations,
78  Algorithm 3.4.1).
79 
80  \future The "swap rows j and i_pivot" section could probably
81  be made more efficient using a "matrix_row"-like object
82  as done in GSL. (7/16/09 - I've tried this, and it doesn't
83  seem to improve the speed significantly.)
84  */
85  template<class mat_t>
86  int LU_decomp(const size_t N, mat_t &A, o2scl::permutation &p,
87  int &signum) {
88 
89  size_t i, j, k;
90 
91  signum=1;
92  p.init();
93 
94  for (j = 0; j < N - 1; j++) {
95 
96  /* Find maximum in the j-th column */
97  double ajj, max = fabs(O2SCL_IX2(A,j,j));
98  size_t i_pivot = j;
99 
100  for (i = j + 1; i < N; i++) {
101  double aij = fabs (O2SCL_IX2(A,i,j));
102 
103  if (aij > max) {
104  max = aij;
105  i_pivot = i;
106  }
107  }
108 
109  if (i_pivot != j) {
110 
111  // Swap rows j and i_pivot
112  double temp;
113  for (k=0;k<N;k++) {
114  temp=O2SCL_IX2(A,j,k);
115  O2SCL_IX2(A,j,k)=O2SCL_IX2(A,i_pivot,k);
116  O2SCL_IX2(A,i_pivot,k)=temp;
117  }
118  p.swap(j,i_pivot);
119  signum=-signum;
120  }
121 
122  ajj = O2SCL_IX2(A,j,j);
123 
124  if (ajj != 0.0) {
125  for (i = j + 1; i < N; i++) {
126  double aij = O2SCL_IX2(A,i,j) / ajj;
127  O2SCL_IX2(A,i,j)=aij;
128  for (k = j + 1; k < N; k++) {
129  double aik = O2SCL_IX2(A,i,k);
130  double ajk = O2SCL_IX2(A,j,k);
131  O2SCL_IX2(A,i,k)=aik - aij * ajk;
132  }
133  }
134  }
135  }
136 
137  return o2scl::success;
138  }
139 
140  /** \brief Solve a linear system after LU decomposition in place
141 
142  These functions solve the square system A x = b in-place using
143  the LU decomposition of A into (LU,p). On input x should contain
144  the right-hand side b, which is replaced by the solution on
145  output.
146  */
147  template<class mat_t, class vec_t>
148  int LU_svx(const size_t N, const mat_t &LU,
149  const o2scl::permutation &p, vec_t &x) {
150 
151  if (diagonal_has_zero(N,LU)) {
152  O2SCL_ERR("LU matrix is singular in LU_svx().",
154  }
155 
156  /* Apply permutation to RHS */
157  p.apply(x);
158 
159  /* Solve for c using forward-substitution, L c = P b */
160  o2scl_cblas::dtrsv(o2scl_cblas::o2cblas_RowMajor,
161  o2scl_cblas::o2cblas_Lower,
162  o2scl_cblas::o2cblas_NoTrans,
163  o2scl_cblas::o2cblas_Unit,
164  N,N,LU,x);
165 
166  /* Perform back-substitution, U x = c */
167  o2scl_cblas::dtrsv(o2scl_cblas::o2cblas_RowMajor,
168  o2scl_cblas::o2cblas_Upper,
169  o2scl_cblas::o2cblas_NoTrans,
170  o2scl_cblas::o2cblas_NonUnit,
171  N,N,LU,x);
172 
173  return o2scl::success;
174  }
175 
176  /** \brief An alternate form of LU decomposition with
177  matrix row objects
178 
179  \comment
180  I was testing this out, but it doesn't seem much faster than
181  the original LU_decomp().
182  \comment
183  */
184  template<class mat_t, class mat_row_t>
185  int LU_decomp_alt(const size_t N, mat_t &A, o2scl::permutation &p,
186  int &signum) {
187 
188  size_t i, j, k;
189 
190  signum=1;
191  p.init();
192 
193  for (j = 0; j < N - 1; j++) {
194 
195  /* Find maximum in the j-th column */
196  double ajj, max = fabs(O2SCL_IX2(A,j,j));
197  size_t i_pivot = j;
198 
199  for (i = j + 1; i < N; i++) {
200  double aij = fabs (O2SCL_IX2(A,i,j));
201 
202  if (aij > max) {
203  max = aij;
204  i_pivot = i;
205  }
206  }
207 
208  if (i_pivot != j) {
209 
210  // Swap rows j and i_pivot
211  double temp;
212  mat_row_t r1(A,j);
213  mat_row_t r2(A,i_pivot);
214  for (k=0;k<N;k++) {
215  temp=O2SCL_IX(r1,k);
216  O2SCL_IX(r1,k)=O2SCL_IX(r2,k);
217  O2SCL_IX(r2,k)=temp;
218  }
219  p.swap(j,i_pivot);
220  signum=-signum;
221  }
222 
223  ajj = O2SCL_IX2(A,j,j);
224 
225  if (ajj != 0.0) {
226  for (i = j + 1; i < N; i++) {
227  double aij = O2SCL_IX2(A,i,j) / ajj;
228  O2SCL_IX2(A,i,j)=aij;
229  for (k = j + 1; k < N; k++) {
230  double aik = O2SCL_IX2(A,i,k);
231  double ajk = O2SCL_IX2(A,j,k);
232  O2SCL_IX2(A,i,k)=aik - aij * ajk;
233  }
234  }
235  }
236  }
237 
238  return o2scl::success;
239  }
240 
241  /** \brief Solve a linear system after LU decomposition
242 
243  This function solve the square system A x = b using the LU
244  decomposition of A into (LU, p) given by gsl_linalg_LU_decomp or
245  gsl_linalg_complex_LU_decomp.
246  */
247  template<class mat_t, class vec_t, class vec2_t>
248  int LU_solve(const size_t N, const mat_t &LU, const o2scl::permutation &p,
249  const vec_t &b, vec2_t &x) {
250 
251  if (diagonal_has_zero(N,LU)) {
252  O2SCL_ERR("LU matrix is singular in LU_solve().",
254  }
255 
256  /* Copy x <- b */
257  o2scl::vector_copy(N,b,x);
258 
259  /* Solve for x */
260  return LU_svx(N,LU,p,x);
261  }
262 
263  /** \brief Refine the solution of a linear system
264 
265  These functions apply an iterative improvement to x, the solution
266  of A x = b, using the LU decomposition of A into (LU,p). The
267  initial residual r = A x - b is also computed and stored in
268  residual.
269  */
270  template<class mat_t, class mat2_t, class vec_t, class vec2_t, class vec3_t>
271  int LU_refine(const size_t N, const mat_t &A, const mat2_t &LU,
272  const o2scl::permutation &p, const vec_t &b, vec2_t &x,
273  vec3_t &residual) {
274 
275  if (diagonal_has_zero(N,LU)) {
276  O2SCL_ERR("LU matrix is singular in LU_refine().",
278  }
279 
280  /* Compute residual, residual = (A * x - b) */
281  o2scl::vector_copy(N,b,residual);
282  o2scl_cblas::dgemv(o2scl_cblas::o2cblas_RowMajor,
283  o2scl_cblas::o2cblas_NoTrans,
284  N,N,1.0,A,x,-1.0,residual);
285 
286  /* Find correction, delta = - (A^-1) * residual, and apply it */
287 
288  int status=LU_svx(N,LU,p,residual);
289  o2scl_cblas::daxpy(-1.0,5,residual,x);
290 
291  return status;
292  }
293 
294  /** \brief Compute the inverse of a matrix from its LU decomposition
295 
296  These functions compute the inverse of a matrix A from its LU
297  decomposition (LU,p), storing the result in the matrix inverse.
298  The inverse is computed by solving the system A x = b for each
299  column of the identity matrix. It is preferable to avoid direct
300  use of the inverse whenever possible, as the linear solver
301  functions can obtain the same result more efficiently and
302  reliably.
303 
304  \future Could rewrite to avoid mat_col_t, (9/16/09 - However,
305  the function may be faster if mat_col_t is left in, so it's
306  unclear what's best.)
307  */
308  template<class mat_t, class mat2_t, class mat_col_t>
309  int LU_invert(const size_t N, const mat_t &LU,
310  const o2scl::permutation &p, mat2_t &inverse) {
311 
312  if (diagonal_has_zero(N,LU)) {
313  O2SCL_ERR("LU matrix is singular in LU_invert().",
315  }
316 
317  size_t i;
318 
319  int status=o2scl::success;
320 
321  // Set matrix 'inverse' to the identity
322  for(i=0;i<N;i++) {
323  for(size_t j=0;j<N;j++) {
324  if (i==j) O2SCL_IX2(inverse,i,j)=1.0;
325  else O2SCL_IX2(inverse,i,j)=0.0;
326  }
327  }
328 
329  for (i = 0; i < N; i++) {
330  mat_col_t c(inverse,i);
331  int status_i=LU_svx(N,LU,p,c);
332 
333  if (status_i) {
334  status = status_i;
335  }
336  }
337 
338  return status;
339  }
340 
341  /** \brief Compute the determinant of a matrix from its LU decomposition
342 
343  Compute the determinant of a matrix A from its LU decomposition,
344  LU. The determinant is computed as the product of the diagonal
345  elements of U and the sign of the row permutation signum.
346  */
347  template<class mat_t>
348  double LU_det(const size_t N, const mat_t &LU, int signum) {
349 
350  size_t i;
351 
352  double det=(double)signum;
353 
354  for (i=0;i<N;i++) {
355  det*=O2SCL_IX2(LU,i,i);
356  }
357 
358  return det;
359  }
360 
361  /** \brief Compute the logarithm of the absolute value of the
362  determinant of a matrix from its LU decomposition
363 
364  Compute the logarithm of the absolute value of the determinant of
365  a matrix A, \f$ \ln|\det(A)| \f$, from its LU decomposition, LU.
366  This function may be useful if the direct computation of the
367  determinant would overflow or underflow.
368  */
369  template<class mat_t>
370  double LU_lndet(const size_t N, const mat_t &LU) {
371 
372  size_t i;
373  double lndet = 0.0;
374 
375  for (i = 0; i < N; i++) {
376  lndet+=log(fabs(O2SCL_IX2(LU,i,i)));
377  }
378 
379  return lndet;
380  }
381 
382  /** \brief Compute the sign of the
383  determinant of a matrix from its LU decomposition
384 
385  Compute the sign or phase factor of the determinant of a matrix
386  A, \f$ \det(A)/|\det(A)| \f$, from its LU decomposition, LU.
387  */
388  template<class mat_t>
389  int LU_sgndet(const size_t N, const mat_t &LU, int signum) {
390 
391  size_t i;
392  int s=signum;
393 
394  for (i=0;i<N;i++) {
395  double u=O2SCL_IX2(LU,i,i);
396  if (u<0) {
397  s*=-1;
398  } else if (u == 0) {
399  s=0;
400  break;
401  }
402  }
403 
404  return s;
405  }
406 
407 #ifdef DOXYGEN
408 }
409 #endif
void dtrsv(const enum o2cblas_order order, const enum o2cblas_uplo Uplo, const enum o2cblas_transpose TransA, const enum o2cblas_diag Diag, const size_t M, const size_t N, const mat_t &A, vec_t &X)
Compute .
Definition: cblas_base.h:308
double LU_lndet(const size_t N, const mat_t &LU)
Compute the logarithm of the absolute value of the determinant of a matrix from its LU decomposition...
Definition: lu_base.h:370
int LU_svx(const size_t N, const mat_t &LU, const o2scl::permutation &p, vec_t &x)
Solve a linear system after LU decomposition in place.
Definition: lu_base.h:148
A class for representing permutations.
Definition: permutation.h:70
int LU_decomp(const size_t N, mat_t &A, o2scl::permutation &p, int &signum)
Compute the LU decomposition of the matrix A.
Definition: lu_base.h:86
double LU_det(const size_t N, const mat_t &LU, int signum)
Compute the determinant of a matrix from its LU decomposition.
Definition: lu_base.h:348
void dgemv(const enum o2cblas_order order, const enum o2cblas_transpose TransA, const size_t M, const size_t N, const double alpha, const mat_t &A, const vec_t &X, const double beta, vec2_t &Y)
Compute .
Definition: cblas_base.h:218
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:127
int LU_sgndet(const size_t N, const mat_t &LU, int signum)
Compute the sign of the determinant of a matrix from its LU decomposition.
Definition: lu_base.h:389
The namespace for linear algebra classes and functions.
Definition: bidiag.h:36
input domain error, e.g sqrt(-1)
Definition: err_hnd.h:53
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
int LU_decomp_alt(const size_t N, mat_t &A, o2scl::permutation &p, int &signum)
An alternate form of LU decomposition with matrix row objects.
Definition: lu_base.h:185
int init()
Initialize permutation to the identity.
Definition: permutation.h:203
void daxpy(const double alpha, const size_t N, const vec_t &X, vec2_t &Y)
Compute .
Definition: cblas_base.h:98
int LU_refine(const size_t N, const mat_t &A, const mat2_t &LU, const o2scl::permutation &p, const vec_t &b, vec2_t &x, vec3_t &residual)
Refine the solution of a linear system.
Definition: lu_base.h:271
int LU_invert(const size_t N, const mat_t &LU, const o2scl::permutation &p, mat2_t &inverse)
Compute the inverse of a matrix from its LU decomposition.
Definition: lu_base.h:309
int LU_solve(const size_t N, const mat_t &LU, const o2scl::permutation &p, const vec_t &b, vec2_t &x)
Solve a linear system after LU decomposition.
Definition: lu_base.h:248
int diagonal_has_zero(const size_t N, mat_t &A)
Return 1 if at least one of the elements in the diagonal is zero.
Definition: lu_base.h:54
int swap(const size_t i, const size_t j)
Swap two elements of a permutation.
Definition: permutation.h:249
Success.
Definition: err_hnd.h:47
int apply(vec_t &v) const
Apply the permutation to a vector.
Definition: permutation.h:290

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).