qr_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/qr.c
24  *
25  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman,
26  * Brian Gough
27  *
28  * This program is free software; you can redistribute it and/or modify
29  * it under the terms of the GNU General Public License as published by
30  * the Free Software Foundation; either version 3 of the License, or (at
31  * your option) any later version.
32  *
33  * This program is distributed in the hope that it will be useful, but
34  * WITHOUT ANY WARRANTY; without even the implied warranty of
35  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
36  * General Public License for more details.
37  *
38  * You should have received a copy of the GNU General Public License
39  * along with this program; if not, write to the Free Software
40  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
41  * 02110-1301, USA.
42  */
43 /** \file qr_base.h
44  \brief File for QR decomposition and associated solver
45 */
46 
47 #ifdef DOXYGEN
48 namespace o2scl_linalg {
49 #endif
50 
51  /** \brief Compute the QR decomposition of matrix \c A
52  */
53  template<class mat_t, class vec_t>
54  void QR_decomp(size_t M, size_t N, mat_t &A, vec_t &tau) {
55 
56  size_t imax;
57  if (M<N) imax=M;
58  else imax=N;
59  for(size_t i=0;i<imax;i++) {
60  O2SCL_IX(tau,i)=householder_transform_subcol(A,i,i,M);
61  if (i+1<N) {
62  householder_hm_subcol(A,i,i+1,M,N,A,i,i,O2SCL_IX(tau,i));
63  }
64  }
65  return;
66  }
67 
68  /** \brief Form the product Q^T v from a QR factorized matrix
69  */
70  template<class mat_t, class vec_t, class vec2_t>
71  void QR_QTvec(const size_t M, const size_t N,
72  const mat_t &QR, const vec_t &tau, vec2_t &v) {
73 
74  // compute Q^T v
75  size_t imax;
76  if (M<N) imax=M;
77  else imax=N;
78  for(size_t i=0;i<imax;i++) {
79  householder_hv_subcol(QR,v,O2SCL_IX(tau,i),i,M);
80  }
81  return;
82  }
83 
84  /** \brief Unpack the QR matrix to the individual Q and R components
85  */
86  template<class mat1_t, class mat2_t, class mat3_t, class vec_t>
87  void QR_unpack(const size_t M, const size_t N,
88  const mat1_t &QR, const vec_t &tau, mat2_t &Q,
89  mat3_t &R) {
90 
91  size_t i, j;
92 
93  // [GSL] Initialize Q to the identity
94 
95  for(i=0;i<M;i++) {
96  for(j=0;j<M;j++) {
97  if (i==j) O2SCL_IX2(Q,i,j)=1.0;
98  else O2SCL_IX2(Q,i,j)=0.0;
99  }
100  }
101 
102  size_t istart;
103  if (M<N) istart=M;
104  else istart=N;
105 
106  for (i=istart; i-- > 0;) {
107  householder_hm_subcol(Q,i,i,M,M,QR,i,i,O2SCL_IX(tau,i));
108  }
109 
110  // [GSL] Form the right triangular matrix R from a packed QR matrix
111 
112  for (i=0; i < M; i++) {
113  for (j=0; j < i && j < N; j++) {
114  O2SCL_IX2(R,i,j)=0.0;
115  }
116  for (j=i; j < N; j++) {
117  O2SCL_IX2(R,i,j)=O2SCL_IX2(QR,i,j);
118  }
119  }
120 
121  return;
122  }
123 
124  /** \brief Solve the system A x = b in place using the QR factorization
125  */
126  template<class mat_t, class vec_t, class vec2_t>
127  void QR_svx(size_t M, size_t N, const mat_t &QR, const vec_t &tau,
128  vec2_t &x) {
129  QR_QTvec(M,N,QR,tau,x);
130  o2scl_cblas::dtrsv(o2scl_cblas::o2cblas_RowMajor,
131  o2scl_cblas::o2cblas_Upper,
132  o2scl_cblas::o2cblas_NoTrans,
133  o2scl_cblas::o2cblas_NonUnit,
134  M,N,QR,x);
135  return;
136  }
137 
138  /** \brief Solve the system A x = b using the QR factorization
139  */
140  template<class mat_t, class vec_t, class vec2_t, class vec3_t>
141  void QR_solve(size_t N, const mat_t &QR, const vec_t &tau,
142  const vec2_t &b, vec3_t &x) {
143  o2scl::vector_copy(N,b,x);
144  QR_svx(N,N,QR,tau,x);
145  return;
146  }
147 
148  /** \brief Update a QR factorisation for A= Q R, A' = A + u v^T,
149 
150  The parameters \c M and \c N are the number of rows and columns
151  of the matrix \c R.
152 
153  \verbatim
154  * Q' R' = QR + u v^T
155  * = Q (R + Q^T u v^T)
156  * = Q (R + w v^T)
157  *
158  * where w = Q^T u.
159  *
160  * Algorithm from Golub and Van Loan, "Matrix Computations", Section
161  * 12.5 (Updating Matrix Factorizations, Rank-One Changes)
162  \endverbatim
163  */
164  template<class mat1_t, class mat2_t, class vec1_t, class vec2_t>
165  void QR_update(size_t M, size_t N, mat1_t &Q, mat2_t &R,
166  vec1_t &w, vec2_t &v) {
167 
168  size_t j;
169  // Integer to avoid problems with decreasing loop below
170  int k;
171  double w0;
172 
173  /* [GSL] Apply Given's rotations to reduce w to (|w|, 0, 0, ... , 0)
174 
175  J_1^T .... J_(n-1)^T w = +/- |w| e_1
176 
177  simultaneously applied to R, H = J_1^T ... J^T_(n-1) R
178  so that H is upper Hessenberg. (12.5.2)
179  */
180 
181  // Loop from k = M-1 to 1
182  for (k=(int)(M-1); k > 0; k--) {
183 
184  double c, s;
185  double wk=O2SCL_IX(w,k);
186  double wkm1=O2SCL_IX(w,k-1);
187 
188  o2scl_linalg::create_givens(wkm1,wk,c,s);
189  apply_givens_vec(w,k-1,k,c,s);
190  apply_givens_qr(M,N,Q,R,k-1,k,c,s);
191  }
192 
193  w0=O2SCL_IX(w,0);
194 
195  // [GSL] Add in w v^T (Equation 12.5.3)
196 
197  for (j=0; j < N; j++) {
198  double r0j=O2SCL_IX2(R,0,j);
199  double vj=O2SCL_IX(v,j);
200  O2SCL_IX2(R,0,j)=r0j+w0*vj;
201  }
202 
203  // [GSL] Apply Givens transformations R' = G_(n-1)^T ... G_1^T H
204  // (Equation 12.5.4)
205  int kmax;
206  if (M<N+1) kmax=M;
207  else kmax=N+1;
208  for (k=1;k<kmax;k++) {
209 
210  double c, s;
211  double diag=O2SCL_IX2(R,k-1,k-1);
212  double offdiag=O2SCL_IX2(R,k,k-1);
213 
214  o2scl_linalg::create_givens(diag,offdiag,c,s);
215  apply_givens_qr(M,N,Q,R,k-1,k,c,s);
216 
217  O2SCL_IX2(R,k,k-1)=0.0;
218  }
219 
220  return;
221  }
222 
223  /** \brief Compute the unpacked QR decomposition of matrix \c A
224 
225  If \o2 is compiled with Armadillo support, this is specialized
226  for <tt>arma::mat</tt> to use <tt>arma::qr_econ</tt>. If \o2 is
227  compiled with Eigen support, this is specialized for
228  <tt>Eigen::MatrixXd</tt>.
229  */
230  template<class mat_t, class mat2_t, class mat3_t>
231  void QR_decomp_unpack(const size_t M, const size_t N,
232  mat_t &A, mat2_t &Q, mat3_t &R) {
233 
235  QR_decomp(M,N,A,tau);
236  QR_unpack(M,N,A,tau,Q,R);
237 
238  return;
239  }
240 
241 #ifdef DOXYGEN
242 }
243 #endif
void QR_solve(size_t N, const mat_t &QR, const vec_t &tau, const vec2_t &b, vec3_t &x)
Solve the system A x = b using the QR factorization.
Definition: qr_base.h:141
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
void QR_unpack(const size_t M, const size_t N, const mat1_t &QR, const vec_t &tau, mat2_t &Q, mat3_t &R)
Unpack the QR matrix to the individual Q and R components.
Definition: qr_base.h:87
void apply_givens_vec(vec_t &v, size_t i, size_t j, double c, double s)
Apply a rotation to a vector, .
Definition: givens_base.h:110
void QR_decomp(size_t M, size_t N, mat_t &A, vec_t &tau)
Compute the QR decomposition of matrix A.
Definition: qr_base.h:54
void QR_decomp_unpack(const size_t M, const size_t N, mat_t &A, mat2_t &Q, mat3_t &R)
Compute the unpacked QR decomposition of matrix A.
Definition: qr_base.h:231
double householder_transform_subcol(mat_t &A, const size_t ir, const size_t ic, const size_t M)
Compute the Householder transform of a vector formed with n rows of a column of a matrix...
void create_givens(const double a, const double b, double &c, double &s)
Create a Givens rotation matrix.
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:127
void householder_hm_subcol(mat_t &M, const size_t ir, const size_t ic, const size_t nr, const size_t nc, const mat_t &M2, const size_t ir2, const size_t ic2, double tau)
Apply a Householder transformation to the lower-right part of M when the transformation is stored in ...
The namespace for linear algebra classes and functions.
Definition: bidiag.h:36
void QR_update(size_t M, size_t N, mat1_t &Q, mat2_t &R, vec1_t &w, vec2_t &v)
Update a QR factorisation for A= Q R, A&#39; = A + u v^T,.
Definition: qr_base.h:165
void apply_givens_qr(size_t M, size_t N, mat1_t &Q, mat2_t &R, size_t i, size_t j, double c, double s)
Apply a rotation to matrices from the QR decomposition.
Definition: givens_base.h:58
void QR_QTvec(const size_t M, const size_t N, const mat_t &QR, const vec_t &tau, vec2_t &v)
Form the product Q^T v from a QR factorized matrix.
Definition: qr_base.h:71
void QR_svx(size_t M, size_t N, const mat_t &QR, const vec_t &tau, vec2_t &x)
Solve the system A x = b in place using the QR factorization.
Definition: qr_base.h:127
void householder_hv_subcol(const mat_t &A, vec_t &w, double tau, const size_t ie, const size_t N)
Apply a Householder transformation v to vector w where v is stored as a column in a matrix A...

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