hh_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/hh.c
24  *
25  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, 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 hh_base.h
43  \brief File for householder solver
44 */
45 #include <o2scl/err_hnd.h>
46 #include <o2scl/cblas.h>
47 #include <o2scl/permutation.h>
48 
49 #ifdef DOXYGEN
50 namespace o2scl_linalg {
51 #endif
52 
53  /** \brief Solve a linear system after Householder decomposition in place
54 
55  \future Handle memory allocation like the tridiagonal
56  functions.
57  */
58  template<class mat_t, class vec_t>
59  int HH_svx(size_t N, size_t M, mat_t &A, vec_t &x) {
60  size_t i, j, k;
61 
63 
64  /* Perform Householder transformation. */
65 
66  for (i = 0; i < N; i++) {
67  const double aii = O2SCL_IX2(A,i,i);
68  double alpha;
69  double f;
70  double ak;
71  double max_norm = 0.0;
72  double r = 0.0;
73 
74  for (k = i; k < M; k++) {
75  double aki = O2SCL_IX2(A,k,i);
76  r += aki * aki;
77  }
78 
79  if (r == 0.0) {
80  /* Rank of matrix is less than size1. */
81  O2SCL_ERR("Matrix is rank deficient in HH_svx().",
83  }
84 
85  alpha = sqrt (r) * GSL_SIGN (aii);
86 
87  ak = 1.0 / (r + alpha * aii);
88  O2SCL_IX2(A,i,i)=aii+alpha;
89 
90  d[i] = -alpha;
91 
92  for (k = i + 1; k < N; k++) {
93  double norm = 0.0;
94  f = 0.0;
95  for (j = i; j < M; j++) {
96  double ajk = O2SCL_IX2(A,j,k);
97  double aji = O2SCL_IX2(A,j,i);
98  norm += ajk * ajk;
99  f += ajk * aji;
100  }
101  max_norm = GSL_MAX (max_norm, norm);
102 
103  f *= ak;
104 
105  for (j = i; j < M; j++) {
106  double ajk = O2SCL_IX2(A,j,k);
107  double aji = O2SCL_IX2(A,j,i);
108  O2SCL_IX2(A,j,k)=ajk-f*aji;
109  }
110  }
111 
112  double dbl_eps=std::numeric_limits<double>::epsilon();
113 
114  if (fabs (alpha) < 2.0*dbl_eps*sqrt(max_norm)) {
115  /* Apparent singularity. */
116  O2SCL_ERR("Apparent singularity in HH_svx().",o2scl::exc_esing);
117  }
118 
119  /* Perform update of RHS. */
120 
121  f = 0.0;
122  for (j = i; j < M; j++) {
123  f += O2SCL_IX(x,j)*O2SCL_IX2(A,j,i);
124  }
125  f *= ak;
126  for (j = i; j < M; j++) {
127  double xj = O2SCL_IX(x,j);
128  double aji = O2SCL_IX2(A,j,i);
129  O2SCL_IX(x,j)=xj - f * aji;
130  }
131  }
132 
133  /* Perform back-substitution. */
134 
135  for (i = N; i-- > 0; ) {
136  double xi = O2SCL_IX(x,i);
137  double sum = 0.0;
138  for (k = i + 1; k < N; k++) {
139  sum += O2SCL_IX2(A,i,k)*O2SCL_IX(x,k);
140  }
141 
142  O2SCL_IX(x,i)=(xi - sum) / d[i];
143  }
144 
145  return o2scl::success;
146  }
147 
148  /** \brief Solve linear system after Householder decomposition
149  */
150  template<class mat_t, class vec_t, class vec2_t>
151  int HH_solve(size_t n, mat_t &A, const vec_t &b, vec2_t &x) {
152  for(size_t i=0;i<n;i++) O2SCL_IX(x,i)=O2SCL_IX(b,i);
153  return HH_svx(n,n,A,x);
154  }
155 
156 #ifdef DOXYGEN
157 }
158 #endif
apparent singularity detected
Definition: err_hnd.h:93
int HH_svx(size_t N, size_t M, mat_t &A, vec_t &x)
Solve a linear system after Householder decomposition in place.
Definition: hh_base.h:59
The namespace for linear algebra classes and functions.
Definition: bidiag.h:36
int HH_solve(size_t n, mat_t &A, const vec_t &b, vec2_t &x)
Solve linear system after Householder decomposition.
Definition: hh_base.h:151
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
Success.
Definition: err_hnd.h:47

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