qrpt_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/qrpt.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 qrpt_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, class vec2_t>
54  void QRPT_decomp(size_t M, size_t N, mat_t &A, vec_t &tau,
55  o2scl::permutation &p, int &signum, vec2_t &norm) {
56 
57  signum=1;
58  p.init();
59 
60  for(size_t i=0;i<N;i++) {
61  norm[i]=O2SCL_CBLAS_NAMESPACE::dnrm2_subcol(A,0,i,M);
62  }
63 
64  size_t imax;
65  if (M<N) imax=M;
66  else imax=N;
67  for(size_t i=0;i<imax;i++) {
68  double max_norm=norm[i];
69  size_t kmax=i;
70  for(size_t j=i+1;j<N;j++) {
71  double x=norm[j];
72  if (x>max_norm) {
73  max_norm=x;
74  kmax=j;
75  }
76  }
77 
78  if (kmax!=i) {
80  p.swap(i,kmax);
81  o2scl::vector_swap_double(norm,i,kmax);
82  signum=-signum;
83  }
84 
85  tau[i]=householder_transform_subcol(A,i,i,M);
86  if (i+1<N) {
87  householder_hm_subcol(A,i,i+1,M,N,A,i,i,tau[i]);
88  }
89 
90  if (i+1<M) {
91  for(size_t j=i+1;j<N;j++) {
92  double x=norm[j];
93  if (x>0.0) {
94  double y=0.0;
95  double temp=O2SCL_IX2(A,i,j)/x;
96  if (fabs(temp)>=1.0) {
97  y=0.0;
98  } else {
99  y=x*sqrt(1.0-temp*temp);
100  }
101 
102  double sqrt_dbl_eps=sqrt(std::numeric_limits<double>::epsilon());
103 
104  if (fabs(y/x)<sqrt(20.0)*sqrt_dbl_eps) {
106  }
107  norm[j]=y;
108 
109  }
110  }
111  }
112 
113  }
114 
115  return;
116  }
117 
118 #ifdef DOXYGEN
119 }
120 #endif
A class for representing permutations.
Definition: permutation.h:70
void vector_swap_double(size_t N, vec_t &v1, vec2_t &v2)
Swap of of the first N elements of two double-precision vectors.
Definition: vector.h:493
void QRPT_decomp(size_t M, size_t N, mat_t &A, vec_t &tau, o2scl::permutation &p, int &signum, vec2_t &norm)
Compute the QR decomposition of matrix A.
Definition: qrpt_base.h:54
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_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
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
int init()
Initialize permutation to the identity.
Definition: permutation.h:203
void matrix_swap_cols_double(size_t M, mat_t &m, size_t j1, size_t j2)
Swap the first M rows of two columns in a double-precision matrix.
Definition: vector.h:615
int swap(const size_t i, const size_t j)
Swap two elements of a permutation.
Definition: permutation.h:249

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