cholesky_base.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2018, 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 /* Cholesky Decomposition
24  *
25  * Copyright (C) 2000 Thomas Walter
26  *
27  * 03 May 2000: Modified for GSL by Brian Gough
28  * 29 Jul 2005: Additions by Gerard Jungman
29  * Copyright (C) 2000,2001, 2002, 2003, 2005, 2007 Brian Gough,
30  * Gerard Jungman
31  *
32  * This is free software; you can redistribute it and/or modify it
33  * under the terms of the GNU General Public License as published by the
34  * Free Software Foundation; either version 3, or (at your option) any
35  * later version.
36  *
37  * This source is distributed in the hope that it will be useful, but
38  * WITHOUT ANY WARRANTY; without even the implied warranty of
39  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
40  * General Public License for more details.
41  */
42 /** \file cholesky_base.h
43  \brief File defining Cholesky decomposition
44 */
45 
46 #ifdef DOXYGEN
47 namespace o2scl_linalg {
48 #endif
49 
50  /** \brief Compute the in-place Cholesky decomposition of a symmetric
51  positive-definite square matrix
52 
53  On input, the upper triangular part of A is ignored (only the
54  lower triangular part and diagonal are used). On output,
55  the diagonal and lower triangular part contain the matrix L
56  and the upper triangular part contains L^T.
57 
58  If the matrix is not positive-definite, the error handler
59  will be called, unless \c err_on_fail is false, in which
60  case a non-zero value will be returned.
61  */
62  template<class mat_t> int cholesky_decomp(const size_t M, mat_t &A,
63  bool err_on_fail=true) {
64 
65  size_t i,j,k;
66 
67  /* [GSL] Do the first 2 rows explicitly. It is simple, and faster.
68  And one can return if the matrix has only 1 or 2 rows.
69  */
70 
71  double A_00=O2SCL_IX2(A,0,0);
72 
73  // AWS: The GSL version stores GSL_NAN in L_00 and then throws
74  // an error if A_00 <= 0. We throw the error first and then
75  // the square root should always be safe?
76 
77  if (A_00<=0.0) {
78  if (err_on_fail) {
79  O2SCL_ERR2("Matrix not positive definite (A[0][0]<=0) in ",
80  "cholesky_decomp().",o2scl::exc_einval);
81  } else {
82  return 1;
83  }
84  }
85 
86  double L_00=sqrt(A_00);
87  O2SCL_IX2(A,0,0)=L_00;
88 
89  if (M>1) {
90  double A_10=O2SCL_IX2(A,1,0);
91  double A_11=O2SCL_IX2(A,1,1);
92 
93  double L_10=A_10/L_00;
94  double diag=A_11-L_10*L_10;
95 
96  if (diag<=0.0) {
97  if (err_on_fail) {
98  O2SCL_ERR2("Matrix not positive definite (diag<=0 for 2x2) in ",
99  "cholesky_decomp().",o2scl::exc_einval);
100  } else {
101  return 2;
102  }
103  }
104  double L_11=sqrt(diag);
105 
106  O2SCL_IX2(A,1,0)=L_10;
107  O2SCL_IX2(A,1,1)=L_11;
108  }
109 
110  for (k=2;k<M;k++) {
111  double A_kk=O2SCL_IX2(A,k,k);
112 
113  for (i=0;i<k;i++) {
114  double sum=0.0;
115 
116  double A_ki=O2SCL_IX2(A,k,i);
117  double A_ii=O2SCL_IX2(A,i,i);
118 
119  // AWS: Should change to use a form of ddot() here
120  if (i>0) {
121  sum=0.0;
122  for(j=0;j<i;j++) {
123  sum+=O2SCL_IX2(A,i,j)*O2SCL_IX2(A,k,j);
124  }
125  }
126 
127  A_ki=(A_ki-sum)/A_ii;
128  O2SCL_IX2(A,k,i)=A_ki;
129  }
130 
131  {
132  double sum=dnrm2_subrow(A,k,0,k);
133  double diag=A_kk-sum*sum;
134 
135  if (diag<=0.0) {
136  if (err_on_fail) {
137  O2SCL_ERR2("Matrix not positive definite (diag<=0) in ",
138  "cholesky_decomp().",o2scl::exc_einval);
139  } else {
140  return 3;
141  }
142  }
143 
144  double L_kk=sqrt(diag);
145 
146  O2SCL_IX2(A,k,k)=L_kk;
147  }
148  }
149 
150  /* [GSL] Now copy the transposed lower triangle to the upper
151  triangle, the diagonal is common.
152  */
153 
154  for (i=1;i<M;i++) {
155  for (j=0;j<i;j++) {
156  double A_ij=O2SCL_IX2(A,i,j);
157  O2SCL_IX2(A,j,i)=A_ij;
158  }
159  }
160 
161  return 0;
162  }
163 
164  /** \brief Solve a symmetric positive-definite linear system after a
165  Cholesky decomposition
166 
167  Given the Cholesky decomposition of a matrix A in \c LLT,
168  this function solves the system <tt>A*x=b</tt>.
169  */
170  template<class mat_t, class vec_t, class vec2_t>
171  void cholesky_solve(const size_t N, const mat_t &LLT,
172  const vec_t &b, vec2_t &x) {
173 
174  // [GSL] Copy x <- b
175  o2scl::vector_copy(N,b,x);
176 
177  // [GSL] Solve for c using forward-substitution, L c=b
178  o2scl_cblas::dtrsv(o2scl_cblas::o2cblas_RowMajor,
179  o2scl_cblas::o2cblas_Lower,
180  o2scl_cblas::o2cblas_NoTrans,
181  o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
182 
183  // [GSL] Perform back-substitution,U x=c
184  o2scl_cblas::dtrsv(o2scl_cblas::o2cblas_RowMajor,
185  o2scl_cblas::o2cblas_Upper,
186  o2scl_cblas::o2cblas_NoTrans,
187  o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
188 
189  return;
190  }
191 
192  /** \brief Solve a linear system in place using a Cholesky decomposition
193  */
194  template<class mat_t, class vec_t>
195  void cholesky_svx(const size_t N, const mat_t &LLT, vec_t &x) {
196 
197  // [GSL] Solve for c using forward-substitution, L c=b
198  o2scl_cblas::dtrsv(o2scl_cblas::o2cblas_RowMajor,
199  o2scl_cblas::o2cblas_Lower,
200  o2scl_cblas::o2cblas_NoTrans,
201  o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
202 
203  // [GSL] Perform back-substitution, U x=c
204  o2scl_cblas::dtrsv(o2scl_cblas::o2cblas_RowMajor,
205  o2scl_cblas::o2cblas_Upper,
206  o2scl_cblas::o2cblas_NoTrans,
207  o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
208 
209  return;
210  }
211 
212  /** \brief Compute the inverse of a symmetric positive definite matrix
213  given the Cholesky decomposition
214 
215  Given a Cholesky decomposition produced by \ref cholesky_decomp(),
216  this function returns the inverse of that matrix in \c LLT.
217  */
218  template<class mat_t> void cholesky_invert(const size_t N, mat_t &LLT) {
219 
220  size_t i, j;
221  double sum;
222 
223  // [GSL] invert the lower triangle of LLT
224  for (i=0;i<N;++i) {
225 
226  j=N-i-1;
227 
228  O2SCL_IX2(LLT,j,j)=1.0/O2SCL_IX2(LLT,j,j);
229  double ajj=-O2SCL_IX2(LLT,j,j);
230 
231  if (j<N-1) {
232 
233  // This section is just the equivalent of dtrmv() for a part of
234  // the matrix LLT.
235  {
236 
237  size_t ix=N-j-2;
238  for (size_t ii=N-j-1;ii>0 && ii--;) {
239  double temp=0.0;
240  const size_t j_min=0;
241  const size_t j_max=ii;
242  size_t jx=j_min;
243  for (size_t jj=j_min;jj<j_max;jj++) {
244  temp+=O2SCL_IX2(LLT,jx+j+1,j)*O2SCL_IX2(LLT,ii+j+1,jj+j+1);
245  jx++;
246  }
247  O2SCL_IX2(LLT,ix+j+1,j)=temp+O2SCL_IX2(LLT,ix+j+1,j)*
248  O2SCL_IX2(LLT,ii+j+1,ii+j+1);
249  ix--;
250  }
251 
252  }
253 
254  o2scl_cblas::dscal_subcol(LLT,j+1,j,N,ajj);
255 
256  }
257  }
258 
259  /*
260  [GSL] The lower triangle of LLT now contains L^{-1}. Now compute
261  A^{-1}=L^{-t} L^{-1}
262 
263  The (ij) element of A^{-1} is column i of L^{-1} dotted into
264  column j of L^{-1}
265  */
266 
267  for (i=0;i<N;++i) {
268 
269  for (j=i+1;j<N;++j) {
270 
271  // [GSL] Compute Ainv_{ij}=sum_k Linv_{ki} Linv_{kj}.
272 
273  // AWS: Should change to use a form of ddot() here
274  sum=0.0;
275  for(size_t k=j;k<N;k++) {
276  sum+=O2SCL_IX2(LLT,k,i)*O2SCL_IX2(LLT,k,j);
277  }
278 
279  // [GSL] Store in upper triangle
280  O2SCL_IX2(LLT,i,j)=sum;
281  }
282 
283  // [GSL] now compute the diagonal element
284 
285  // AWS: Should change to use a form of ddot() here
286  sum=0.0;
287  for(size_t k=i;k<N;k++) {
288  sum+=O2SCL_IX2(LLT,k,i)*O2SCL_IX2(LLT,k,i);
289  }
290 
291  O2SCL_IX2(LLT,i,i)=sum;
292  }
293 
294  // [GSL] Copy the transposed upper triangle to the lower triangle
295 
296  for (j=1;j<N;j++) {
297  for (i=0;i<j;i++) {
298  O2SCL_IX2(LLT,j,i)=O2SCL_IX2(LLT,i,j);
299  }
300  }
301 
302  return;
303  }
304 
305  /** \brief Cholesky decomposition with unit-diagonal triangular parts.
306 
307  A = L D L^T, where diag(L) = (1,1,...,1).
308  Upon exit, A contains L and L^T as for Cholesky, and
309  the diagonal of A is (1,1,...,1). The vector Dis set
310  to the diagonal elements of the diagonal matrix D.
311  */
312  template<class mat_t, class vec_t>
313  int cholesky_decomp_unit(const size_t N, mat_t &A, vec_t &D) {
314 
315  size_t i, j;
316 
317  // [GSL] Initial Cholesky decomposition
318  int stat_chol=cholesky_decomp(N,A);
319 
320  // [GSL] Calculate D from diagonal part of initial Cholesky
321  for(i=0;i<N;++i) {
322  const double C_ii=O2SCL_IX2(A,i,i);
323  O2SCL_IX(D,i)=C_ii*C_ii;
324  }
325 
326  // [GSL] Multiply initial Cholesky by 1/sqrt(D) on the right
327  for(i=0;i<N;++i) {
328  for(j=0;j<N;++j) {
329  O2SCL_IX2(A,i,j)=O2SCL_IX2(A,i,j)/sqrt(O2SCL_IX(D,j));
330  }
331  }
332 
333  /* [GSL] Because the initial Cholesky contained both L and
334  transpose(L), the result of the multiplication is not symmetric
335  anymore; but the lower triangle _is_ correct. Therefore we
336  reflect it to the upper triangle and declare victory.
337  */
338  for(i=0;i<N;++i) {
339  for(j=i+1;j<N;++j) {
340  O2SCL_IX2(A,i,j)=O2SCL_IX2(A,j,i);
341  }
342  }
343 
344  return stat_chol;
345  }
346 
347 #ifdef DOXYGEN
348 }
349 #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:335
invalid argument supplied by user
Definition: err_hnd.h:59
void cholesky_svx(const size_t N, const mat_t &LLT, vec_t &x)
Solve a linear system in place using a Cholesky decomposition.
int cholesky_decomp_unit(const size_t N, mat_t &A, vec_t &D)
Cholesky decomposition with unit-diagonal triangular parts.
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:127
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
int cholesky_decomp(const size_t M, mat_t &A, bool err_on_fail=true)
Compute the in-place Cholesky decomposition of a symmetric positive-definite square matrix...
Definition: cholesky_base.h:62
The namespace for linear algebra classes and functions.
Definition: bidiag.h:36
double dnrm2_subrow(const mat_t &M, const size_t ir, const size_t ic, const size_t N)
Compute the norm of a subrow of a matrix.
Definition: cblas_base.h:1289
void cholesky_solve(const size_t N, const mat_t &LLT, const vec_t &b, vec2_t &x)
Solve a symmetric positive-definite linear system after a Cholesky decomposition. ...
void cholesky_invert(const size_t N, mat_t &LLT)
Compute the inverse of a symmetric positive definite matrix given the Cholesky decomposition.
void dscal_subcol(mat_t &A, const size_t ir, const size_t ic, const size_t M, const double alpha)
Compute for a subcolumn of a matrix.
Definition: cblas_base.h:1120

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