45 #include <o2scl/err_hnd.h> 46 #include <o2scl/cblas.h> 47 #include <o2scl/permutation.h> 58 template<
class mat_t,
class vec_t>
59 int HH_svx(
size_t N,
size_t M, mat_t &A, vec_t &x) {
66 for (i = 0; i < N; i++) {
67 const double aii = O2SCL_IX2(A,i,i);
71 double max_norm = 0.0;
74 for (k = i; k < M; k++) {
75 double aki = O2SCL_IX2(A,k,i);
81 O2SCL_ERR(
"Matrix is rank deficient in HH_svx().",
85 alpha = sqrt (r) * GSL_SIGN (aii);
87 ak = 1.0 / (r + alpha * aii);
88 O2SCL_IX2(A,i,i)=aii+alpha;
92 for (k = i + 1; k < N; k++) {
95 for (j = i; j < M; j++) {
96 double ajk = O2SCL_IX2(A,j,k);
97 double aji = O2SCL_IX2(A,j,i);
101 max_norm = GSL_MAX (max_norm, norm);
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;
112 double dbl_eps=std::numeric_limits<double>::epsilon();
114 if (fabs (alpha) < 2.0*dbl_eps*sqrt(max_norm)) {
122 for (j = i; j < M; j++) {
123 f += O2SCL_IX(x,j)*O2SCL_IX2(A,j,i);
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;
135 for (i = N; i-- > 0; ) {
136 double xi = O2SCL_IX(x,i);
138 for (k = i + 1; k < N; k++) {
139 sum += O2SCL_IX2(A,i,k)*O2SCL_IX(x,k);
142 O2SCL_IX(x,i)=(xi - sum) / d[i];
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);
apparent singularity detected
int HH_svx(size_t N, size_t M, mat_t &A, vec_t &x)
Solve a linear system after Householder decomposition in place.
The namespace for linear algebra classes and functions.
int HH_solve(size_t n, mat_t &A, const vec_t &b, vec2_t &x)
Solve linear system after Householder decomposition.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.