49 #include <gsl/gsl_machine.h> 51 #include <o2scl/err_hnd.h> 52 #include <o2scl/cblas.h> 53 #include <o2scl/permutation.h> 77 double alpha, beta, tau;
86 beta=-(alpha >= 0.0 ? +1.0 : -1.0)*hypot(alpha,xnorm);
87 tau=(beta-alpha)/beta;
89 double dbl_eps=std::numeric_limits<double>::epsilon();
90 double dbl_min=std::numeric_limits<double>::min();
92 double s=(alpha-beta);
93 if (fabs(s)>dbl_min) {
117 template<
class mat_t>
119 const size_t ic,
const size_t M) {
125 double alpha, beta, tau;
134 alpha=O2SCL_IX2(A,ir,ic);
135 beta=-(alpha >= 0.0 ? +1.0 : -1.0)*hypot(alpha,xnorm);
136 tau=(beta-alpha)/beta;
138 double dbl_eps=std::numeric_limits<double>::epsilon();
139 double dbl_min=std::numeric_limits<double>::min();
141 double s=(alpha-beta);
142 if (fabs(s)>dbl_min) {
144 O2SCL_IX2(A,ir,ic)=beta;
148 O2SCL_IX2(A,ir,ic)=beta;
164 template<
class mat_t>
166 const size_t ic,
const size_t N) {
172 double alpha, beta, tau;
181 alpha=O2SCL_IX2(A,ir,ic);
182 beta=-(alpha >= 0.0 ? +1.0 : -1.0)*hypot(alpha,xnorm);
183 tau=(beta-alpha)/beta;
185 double dbl_eps=std::numeric_limits<double>::epsilon();
186 double dbl_min=std::numeric_limits<double>::min();
188 double s=(alpha-beta);
189 if (fabs(s)>dbl_min) {
191 O2SCL_IX2(A,ir,ic)=beta;
195 O2SCL_IX2(A,ir,ic)=beta;
209 template<
class vec_t,
class mat_t>
211 double tau,
const vec_t &v, mat_t &A) {
217 for (
size_t j=0;j<N;j++) {
220 double wj=O2SCL_IX2(A,0,j);
222 for (
size_t i=1;i<M;i++) {
223 wj+=O2SCL_IX2(A,i,j)*O2SCL_IX(v,i);
229 O2SCL_IX2(A,0,j)-=tau*wj;
232 for (
size_t i=1;i<M;i++) {
233 O2SCL_IX2(A,i,j)-=tau*wj*O2SCL_IX(v,i);
254 template<
class mat_t>
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,
265 for (
size_t j=ic;j<nc;j++) {
268 double wj=O2SCL_IX2(M,ir,j);
270 for (
size_t i=ir+1;i<nr;i++) {
271 wj+=O2SCL_IX2(M,i,j)*O2SCL_IX2(M2,i-ir+ir2,ic2);
277 O2SCL_IX2(M,ir,j)-=tau*wj;
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);
300 template<
class mat_t>
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,
311 for (
size_t j=ic;j<nc;j++) {
314 double wj=O2SCL_IX2(M,ir,j);
316 for (
size_t i=ir+1;i<nr;i++) {
317 wj+=O2SCL_IX2(M,i,j)*O2SCL_IX2(M2,ir2,i-ir+ic2);
323 O2SCL_IX2(M,ir,j)-=tau*wj;
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);
335 template<
class vec_t,
class vec2_t>
343 double d0=O2SCL_IX(w,0);
351 O2SCL_IX(w,0)-=tau*d;
363 template<
class mat_t,
class vec_t>
365 const size_t ie,
const size_t N) {
371 double d0=O2SCL_IX(w,ie);
379 O2SCL_IX(w,ie)-=tau*d;
390 template<
class mat_t>
392 double tau, mat_t &A) {
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;
399 for (
size_t i=1;i<M;i++) {
400 O2SCL_IX2(A,i,0)=0.0;
407 for (
size_t j=1;j<N;j++) {
409 for (
size_t i=1;i<M;i++) {
410 wj+=O2SCL_IX2(A,i,j)*O2SCL_IX2(A,i,0);
414 O2SCL_IX2(A,0,j)=-tau*wj;
416 for (
size_t i=1;i<M;i++) {
417 O2SCL_IX2(A,i,j)-=tau*wj*O2SCL_IX2(A,i,0);
421 for (
size_t i=1;i<M;i++) {
422 O2SCL_IX2(A,i,0)=-tau*O2SCL_IX2(A,i,0);
424 O2SCL_IX2(A,0,0)=1.0-tau;
435 template<
class mat_t>
437 double tau, mat_t &A,
size_t ir,
size_t ic) {
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;
447 for (
size_t i=irp1;i<M;i++) {
448 O2SCL_IX2(A,i,ic)=0.0;
455 for (
size_t j=icp1;j<N;j++) {
457 for (
size_t i=irp1;i<M;i++) {
458 wj+=O2SCL_IX2(A,i,j)*O2SCL_IX2(A,i,ic);
462 O2SCL_IX2(A,ir,j)=-tau*wj;
464 for (
size_t i=irp1;i<M;i++) {
465 O2SCL_IX2(A,i,j)-=tau*wj*O2SCL_IX2(A,i,ic);
469 for (
size_t i=irp1;i<M;i++) {
470 O2SCL_IX2(A,i,ic)=-tau*O2SCL_IX2(A,i,ic);
472 O2SCL_IX2(A,ir,ic)=1.0-tau;
480 template<
class vec_t,
class mat_t>
482 double tau,
const vec_t &v, mat_t &A) {
493 double wi=O2SCL_IX2(A,i,0);
497 wi+=O2SCL_IX2(A,i,j)*O2SCL_IX(v,j);
501 O2SCL_IX2(A,i,0)-=tau*wi;
505 O2SCL_IX2(A,i,j)-=tau*wi*O2SCL_IX(v,j);
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) {
528 size_t i, j, last=nc;
530 for (i=ir;i<nr;i++) {
532 double wi=O2SCL_IX2(M,i,ic);
535 for (j=ic+1;j<last;j++) {
536 wi+=O2SCL_IX2(M,i,j)*O2SCL_IX2(M2,ir2,j-ic+ic2);
540 O2SCL_IX2(M,i,ic)-=tau*wi;
543 for (j=ic+1;j<last;j++) {
544 O2SCL_IX2(M,i,j)-=tau*wi*O2SCL_IX2(M2,ir2,j-ic+ic2);
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.