cholesky_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 /* 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.
60  */
61  template<class mat_t> void cholesky_decomp(const size_t M, mat_t &A) {
62 
63  size_t i,j,k;
64 
65  /* [GSL] Do the first 2 rows explicitly. It is simple, and faster.
66  And one can return if the matrix has only 1 or 2 rows.
67  */
68 
69  double A_00=O2SCL_IX2(A,0,0);
70 
71  // AWS: The GSL version stores GSL_NAN in L_00 and then throws
72  // an error if A_00 <= 0. We throw the error first and then
73  // the square root should always be safe?
74 
75  if (A_00<=0) {
76  O2SCL_ERR2("Matrix not positive definite (A[0][0]<=0) in ",
77  "cholesky_decomp().",o2scl::exc_einval);
78  }
79 
80  double L_00=sqrt(A_00);
81  O2SCL_IX2(A,0,0)=L_00;
82 
83  if (M>1) {
84  double A_10=O2SCL_IX2(A,1,0);
85  double A_11=O2SCL_IX2(A,1,1);
86 
87  double L_10=A_10/L_00;
88  double diag=A_11-L_10*L_10;
89 
90  if (diag<=0) {
91  O2SCL_ERR2("Matrix not positive definite (diag<=0 for 2x2) in ",
92  "cholesky_decomp().",o2scl::exc_einval);
93  }
94  double L_11=sqrt(diag);
95 
96  O2SCL_IX2(A,1,0)=L_10;
97  O2SCL_IX2(A,1,1)=L_11;
98  }
99 
100  for (k=2;k<M;k++) {
101  double A_kk=O2SCL_IX2(A,k,k);
102 
103  for (i=0;i<k;i++) {
104  double sum=0.0;
105 
106  double A_ki=O2SCL_IX2(A,k,i);
107  double A_ii=O2SCL_IX2(A,i,i);
108 
109  // AWS: Should change to use a form of ddot() here
110  if (i>0) {
111  sum=0.0;
112  for(j=0;j<i;j++) {
113  sum+=O2SCL_IX2(A,i,j)*O2SCL_IX2(A,k,j);
114  }
115  }
116 
117  A_ki=(A_ki-sum)/A_ii;
118  O2SCL_IX2(A,k,i)=A_ki;
119  }
120 
121  {
122  double sum=dnrm2_subrow(A,k,0,k);
123  double diag=A_kk-sum*sum;
124 
125  if(diag<=0) {
126  O2SCL_ERR2("Matrix not positive definite (diag<=0) in ",
127  "cholesky_decomp().",o2scl::exc_einval);
128  }
129 
130  double L_kk=sqrt(diag);
131 
132  O2SCL_IX2(A,k,k)=L_kk;
133  }
134  }
135 
136  /* [GSL] Now copy the transposed lower triangle to the upper
137  triangle, the diagonal is common.
138  */
139 
140  for (i=1;i<M;i++) {
141  for (j=0;j<i;j++) {
142  double A_ij=O2SCL_IX2(A,i,j);
143  O2SCL_IX2(A,j,i)=A_ij;
144  }
145  }
146 
147  return;
148  }
149 
150  /** \brief Solve a symmetric positive-definite linear system after a
151  Cholesky decomposition
152 
153  Given the Cholesky decomposition of a matrix A in \c LLT,
154  this function solves the system <tt>A*x=b</tt>.
155  */
156  template<class mat_t, class vec_t, class vec2_t>
157  void cholesky_solve(const size_t N, const mat_t &LLT,
158  const vec_t &b, vec2_t &x) {
159 
160  // [GSL] Copy x <- b
161  o2scl::vector_copy(N,b,x);
162 
163  // [GSL] Solve for c using forward-substitution, L c=b
164  o2scl_cblas::dtrsv(o2scl_cblas::o2cblas_RowMajor,
165  o2scl_cblas::o2cblas_Lower,
166  o2scl_cblas::o2cblas_NoTrans,
167  o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
168 
169  // [GSL] Perform back-substitution,U x=c
170  o2scl_cblas::dtrsv(o2scl_cblas::o2cblas_RowMajor,
171  o2scl_cblas::o2cblas_Upper,
172  o2scl_cblas::o2cblas_NoTrans,
173  o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
174 
175  return;
176  }
177 
178  /** \brief Solve a linear system in place using a Cholesky decomposition
179  */
180  template<class mat_t, class vec_t>
181  void cholesky_svx(const size_t N, const mat_t &LLT, vec_t &x) {
182 
183  // [GSL] Solve for c using forward-substitution, L c=b
184  o2scl_cblas::dtrsv(o2scl_cblas::o2cblas_RowMajor,
185  o2scl_cblas::o2cblas_Lower,
186  o2scl_cblas::o2cblas_NoTrans,
187  o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
188 
189  // [GSL] Perform back-substitution, U x=c
190  o2scl_cblas::dtrsv(o2scl_cblas::o2cblas_RowMajor,
191  o2scl_cblas::o2cblas_Upper,
192  o2scl_cblas::o2cblas_NoTrans,
193  o2scl_cblas::o2cblas_NonUnit,N,N,LLT,x);
194 
195  return;
196  }
197 
198  /** \brief Compute the inverse of a symmetric positive definite matrix
199  given the Cholesky decomposition
200 
201  Given a Cholesky decomposition produced by \ref cholesky_decomp(),
202  this function returns the inverse of that matrix in \c LLT.
203  */
204  template<class mat_t> void cholesky_invert(const size_t N, mat_t &LLT) {
205 
206  size_t i, j;
207  double sum;
208 
209  // [GSL] invert the lower triangle of LLT
210  for (i=0;i<N;++i) {
211 
212  j=N-i-1;
213 
214  O2SCL_IX2(LLT,j,j)=1.0/O2SCL_IX2(LLT,j,j);
215  double ajj=-O2SCL_IX2(LLT,j,j);
216 
217  if (j<N-1) {
218 
219  // This section is just the equivalent of dtrmv() for a part of
220  // the matrix LLT.
221  {
222 
223  size_t ix=N-j-2;
224  for (size_t ii=N-j-1;ii>0 && ii--;) {
225  double temp=0.0;
226  const size_t j_min=0;
227  const size_t j_max=ii;
228  size_t jx=j_min;
229  for (size_t jj=j_min;jj<j_max;jj++) {
230  temp+=O2SCL_IX2(LLT,jx+j+1,j)*O2SCL_IX2(LLT,ii+j+1,jj+j+1);
231  jx++;
232  }
233  O2SCL_IX2(LLT,ix+j+1,j)=temp+O2SCL_IX2(LLT,ix+j+1,j)*
234  O2SCL_IX2(LLT,ii+j+1,ii+j+1);
235  ix--;
236  }
237 
238  }
239 
240  o2scl_cblas::dscal_subcol(LLT,j+1,j,N,ajj);
241 
242  }
243  }
244 
245  /*
246  [GSL] The lower triangle of LLT now contains L^{-1}. Now compute
247  A^{-1}=L^{-t} L^{-1}
248 
249  The (ij) element of A^{-1} is column i of L^{-1} dotted into
250  column j of L^{-1}
251  */
252 
253  for (i=0;i<N;++i) {
254 
255  for (j=i+1;j<N;++j) {
256 
257  // [GSL] Compute Ainv_{ij}=sum_k Linv_{ki} Linv_{kj}.
258 
259  // AWS: Should change to use a form of ddot() here
260  sum=0.0;
261  for(size_t k=j;k<N;k++) {
262  sum+=O2SCL_IX2(LLT,k,i)*O2SCL_IX2(LLT,k,j);
263  }
264 
265  // [GSL] Store in upper triangle
266  O2SCL_IX2(LLT,i,j)=sum;
267  }
268 
269  // [GSL] now compute the diagonal element
270 
271  // AWS: Should change to use a form of ddot() here
272  sum=0.0;
273  for(size_t k=i;k<N;k++) {
274  sum+=O2SCL_IX2(LLT,k,i)*O2SCL_IX2(LLT,k,i);
275  }
276 
277  O2SCL_IX2(LLT,i,i)=sum;
278  }
279 
280  // [GSL] Copy the transposed upper triangle to the lower triangle
281 
282  for (j=1;j<N;j++) {
283  for (i=0;i<j;i++) {
284  O2SCL_IX2(LLT,j,i)=O2SCL_IX2(LLT,i,j);
285  }
286  }
287 
288  return;
289  }
290 
291  /** \brief Cholesky decomposition with unit-diagonal triangular parts.
292 
293  A = L D L^T, where diag(L) = (1,1,...,1).
294  Upon exit, A contains L and L^T as for Cholesky, and
295  the diagonal of A is (1,1,...,1). The vector Dis set
296  to the diagonal elements of the diagonal matrix D.
297  */
298  template<class mat_t, class vec_t>
299  int cholesky_decomp_unit(const size_t N, mat_t &A, vec_t &D) {
300 
301  size_t i, j;
302 
303  // [GSL] Initial Cholesky decomposition
304  int stat_chol=cholesky_decomp(N,A);
305 
306  // [GSL] Calculate D from diagonal part of initial Cholesky
307  for(i=0;i<N;++i) {
308  const double C_ii=O2SCL_IX2(A,i,i);
309  O2SCL_IX(D,i)=C_ii*C_ii;
310  }
311 
312  // [GSL] Multiply initial Cholesky by 1/sqrt(D) on the right
313  for(i=0;i<N;++i) {
314  for(j=0;j<N;++j) {
315  O2SCL_IX2(A,i,j)=O2SCL_IX2(A,i,j)/sqrt(O2SCL_IX(D,j));
316  }
317  }
318 
319  /* [GSL] Because the initial Cholesky contained both L and
320  transpose(L), the result of the multiplication is not symmetric
321  anymore; but the lower triangle _is_ correct. Therefore we
322  reflect it to the upper triangle and declare victory.
323  */
324  for(i=0;i<N;++i) {
325  for(j=i+1;j<N;++j) {
326  O2SCL_IX2(A,i,j)=O2SCL_IX2(A,j,i);
327  }
328  }
329 
330  return stat_chol;
331  }
332 
333 #ifdef DOXYGEN
334 }
335 #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
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
void cholesky_decomp(const size_t M, mat_t &A)
Compute the in-place Cholesky decomposition of a symmetric positive-definite square matrix...
Definition: cholesky_base.h:61
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
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:1262
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:1093

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