householder_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/householder.c
24  *
25  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004,
26  * 2007 Gerard Jungman, 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 householder_base.h
44  \brief File for Householder transformations
45 
46  \todo Better documentation for the Householder functions.
47 */
48 
49 #include <gsl/gsl_machine.h>
50 
51 #include <o2scl/err_hnd.h>
52 #include <o2scl/cblas.h>
53 #include <o2scl/permutation.h>
54 
55 #ifdef DOXYGEN
56 namespace o2scl_linalg {
57 #endif
58 
59  /** \brief Replace the vector \c v with a Householder vector and a
60  coefficient tau that annihilates <tt>v[1]</tt> through
61  <tt>v[n-1]</tt> (inclusive)
62 
63  On exit, this function returns the value of \f$ \tau = 2/ (v^{T}
64  v) \f$. If \c n is less than or equal to 1 then this function
65  returns zero without calling the error handler.
66  */
67  template<class vec_t>
68  double householder_transform(const size_t n, vec_t &v) {
69 
70  if (n <= 1) {
71 
72  // tau=0
73  return 0.0;
74 
75  } else {
76 
77  double alpha, beta, tau;
78  double xnorm=O2SCL_CBLAS_NAMESPACE::dnrm2_subvec(n,v,1);
79 
80  if (xnorm == 0) {
81  // tau=0
82  return 0.0;
83  }
84 
85  alpha=O2SCL_IX(v,0);
86  beta=-(alpha >= 0.0 ? +1.0 : -1.0)*hypot(alpha,xnorm);
87  tau=(beta-alpha)/beta;
88 
89  double dbl_eps=std::numeric_limits<double>::epsilon();
90  double dbl_min=std::numeric_limits<double>::min();
91 
92  double s=(alpha-beta);
93  if (fabs(s)>dbl_min) {
95  O2SCL_IX(v,0)=beta;
96  } else {
97  O2SCL_CBLAS_NAMESPACE::dscal_subvec(dbl_eps/s,n,v,1);
98  O2SCL_CBLAS_NAMESPACE::dscal_subvec(1.0/dbl_eps,n,v,1);
99  O2SCL_IX(v,0)=beta;
100  }
101 
102  return tau;
103  }
104  }
105 
106  /** \brief Compute the Householder transform of a vector
107  formed with \c n rows of a column of a matrix
108 
109  This performs a Householder transform of a vector defined by a
110  column of a matrix \c A which starts at element
111  <tt>A(ir,ic)</tt> and ends at element <tt>A(M-1,ic)</tt>.
112  If <tt>M-1</tt> is equal to <tt>ir+1</tt>, this function quietly
113  does nothing.
114 
115  Used in \ref QR_decomp() and \ref SV_decomp_mod().
116  */
117  template<class mat_t>
118  double householder_transform_subcol(mat_t &A, const size_t ir,
119  const size_t ic, const size_t M) {
120 
121  if (M == ir+1) {
122  // [GSL] tau=0
123  return 0.0;
124  } else {
125  double alpha, beta, tau;
126 
127  double xnorm=O2SCL_CBLAS_NAMESPACE::dnrm2_subcol(A,ir+1,ic,M);
128 
129  if (xnorm == 0) {
130  // [GSL] tau=0
131  return 0.0;
132  }
133 
134  alpha=O2SCL_IX2(A,ir,ic);
135  beta=-(alpha >= 0.0 ? +1.0 : -1.0)*hypot(alpha,xnorm);
136  tau=(beta-alpha)/beta;
137 
138  double dbl_eps=std::numeric_limits<double>::epsilon();
139  double dbl_min=std::numeric_limits<double>::min();
140 
141  double s=(alpha-beta);
142  if (fabs(s)>dbl_min) {
143  O2SCL_CBLAS_NAMESPACE::dscal_subcol(A,ir+1,ic,M,1.0/s);
144  O2SCL_IX2(A,ir,ic)=beta;
145  } else {
146  O2SCL_CBLAS_NAMESPACE::dscal_subcol(A,ir+1,ic,M,dbl_eps/s);
147  O2SCL_CBLAS_NAMESPACE::dscal_subcol(A,ir+1,ic,M,1.0/dbl_eps);
148  O2SCL_IX2(A,ir,ic)=beta;
149  }
150 
151  return tau;
152  }
153  }
154 
155  /** \brief Compute the Householder transform of a vector
156  formed with the last \c n columns of a row of a matrix
157 
158  This performs a Householder transform of a vector
159  defined by a row of a matrix \c A which starts at element
160  <tt>A(ir,ic)</tt> and ends at element <tt>A(ir,N-1)</tt>
161  If <tt>N-1</tt> is equal to <tt>ic</tt>, this function quietly
162  does nothing.
163  */
164  template<class mat_t>
165  double householder_transform_subrow(mat_t &A, const size_t ir,
166  const size_t ic, const size_t N) {
167 
168  if (N == ic+1) {
169  // [GSL] tau=0
170  return 0.0;
171  } else {
172  double alpha, beta, tau;
173 
174  double xnorm=O2SCL_CBLAS_NAMESPACE::dnrm2_subrow(A,ir,ic+1,N);
175 
176  if (xnorm == 0) {
177  // [GSL] tau=0
178  return 0.0;
179  }
180 
181  alpha=O2SCL_IX2(A,ir,ic);
182  beta=-(alpha >= 0.0 ? +1.0 : -1.0)*hypot(alpha,xnorm);
183  tau=(beta-alpha)/beta;
184 
185  double dbl_eps=std::numeric_limits<double>::epsilon();
186  double dbl_min=std::numeric_limits<double>::min();
187 
188  double s=(alpha-beta);
189  if (fabs(s)>dbl_min) {
190  O2SCL_CBLAS_NAMESPACE::dscal_subrow(A,ir,ic+1,N,1.0/s);
191  O2SCL_IX2(A,ir,ic)=beta;
192  } else {
193  O2SCL_CBLAS_NAMESPACE::dscal_subrow(A,ir,ic+1,N,dbl_eps/s);
194  O2SCL_CBLAS_NAMESPACE::dscal_subrow(A,ir,ic+1,N,1.0/dbl_eps);
195  O2SCL_IX2(A,ir,ic)=beta;
196  }
197 
198  return tau;
199  }
200  }
201 
202  /** \brief Apply a Householder transformation \f$ (v,\tau) \f$ to
203  matrix \f$ A \f$ of size <tt>(M,N)</tt>
204 
205  The vector \c v must have at least \c N entries, with the
206  exception that the vector element <tt>v[0]</tt> is never
207  referenced by this function.
208  */
209  template<class vec_t, class mat_t>
210  void householder_hm(const size_t M, const size_t N,
211  double tau, const vec_t &v, mat_t &A) {
212 
213  if (tau == 0.0) {
214  return;
215  }
216 
217  for (size_t j=0;j<N;j++) {
218 
219  // [GSL] Compute wj=Akj vk
220  double wj=O2SCL_IX2(A,0,j);
221 
222  for (size_t i=1;i<M;i++) {
223  wj+=O2SCL_IX2(A,i,j)*O2SCL_IX(v,i);
224  }
225 
226  // [GSL] Aij=Aij-tau vi wj
227 
228  // [GSL] i=0
229  O2SCL_IX2(A,0,j)-=tau*wj;
230 
231  // [GSL] i=1 .. M-1
232  for (size_t i=1;i<M;i++) {
233  O2SCL_IX2(A,i,j)-=tau*wj*O2SCL_IX(v,i);
234  }
235  }
236  return;
237  }
238 
239  /** \brief Apply a Householder transformation to the lower-right
240  part of \c M when the transformation is stored in a column of
241  \c M2
242 
243  This applies a householder transformation <tt>(v,tau)</tt> to a
244  lower-right submatrix of \c M. The submatrix has <tt>nr-ir</tt>
245  rows and <tt>nc-ic</tt> columns and starts at row \c ir of
246  column \c ic of the original matrix \c M. The vector containing
247  the transformation is taken from a column of \c M2 starting at
248  row \c ir2 and column \c ic2. The matrix \c M2 must have at
249  least <tt>ic2+1</tt> columns and at least <tt>nr-ir+ir2</tt>
250  rows.
251 
252  This function is used in \ref QR_decomp() and \ref QR_unpack() .
253  */
254  template<class mat_t>
255  void householder_hm_subcol(mat_t &M, const size_t ir,
256  const size_t ic, const size_t nr,
257  const size_t nc, const mat_t &M2,
258  const size_t ir2, const size_t ic2,
259  double tau) {
260 
261  if (tau == 0.0) {
262  return;
263  }
264 
265  for (size_t j=ic;j<nc;j++) {
266 
267  // [GSL] Compute wj=Akj vk
268  double wj=O2SCL_IX2(M,ir,j);
269 
270  for (size_t i=ir+1;i<nr;i++) {
271  wj+=O2SCL_IX2(M,i,j)*O2SCL_IX2(M2,i-ir+ir2,ic2);
272  }
273 
274  // [GSL] Aij=Aij-tau vi wj
275 
276  // [GSL] i=0
277  O2SCL_IX2(M,ir,j)-=tau*wj;
278 
279  // [GSL] i=1 .. M-1
280  for (size_t i=ir+1;i<nr;i++) {
281  O2SCL_IX2(M,i,j)-=tau*wj*O2SCL_IX2(M2,i-ir+ir2,ic2);
282  }
283  }
284  return;
285  }
286 
287  /** \brief Apply a Householder transformation to the lower-right
288  part of \c M when the transformation is stored in a row of \c M2
289 
290  This applies a householder transformation <tt>(v,tau)</tt> to a
291  lower-right submatrix of \c M. The submatrix has <tt>nr-ir</tt>
292  rows and <tt>nc-ic</tt> columns and starts at row \c ir of
293  column \c ic of the original matrix \c M. The vector containing
294  the transformation is taken from a row of \c M2 starting at row
295  \c ir2 and column \c ic2. The matrix \c M2 must have
296  at least <tt>ir2+1</tt> rows and <tt>nr-ir+ic2</tt> columns.
297 
298  Used in \ref bidiag_unpack().
299  */
300  template<class mat_t>
301  void householder_hm_subrow(mat_t &M, const size_t ir,
302  const size_t ic, const size_t nr,
303  const size_t nc, const mat_t &M2,
304  const size_t ir2, const size_t ic2,
305  double tau) {
306 
307  if (tau == 0.0) {
308  return;
309  }
310 
311  for (size_t j=ic;j<nc;j++) {
312 
313  // [GSL] Compute wj=Akj vk
314  double wj=O2SCL_IX2(M,ir,j);
315 
316  for (size_t i=ir+1;i<nr;i++) {
317  wj+=O2SCL_IX2(M,i,j)*O2SCL_IX2(M2,ir2,i-ir+ic2);
318  }
319 
320  // [GSL] Aij=Aij-tau vi wj
321 
322  // [GSL] i=0
323  O2SCL_IX2(M,ir,j)-=tau*wj;
324 
325  // [GSL] i=1 .. M-1
326  for (size_t i=ir+1;i<nr;i++) {
327  O2SCL_IX2(M,i,j)-=tau*wj*O2SCL_IX2(M2,ir2,i-ic+ic2);
328  }
329  }
330  return;
331  }
332 
333  /** \brief Apply a Householder transformation \c v to vector \c w
334  */
335  template<class vec_t, class vec2_t>
336  void householder_hv(const size_t N, double tau, const vec_t &v,
337  vec2_t &w) {
338 
339  if (tau==0) return;
340 
341  // compute d=v'w
342 
343  double d0=O2SCL_IX(w,0);
344  double d1, d;
345 
347 
348  d=d0+d1;
349 
350  // compute w=w-tau(v)(v'w)
351  O2SCL_IX(w,0)-=tau*d;
352 
353  O2SCL_CBLAS_NAMESPACE::daxpy_subvec(-tau*d,N,v,w,1);
354 
355  return;
356  }
357 
358  /** \brief Apply a Householder transformation \c v to vector \c w
359  where \c v is stored as a column in a matrix \c A
360 
361  Used in \ref QR_QTvec().
362  */
363  template<class mat_t, class vec_t>
364  void householder_hv_subcol(const mat_t &A, vec_t &w, double tau,
365  const size_t ie, const size_t N) {
366 
367  if (tau==0) return;
368 
369  // compute d=v'w
370 
371  double d0=O2SCL_IX(w,ie);
372  double d1, d;
373 
374  d1=O2SCL_CBLAS_NAMESPACE::ddot_subcol(N,A,ie+1,ie,w);
375 
376  d=d0+d1;
377 
378  // compute w=w-tau(v)(v'w)
379  O2SCL_IX(w,ie)-=tau*d;
380 
381  O2SCL_CBLAS_NAMESPACE::daxpy_subcol(-tau*d,N,A,ie+1,ie,w);
382 
383  return;
384  }
385 
386  /** \brief Apply a Householder transformation \f$ (v,\tau) \f$ to a
387  matrix being build up from the identity matrix, using the first
388  column of A as a Householder vector
389  */
390  template<class mat_t>
391  void householder_hm1(const size_t M, const size_t N,
392  double tau, mat_t &A) {
393 
394  if (tau == 0) {
395  O2SCL_IX2(A,0,0)=1.0;
396  for (size_t j=1;j<N;j++) {
397  O2SCL_IX2(A,0,j)=0.0;
398  }
399  for (size_t i=1;i<M;i++) {
400  O2SCL_IX2(A,i,0)=0.0;
401  }
402  return;
403  }
404 
405  // [GSL] w=A' v
406 
407  for (size_t j=1;j<N;j++) {
408  double wj=0.0;
409  for (size_t i=1;i<M;i++) {
410  wj+=O2SCL_IX2(A,i,j)*O2SCL_IX2(A,i,0);
411  }
412 
413  // [GSL] A=A-tau v w'
414  O2SCL_IX2(A,0,j)=-tau*wj;
415 
416  for (size_t i=1;i<M;i++) {
417  O2SCL_IX2(A,i,j)-=tau*wj*O2SCL_IX2(A,i,0);
418  }
419  }
420 
421  for (size_t i=1;i<M;i++) {
422  O2SCL_IX2(A,i,0)=-tau*O2SCL_IX2(A,i,0);
423  }
424  O2SCL_IX2(A,0,0)=1.0-tau;
425 
426  return;
427  }
428 
429  /** \brief Apply a Householder transformation \f$ (v,\tau) \f$ to a
430  matrix being build up from the identity matrix, using the first
431  column of A as a Householder vector
432 
433  Used in \ref SV_decomp_mod() and \ref bidiag_unpack2().
434  */
435  template<class mat_t>
436  void householder_hm1_sub(const size_t M, const size_t N,
437  double tau, mat_t &A, size_t ir, size_t ic) {
438 
439  size_t irp1=ir+1;
440  size_t icp1=ic+1;
441 
442  if (tau == 0) {
443  O2SCL_IX2(A,ir,ic)=1.0;
444  for (size_t j=icp1;j<N;j++) {
445  O2SCL_IX2(A,ir,j)=0.0;
446  }
447  for (size_t i=irp1;i<M;i++) {
448  O2SCL_IX2(A,i,ic)=0.0;
449  }
450  return;
451  }
452 
453  // [GSL] w=A' v
454 
455  for (size_t j=icp1;j<N;j++) {
456  double wj=0.0;
457  for (size_t i=irp1;i<M;i++) {
458  wj+=O2SCL_IX2(A,i,j)*O2SCL_IX2(A,i,ic);
459  }
460 
461  // [GSL] A=A-tau v w'
462  O2SCL_IX2(A,ir,j)=-tau*wj;
463 
464  for (size_t i=irp1;i<M;i++) {
465  O2SCL_IX2(A,i,j)-=tau*wj*O2SCL_IX2(A,i,ic);
466  }
467  }
468 
469  for (size_t i=irp1;i<M;i++) {
470  O2SCL_IX2(A,i,ic)=-tau*O2SCL_IX2(A,i,ic);
471  }
472  O2SCL_IX2(A,ir,ic)=1.0-tau;
473 
474  return;
475  }
476 
477  /** \brief Apply the Householder transformation <tt>(v,tau)</tt> to
478  the right-hand side of the matrix \c A.
479  */
480  template<class vec_t, class mat_t>
481  void householder_mh(const size_t M, const size_t N,
482  double tau, const vec_t &v, mat_t &A) {
483 
484  if (tau==0.0) {
485  return;
486  }
487 
488  // [GSL] A=A-tau w v'
489  size_t i, j;
490 
491  for (i=0;i<M;i++) {
492 
493  double wi=O2SCL_IX2(A,i,0);
494 
495  // [GSL] note, computed for v(0)=1 above
496  for (j=1;j<N;j++) {
497  wi+=O2SCL_IX2(A,i,j)*O2SCL_IX(v,j);
498  }
499 
500  // [GSL] j=0
501  O2SCL_IX2(A,i,0)-=tau*wi;
502 
503  // [GSL] j=1 .. N-1
504  for (j=1;j<N;j++) {
505  O2SCL_IX2(A,i,j)-=tau*wi*O2SCL_IX(v,j);
506  }
507  }
508 
509  return;
510  }
511 
512  /** \brief Apply the Householder transformation <tt>(v,tau)</tt> to
513  the right-hand side of the matrix \c A.
514 
515  Used in \ref bidiag_decomp().
516  */
517  template<class mat_t, class mat2_t>
519  (mat_t &M, const size_t ir, const size_t ic,
520  const size_t nr, const size_t nc, const mat2_t &M2,
521  const size_t ir2, const size_t ic2, double tau) {
522 
523  if (tau==0.0) {
524  return;
525  }
526 
527  // [GSL] A=A-tau w v'
528  size_t i, j, last=nc;
529 
530  for (i=ir;i<nr;i++) {
531 
532  double wi=O2SCL_IX2(M,i,ic);
533 
534  // [GSL] note, computed for v(0)=1 above
535  for (j=ic+1;j<last;j++) {
536  wi+=O2SCL_IX2(M,i,j)*O2SCL_IX2(M2,ir2,j-ic+ic2);
537  }
538 
539  // [GSL] j=0
540  O2SCL_IX2(M,i,ic)-=tau*wi;
541 
542  // [GSL] j=1 .. N-1
543  for (j=ic+1;j<last;j++) {
544  O2SCL_IX2(M,i,j)-=tau*wi*O2SCL_IX2(M2,ir2,j-ic+ic2);
545  }
546  }
547 
548  return;
549  }
550 
551 #ifdef DOXYGEN
552 }
553 #endif
void daxpy_subvec(const double alpha, const size_t N, const vec_t &X, vec2_t &Y, const size_t ie)
Compute beginning with index ie and ending with index N-1.
Definition: cblas_base.h:797
void householder_hm1(const size_t M, const size_t N, double tau, mat_t &A)
Apply a Householder transformation to a matrix being build up from the identity matrix, using the first column of A as a Householder vector.
void daxpy_subcol(const double alpha, const size_t M, const mat_t &X, const size_t ir, const size_t ic, vec_t &y)
Compute for a subcolumn of a matrix.
Definition: cblas_base.h:959
void householder_hv(const size_t N, double tau, const vec_t &v, vec2_t &w)
Apply a Householder transformation v to vector w.
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 ...
void dscal_subrow(mat_t &A, const size_t ir, const size_t ic, const size_t N, const double alpha)
Compute for a subrow of a matrix.
Definition: cblas_base.h:1305
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...
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 householder_mh(const size_t M, const size_t N, double tau, const vec_t &v, mat_t &A)
Apply the Householder transformation (v,tau) to the right-hand side of the matrix A...
double dnrm2_subcol(const mat_t &A, const size_t ir, const size_t ic, const size_t M)
Compute the norm of a subcolumn of a matrix.
Definition: cblas_base.h:1044
double dnrm2_subvec(const size_t N, const vec_t &X, const size_t ie)
Compute the norm of the vector X beginning with index ie and ending with index N-1.
Definition: cblas_base.h:875
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
double ddot_subvec(const size_t N, const vec_t &X, const vec2_t &Y, const size_t ie)
Compute beginning with index ie and ending with index N-1.
Definition: cblas_base.h:834
double ddot_subcol(const size_t M, const mat_t &X, const size_t ir, const size_t ic, const vec_t &y)
Compute for a subcolumn of a matrix.
Definition: cblas_base.h:1001
double householder_transform(const size_t n, vec_t &v)
Replace the vector v with a Householder vector and a coefficient tau that annihilates v[1] through v[...
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...
void householder_hm(const size_t M, const size_t N, double tau, const vec_t &v, mat_t &A)
Apply a Householder transformation to matrix of size (M,N)
void dscal_subvec(const double alpha, const size_t N, vec_t &X, const size_t ie)
Compute beginning with index ie and ending with index N-1.
Definition: cblas_base.h:920
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
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...
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).