bidiag_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/bidiag.c
24  *
25  * Copyright (C) 2001, 2007 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 bidiag_base.h
43  \brief File defining bidiagonalization functions
44 */
45 #ifdef DOXYGEN
46 namespace o2scl_linalg {
47 #endif
48 
49  /** \brief Factor a matrix into bidiagonal form
50 
51  Factor matrix A of size <tt>(M,N)</tt> with \f$ M\geq N \f$ into
52  \f$ A = U B V^T \f$ where U and V are orthogonal and B is upper
53  bidiagonal.
54 
55  After the function call, the matrix \f$ B \f$ is stored the
56  diagonal and first superdiagonal of \c A. The matrices \f$ U \f$
57  and \f$ V \f$ are stored as packed sets of Householder
58  transformations in the lower and upper triangular parts of \c A,
59  respectively.
60 
61  Adapted from the GSL version which was based on algorithm 5.4.2
62  in \ref Golub96.
63  */
64  template<class mat_t, class vec_t, class vec2_t>
65  int bidiag_decomp(size_t M, size_t N, mat_t &A,
66  vec_t &tau_U, vec2_t &tau_V) {
67 
68  if (M<N) {
69  O2SCL_ERR2("Bidiagonal decomposition requires M>=N in ",
70  "bidiag_comp().",o2scl::exc_ebadlen);
71  }
72 
73  for (size_t i=0;i<N;i++) {
74 
75  // [GSL] Apply Householder transformation to current column
76  double tau_i=householder_transform_subcol(A,i,i,M);
77 
78  // [GSL] Apply the transformation to the remaining columns
79  if (i+1<N) {
80  householder_hm_subcol(A,i,i+1,M,N,A,i,i,tau_i);
81  }
82  O2SCL_IX(tau_U,i)=tau_i;
83 
84  // [GSL] Apply Householder transformation to current row
85 
86  if (i+1<N) {
87  double tau2_i=householder_transform_subrow(A,i,i+1,N);
88 
89  // [GSL] Apply the transformation to the remaining rows
90  if (i+1<M) {
91  householder_mh_subrow(A,i+1,i+1,M,N,A,i,i+1,tau2_i);
92  }
93  O2SCL_IX(tau_V,i)=tau2_i;
94  }
95 
96  }
97 
98  return o2scl::success;
99  }
100 
101  /** \brief Unpack a matrix \c A with the bidiagonal decomposition
102  and create matrices \c U, \c V, diagonal \c diag and
103  superdiagonal \c superdiag
104 
105  Given a matrix \c A of size <tt>(M,N)</tt> with \f$ M \geq N \f$
106  created by \ref bidiag_decomp(), this function creates the
107  matrix \c U of size <tt>(M,N)</tt>, the matrix \c V of size
108  <tt>(N,N)</tt>, the diagonal \c diag of size \c N and the
109  super-diagonal \c superdiag of size \c N-1.
110  */
111  template<class mat_t, class vec_t, class mat2_t, class vec2_t,
112  class mat3_t, class vec3_t, class vec4_t>
113  int bidiag_unpack(size_t M, size_t N, const mat_t &A,
114  const vec_t &tau_U, mat2_t &U, const vec2_t &tau_V,
115  mat3_t &V, vec3_t &diag, vec4_t &superdiag) {
116 
117  if (M<N) {
118  O2SCL_ERR2("Matrix A must have M >= N in ",
119  "bidiag_unpack().",o2scl::exc_ebadlen);
120  }
121 
122  const size_t K=N;
123 
124  // [GSL] Copy diagonal into diag
125 
126  for (size_t i=0;i<N;i++) {
127  O2SCL_IX(diag,i)=O2SCL_IX2(A,i,i);
128  }
129 
130  // [GSL] Copy superdiagonal into superdiag
131 
132  for (size_t i=0;i<N-1;i++) {
133  O2SCL_IX(superdiag,i)=O2SCL_IX2(A,i,i+1);
134  }
135 
136  // [GSL] Initialize V to the identity
137  for(size_t i=0;i<N;i++) {
138  for(size_t j=0;j<N;j++) {
139  if (i==j) O2SCL_IX2(V,i,j)=1.0;
140  else O2SCL_IX2(V,i,j)=0.0;
141  }
142  }
143 
144  for (size_t i=N-1;i-- > 0;) {
145 
146  // [GSL] Householder row transformation to accumulate V
147  householder_hm_subrow(V,i+1,i+1,N,N,A,i,i+1,O2SCL_IX(tau_V,i));
148  }
149 
150  // [GSL] Initialize U to the identity
151 
152  for(size_t i=0;i<M;i++) {
153  for(size_t j=0;j<N;j++) {
154  if (i==j) O2SCL_IX2(U,i,j)=1.0;
155  else O2SCL_IX2(U,i,j)=0.0;
156  }
157  }
158 
159  for (size_t j=N;j-- > 0;) {
160  householder_hm_subcol(U,j,j,M,N,A,j,j,O2SCL_IX(tau_U,j));
161  }
162 
163  return o2scl::success;
164  }
165 
166  /** \brief Unpack a matrix \c A with the bidiagonal decomposition
167  and create matrix \c V
168  */
169  template<class mat_t, class vec_t, class vec2_t, class mat2_t>
170  int bidiag_unpack2(size_t M, size_t N, mat_t &A, vec_t &tau_U,
171  vec2_t &tau_V, mat2_t &V) {
172 
173  if (M<N) {
174  O2SCL_ERR2("Matrix A must have M >= N in ",
175  "bidiag_unpack2().",o2scl::exc_ebadlen);
176  }
177 
178  const size_t K=M;
179 
180  // [GSL] Initialize V to the identity
181 
182  for(size_t i=0;i<N;i++) {
183  for(size_t j=0;j<N;j++) {
184  if (i==j) O2SCL_IX2(V,i,j)=1.0;
185  else O2SCL_IX2(V,i,j)=0.0;
186  }
187  }
188 
189  for (size_t i=N-1;i-- > 0;) {
190 
191  // [GSL] Householder row transformation to accumulate V
192  householder_hm_subrow(V,i+1,i+1,N,N,A,i,i+1,O2SCL_IX(tau_V,i));
193 
194  }
195 
196  // [GSL] Copy superdiagonal into tau_v
197 
198  for (size_t i=0;i<N-1;i++) {
199  O2SCL_IX(tau_V,i)=O2SCL_IX2(A,i,i+1);
200  }
201 
202  // [GSL] Allow U to be unpacked into the same memory as A, copy
203  // diagonal into tau_U
204 
205  for (size_t j=N; j-- > 0;) {
206  // [GSL] Householder column transformation to accumulate U
207  double tj=O2SCL_IX(tau_U,j);
208  O2SCL_IX(tau_U,j)=O2SCL_IX2(A,j,j);
209  householder_hm1_sub(M,N,tj,A,j,j);
210  }
211 
212  return o2scl::success;
213  }
214 
215  /** \brief Unpack the diagonal and superdiagonal of the bidiagonal
216  decomposition of \c A into \c diag and \c superdiag
217  */
218  template<class mat_t, class vec_t, class vec2_t>
219  int bidiag_unpack_B(size_t M, size_t N, const mat_t &A,
220  vec_t &diag, vec2_t &superdiag) {
221 
222  size_t K=N;
223  if (M<=N) K=M;
224 
225  // [GSL] Copy diagonal into diag
226  for (size_t i=0;i<K;i++) {
227  O2SCL_IX(diag,i)=O2SCL_IX2(A,i,i);
228  }
229 
230  // [GSL] Copy superdiagonal into superdiag
231  for (size_t i=0;i<K-1;i++) {
232  O2SCL_IX(superdiag,i)=O2SCL_IX2(A,i,i+1);
233  }
234 
235  return o2scl::success;
236  }
237 
238 #ifdef DOXYGEN
239 }
240 #endif
241 
int bidiag_unpack2(size_t M, size_t N, mat_t &A, vec_t &tau_U, vec2_t &tau_V, mat2_t &V)
Unpack a matrix A with the bidiagonal decomposition and create matrix V.
Definition: bidiag_base.h:170
int bidiag_unpack(size_t M, size_t N, const mat_t &A, const vec_t &tau_U, mat2_t &U, const vec2_t &tau_V, mat3_t &V, vec3_t &diag, vec4_t &superdiag)
Unpack a matrix A with the bidiagonal decomposition and create matrices U, V, diagonal diag and super...
Definition: bidiag_base.h:113
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 householder_hm_subrow(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 ...
double householder_transform_subrow(mat_t &A, const size_t ir, const size_t ic, const size_t N)
Compute the Householder transform of a vector formed with the last n columns of a row of a matrix...
matrix, vector lengths are not conformant
Definition: err_hnd.h:89
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 ...
#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
int bidiag_unpack_B(size_t M, size_t N, const mat_t &A, vec_t &diag, vec2_t &superdiag)
Unpack the diagonal and superdiagonal of the bidiagonal decomposition of A into diag and superdiag...
Definition: bidiag_base.h:219
int bidiag_decomp(size_t M, size_t N, mat_t &A, vec_t &tau_U, vec2_t &tau_V)
Factor a matrix into bidiagonal form.
Definition: bidiag_base.h:65
void householder_mh_subrow(mat_t &M, const size_t ir, const size_t ic, const size_t nr, const size_t nc, const mat2_t &M2, const size_t ir2, const size_t ic2, double tau)
Apply the Householder transformation (v,tau) to the right-hand side of the matrix A...
Success.
Definition: err_hnd.h:47
void householder_hm1_sub(const size_t M, const size_t N, double tau, mat_t &A, size_t ir, size_t ic)
Apply a Householder transformation to a matrix being build up from the identity matrix, using the first column of A as a Householder vector.

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