lu_base.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, 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 #ifdef O2SCL_NEVER_DEFINED
141 
142  template<class mat_t, class vec_size_t>
143  int LU_decomp_L3_sub(const size_t M, const size_t N, mat_t &A,
144  vec_size_t &ipiv, size_t istart,
145  size_t jstart) {
146 
147  if (M < N) {
148  O2SCL_ERR("matrix must have M >= N",o2scl::exc_ebadlen);
149  } else if (N <= 24) {
150  /* use Level 2 algorithm */
151  return LU_decomp_L2(A, ipiv);
152  } else {
153  /*
154  * partition matrix:
155  *
156  * N1 N2
157  * N1 [ A11 A12 ]
158  * M2 [ A21 A22 ]
159  *
160  * and
161  * N1 N2
162  * M [ AL AR ]
163  */
164  int status;
165  const size_t N1=((N >= 16) ? ((N + 8) / 16) * 8 : N / 2);
166  const size_t N2 = N - N1;
167  const size_t M2 = M - N1;
168  /*
169  gsl_matrix_view A11 = gsl_matrix_submatrix(A, 0, 0, N1, N1);
170  gsl_matrix_view A12 = gsl_matrix_submatrix(A, 0, N1, N1, N2);
171  gsl_matrix_view A21 = gsl_matrix_submatrix(A, N1, 0, M2, N1);
172  gsl_matrix_view A22 = gsl_matrix_submatrix(A, N1, N1, M2, N2);
173 
174  gsl_matrix_view AL = gsl_matrix_submatrix(A, 0, 0, M, N1);
175  gsl_matrix_view AR = gsl_matrix_submatrix(A, 0, N1, M, N2);
176  */
177 
178  /*
179  * partition ipiv = [ ipiv1 ] N1
180  * [ ipiv2 ] N2
181  */
182  //gsl_vector_uint_view ipiv1 = gsl_vector_uint_subvector(ipiv, 0, N1);
183  //gsl_vector_uint_view ipiv2 = gsl_vector_uint_subvector(ipiv, N1, N2);
184 
185  size_t i;
186 
187  /* recursion on (AL, ipiv1) */
188  status=LU_decomp_L3(M,N1,AL,ipiv,0,0);
189  if (status) {
190  return status;
191  }
192 
193  /* apply ipiv1 to AR */
194  apply_pivots(&AR.matrix, &ipiv1.vector);
195 
196  /* A12 = A11^{-1} A12 */
197  dtrsm_submat(CblasLeft,CblasLower,CblasNoTrans,CblasUnit,
198  N1,N2,1.0,0,0,0,N1);
199  //gsl_blas_dtrsm(CblasLeft, CblasLower, CblasNoTrans, CblasUnit,
200  //1.0, &A11.matrix, &A12.matrix);
201 
202  /* A22 = A22 - A21 * A12 */
203  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, -1.0, &A21.matrix, &A12.matrix, 1.0, &A22.matrix);
204 
205  /* recursion on (A22, ipiv2) */
206  status = LU_decomp_L3(&A22.matrix, &ipiv2.vector);
207  if (status)
208  return status;
209 
210  /* apply pivots to A21 */
211  apply_pivots(&A21.matrix, &ipiv2.vector);
212 
213  /* shift pivots */
214  for (i = 0; i < N2; ++i)
215  {
216  unsigned int * ptr = gsl_vector_uint_ptr(&ipiv2.vector, i);
217  *ptr += N1;
218  }
219 
220  return GSL_SUCCESS;
221  }
222 
223  return 0;
224  }
225 #endif
226 
227  /** \brief Solve a linear system after LU decomposition in place
228 
229  These functions solve the square system A x = b in-place using
230  the LU decomposition of A into (LU,p). On input x should contain
231  the right-hand side b, which is replaced by the solution on
232  output.
233  */
234  template<class mat_t, class vec_t>
235  int LU_svx(const size_t N, const mat_t &LU,
236  const o2scl::permutation &p, vec_t &x) {
237 
238  if (diagonal_has_zero(N,LU)) {
239  O2SCL_ERR("LU matrix is singular in LU_svx().",
241  }
242 
243  /* Apply permutation to RHS */
244  p.apply(x);
245 
246  /* Solve for c using forward-substitution, L c = P b */
247  o2scl_cblas::dtrsv(o2scl_cblas::o2cblas_RowMajor,
248  o2scl_cblas::o2cblas_Lower,
249  o2scl_cblas::o2cblas_NoTrans,
250  o2scl_cblas::o2cblas_Unit,
251  N,N,LU,x);
252 
253  /* Perform back-substitution, U x = c */
254  o2scl_cblas::dtrsv(o2scl_cblas::o2cblas_RowMajor,
255  o2scl_cblas::o2cblas_Upper,
256  o2scl_cblas::o2cblas_NoTrans,
257  o2scl_cblas::o2cblas_NonUnit,
258  N,N,LU,x);
259 
260  return o2scl::success;
261  }
262 
263  /** \brief An alternate form of LU decomposition with
264  matrix row objects
265 
266  \comment
267  I was testing this out, but it doesn't seem much faster than
268  the original LU_decomp().
269  \comment
270  */
271  template<class mat_t, class mat_row_t>
272  int LU_decomp_alt(const size_t N, mat_t &A, o2scl::permutation &p,
273  int &signum) {
274 
275  size_t i, j, k;
276 
277  signum=1;
278  p.init();
279 
280  for (j = 0; j < N - 1; j++) {
281 
282  /* Find maximum in the j-th column */
283  double ajj, max = fabs(O2SCL_IX2(A,j,j));
284  size_t i_pivot = j;
285 
286  for (i = j + 1; i < N; i++) {
287  double aij = fabs (O2SCL_IX2(A,i,j));
288 
289  if (aij > max) {
290  max = aij;
291  i_pivot = i;
292  }
293  }
294 
295  if (i_pivot != j) {
296 
297  // Swap rows j and i_pivot
298  double temp;
299  mat_row_t r1(A,j);
300  mat_row_t r2(A,i_pivot);
301  for (k=0;k<N;k++) {
302  temp=O2SCL_IX(r1,k);
303  O2SCL_IX(r1,k)=O2SCL_IX(r2,k);
304  O2SCL_IX(r2,k)=temp;
305  }
306  p.swap(j,i_pivot);
307  signum=-signum;
308  }
309 
310  ajj = O2SCL_IX2(A,j,j);
311 
312  if (ajj != 0.0) {
313  for (i = j + 1; i < N; i++) {
314  double aij = O2SCL_IX2(A,i,j) / ajj;
315  O2SCL_IX2(A,i,j)=aij;
316  for (k = j + 1; k < N; k++) {
317  double aik = O2SCL_IX2(A,i,k);
318  double ajk = O2SCL_IX2(A,j,k);
319  O2SCL_IX2(A,i,k)=aik - aij * ajk;
320  }
321  }
322  }
323  }
324 
325  return o2scl::success;
326  }
327 
328  /** \brief Solve a linear system after LU decomposition
329 
330  This function solve the square system A x = b using the LU
331  decomposition of A into (LU, p) given by gsl_linalg_LU_decomp or
332  gsl_linalg_complex_LU_decomp.
333  */
334  template<class mat_t, class vec_t, class vec2_t>
335  int LU_solve(const size_t N, const mat_t &LU, const o2scl::permutation &p,
336  const vec_t &b, vec2_t &x) {
337 
338  if (diagonal_has_zero(N,LU)) {
339  O2SCL_ERR("LU matrix is singular in LU_solve().",
341  }
342 
343  /* Copy x <- b */
344  o2scl::vector_copy(N,b,x);
345 
346  /* Solve for x */
347  return LU_svx(N,LU,p,x);
348  }
349 
350  /** \brief Refine the solution of a linear system
351 
352  These functions apply an iterative improvement to x, the solution
353  of A x = b, using the LU decomposition of A into (LU,p). The
354  initial residual r = A x - b is also computed and stored in
355  residual.
356  */
357  template<class mat_t, class mat2_t, class vec_t, class vec2_t, class vec3_t>
358  int LU_refine(const size_t N, const mat_t &A, const mat2_t &LU,
359  const o2scl::permutation &p, const vec_t &b, vec2_t &x,
360  vec3_t &residual) {
361 
362  if (diagonal_has_zero(N,LU)) {
363  O2SCL_ERR("LU matrix is singular in LU_refine().",
365  }
366 
367  /* Compute residual, residual = (A * x - b) */
368  o2scl::vector_copy(N,b,residual);
369  o2scl_cblas::dgemv(o2scl_cblas::o2cblas_RowMajor,
370  o2scl_cblas::o2cblas_NoTrans,
371  N,N,1.0,A,x,-1.0,residual);
372 
373  /* Find correction, delta = - (A^-1) * residual, and apply it */
374 
375  int status=LU_svx(N,LU,p,residual);
376  o2scl_cblas::daxpy(-1.0,5,residual,x);
377 
378  return status;
379  }
380 
381  /** \brief Compute the inverse of a matrix from its LU decomposition
382 
383  These functions compute the inverse of a matrix A from its LU
384  decomposition (LU,p), storing the result in the matrix \c
385  inverse. The inverse is computed by solving the system A x = b
386  for each column of the identity matrix. It is preferable to
387  avoid direct use of the inverse whenever possible, as the linear
388  solver functions can obtain the same result more efficiently and
389  reliably.
390 
391  \future Could rewrite to avoid mat_col_t, (9/16/09 - However,
392  the function may be faster if mat_col_t is left in, so it's
393  unclear what's best.)
394  */
395  template<class mat_t, class mat2_t, class mat_col_t>
396  int LU_invert(const size_t N, const mat_t &LU,
397  const o2scl::permutation &p, mat2_t &inverse) {
398 
399  if (diagonal_has_zero(N,LU)) {
400  O2SCL_ERR("LU matrix is singular in LU_invert().",
402  }
403 
404  size_t i;
405 
406  int status=o2scl::success;
407 
408  // Set matrix 'inverse' to the identity
409  for(i=0;i<N;i++) {
410  for(size_t j=0;j<N;j++) {
411  if (i==j) O2SCL_IX2(inverse,i,j)=1.0;
412  else O2SCL_IX2(inverse,i,j)=0.0;
413  }
414  }
415 
416  for (i = 0; i < N; i++) {
417  mat_col_t c(inverse,i);
418  int status_i=LU_svx(N,LU,p,c);
419 
420  if (status_i) {
421  status = status_i;
422  }
423  }
424 
425  return status;
426  }
427 
428  /** \brief Compute the determinant of a matrix from its LU decomposition
429 
430  Compute the determinant of a matrix A from its LU decomposition,
431  LU. The determinant is computed as the product of the diagonal
432  elements of U and the sign of the row permutation signum.
433  */
434  template<class mat_t>
435  double LU_det(const size_t N, const mat_t &LU, int signum) {
436 
437  size_t i;
438 
439  double det=(double)signum;
440 
441  for (i=0;i<N;i++) {
442  det*=O2SCL_IX2(LU,i,i);
443  }
444 
445  return det;
446  }
447 
448  /** \brief Compute the logarithm of the absolute value of the
449  determinant of a matrix from its LU decomposition
450 
451  Compute the logarithm of the absolute value of the determinant of
452  a matrix A, \f$ \ln|\det(A)| \f$, from its LU decomposition, LU.
453  This function may be useful if the direct computation of the
454  determinant would overflow or underflow.
455  */
456  template<class mat_t>
457  double LU_lndet(const size_t N, const mat_t &LU) {
458 
459  size_t i;
460  double lndet = 0.0;
461 
462  for (i = 0; i < N; i++) {
463  lndet+=log(fabs(O2SCL_IX2(LU,i,i)));
464  }
465 
466  return lndet;
467  }
468 
469  /** \brief Compute the sign of the
470  determinant of a matrix from its LU decomposition
471 
472  Compute the sign or phase factor of the determinant of a matrix
473  A, \f$ \det(A)/|\det(A)| \f$, from its LU decomposition, LU.
474  */
475  template<class mat_t>
476  int LU_sgndet(const size_t N, const mat_t &LU, int signum) {
477 
478  size_t i;
479  int s=signum;
480 
481  for (i=0;i<N;i++) {
482  double u=O2SCL_IX2(LU,i,i);
483  if (u<0) {
484  s*=-1;
485  } else if (u == 0) {
486  s=0;
487  break;
488  }
489  }
490 
491  return s;
492  }
493 
494 #ifdef DOXYGEN
495 }
496 #endif
o2scl_cblas::dgemv
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:245
o2scl::permutation::apply
int apply(vec_t &v) const
Apply the permutation to a vector.
Definition: permutation.h:290
o2scl::permutation::swap
int swap(const size_t i, const size_t j)
Swap two elements of a permutation.
Definition: permutation.h:249
o2scl_cblas::daxpy
void daxpy(const double alpha, const size_t N, const vec_t &X, vec2_t &Y)
Compute .
Definition: cblas_base.h:125
o2scl_linalg::LU_svx
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:235
o2scl_linalg
The namespace for linear algebra classes and functions.
Definition: bidiag.h:36
o2scl::permutation
A class for representing permutations.
Definition: permutation.h:70
o2scl::permutation::init
int init()
Initialize permutation to the identity.
Definition: permutation.h:203
o2scl_linalg::LU_lndet
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:457
o2scl_linalg::LU_decomp
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
o2scl::vector_copy
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:124
o2scl::exc_edom
@ exc_edom
input domain error, e.g sqrt(-1)
Definition: err_hnd.h:53
o2scl::success
@ success
Success.
Definition: err_hnd.h:47
o2scl_linalg::LU_sgndet
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:476
o2scl_linalg::LU_det
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:435
o2scl_linalg::LU_decomp_alt
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:272
o2scl_linalg::LU_invert
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:396
o2scl::exc_ebadlen
@ exc_ebadlen
matrix, vector lengths are not conformant
Definition: err_hnd.h:89
o2scl_cblas::dtrsv
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:335
O2SCL_ERR
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
o2scl_linalg::LU_refine
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:358
o2scl_linalg::LU_solve
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:335
o2scl_linalg::diagonal_has_zero
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

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