diff options
Diffstat (limited to 'gsl-1.9/linalg/householdercomplex.c')
| -rw-r--r-- | gsl-1.9/linalg/householdercomplex.c | 207 |
1 files changed, 207 insertions, 0 deletions
diff --git a/gsl-1.9/linalg/householdercomplex.c b/gsl-1.9/linalg/householdercomplex.c new file mode 100644 index 0000000..bc182a6 --- /dev/null +++ b/gsl-1.9/linalg/householdercomplex.c @@ -0,0 +1,207 @@ +/* linalg/householdercomplex.c + * + * Copyright (C) 2001 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/* Computes a householder transformation matrix H such that + * + * H' v = -/+ |v| e_1 + * + * where e_1 is the first unit vector. On exit the matrix H can be + * computed from the return values (tau, v) + * + * H = I - tau * w * w' + * + * where w = (1, v(2), ..., v(N)). The nonzero element of the result + * vector -/+|v| e_1 is stored in v(1). + * + * Note that the matrix H' in the householder transformation is the + * hermitian conjugate of H. To compute H'v, pass the conjugate of + * tau as the first argument to gsl_linalg_householder_hm() rather + * than tau itself. See the LAPACK function CLARFG for details of this + * convention. */ + +#include <config.h> +#include <stdlib.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_matrix.h> +#include <gsl/gsl_blas.h> +#include <gsl/gsl_complex_math.h> + +#include <gsl/gsl_linalg.h> + +gsl_complex +gsl_linalg_complex_householder_transform (gsl_vector_complex * v) +{ + /* replace v[0:n-1] with a householder vector (v[0:n-1]) and + coefficient tau that annihilate v[1:n-1] */ + + const size_t n = v->size ; + + if (n == 1) + { + gsl_complex alpha = gsl_vector_complex_get (v, 0) ; + double absa = gsl_complex_abs (alpha); + double beta_r = - (GSL_REAL(alpha) >= 0 ? +1 : -1) * absa ; + + gsl_complex tau; + + if (beta_r == 0.0) + { + GSL_REAL(tau) = 0.0; + GSL_IMAG(tau) = 0.0; + } + else + { + GSL_REAL(tau) = (beta_r - GSL_REAL(alpha)) / beta_r ; + GSL_IMAG(tau) = - GSL_IMAG(alpha) / beta_r ; + + { + gsl_complex beta = gsl_complex_rect (beta_r, 0.0); + gsl_vector_complex_set (v, 0, beta) ; + } + } + + return tau; + } + else + { + gsl_complex tau ; + double beta_r; + + gsl_vector_complex_view x = gsl_vector_complex_subvector (v, 1, n - 1) ; + gsl_complex alpha = gsl_vector_complex_get (v, 0) ; + double absa = gsl_complex_abs (alpha); + double xnorm = gsl_blas_dznrm2 (&x.vector); + + if (xnorm == 0 && GSL_IMAG(alpha) == 0) + { + gsl_complex zero = gsl_complex_rect(0.0, 0.0); + return zero; /* tau = 0 */ + } + + beta_r = - (GSL_REAL(alpha) >= 0 ? +1 : -1) * hypot(absa, xnorm) ; + + GSL_REAL(tau) = (beta_r - GSL_REAL(alpha)) / beta_r ; + GSL_IMAG(tau) = - GSL_IMAG(alpha) / beta_r ; + + { + gsl_complex amb = gsl_complex_sub_real(alpha, beta_r); + gsl_complex s = gsl_complex_inverse(amb); + gsl_blas_zscal (s, &x.vector); + } + + { + gsl_complex beta = gsl_complex_rect (beta_r, 0.0); + gsl_vector_complex_set (v, 0, beta) ; + } + + return tau; + } +} + +int +gsl_linalg_complex_householder_hm (gsl_complex tau, const gsl_vector_complex * v, gsl_matrix_complex * A) +{ + /* applies a householder transformation v,tau to matrix m */ + + size_t i, j; + + if (GSL_REAL(tau) == 0.0 && GSL_IMAG(tau) == 0.0) + { + return GSL_SUCCESS; + } + + /* w = (v' A)^T */ + + for (j = 0; j < A->size2; j++) + { + gsl_complex tauwj; + gsl_complex wj = gsl_matrix_complex_get(A,0,j); + + for (i = 1; i < A->size1; i++) /* note, computed for v(0) = 1 above */ + { + gsl_complex Aij = gsl_matrix_complex_get(A,i,j); + gsl_complex vi = gsl_vector_complex_get(v,i); + gsl_complex Av = gsl_complex_mul (Aij, gsl_complex_conjugate(vi)); + wj = gsl_complex_add (wj, Av); + } + + tauwj = gsl_complex_mul (tau, wj); + + /* A = A - v w^T */ + + { + gsl_complex A0j = gsl_matrix_complex_get (A, 0, j); + gsl_complex Atw = gsl_complex_sub (A0j, tauwj); + /* store A0j - tau * wj */ + gsl_matrix_complex_set (A, 0, j, Atw); + } + + for (i = 1; i < A->size1; i++) + { + gsl_complex vi = gsl_vector_complex_get (v, i); + gsl_complex tauvw = gsl_complex_mul(vi, tauwj); + gsl_complex Aij = gsl_matrix_complex_get (A, i, j); + gsl_complex Atwv = gsl_complex_sub (Aij, tauvw); + /* store Aij - tau * vi * wj */ + gsl_matrix_complex_set (A, i, j, Atwv); + } + } + + return GSL_SUCCESS; +} + +int +gsl_linalg_complex_householder_hv (gsl_complex tau, const gsl_vector_complex * v, gsl_vector_complex * w) +{ + const size_t N = v->size; + + if (GSL_REAL(tau) == 0.0 && GSL_IMAG(tau) == 0.0) + return GSL_SUCCESS; + + { + /* compute z = v'w */ + + gsl_complex z0 = gsl_vector_complex_get(w,0); + gsl_complex z1, z; + gsl_complex tz, ntz; + + gsl_vector_complex_const_view v1 = gsl_vector_complex_const_subvector(v, 1, N-1); + gsl_vector_complex_view w1 = gsl_vector_complex_subvector(w, 1, N-1); + + gsl_blas_zdotc(&v1.vector, &w1.vector, &z1); + + z = gsl_complex_add (z0, z1); + + tz = gsl_complex_mul(tau, z); + ntz = gsl_complex_negative (tz); + + /* compute w = w - tau * (v'w) * v */ + + { + gsl_complex w0 = gsl_vector_complex_get(w, 0); + gsl_complex w0ntz = gsl_complex_add (w0, ntz); + gsl_vector_complex_set (w, 0, w0ntz); + } + + gsl_blas_zaxpy(ntz, &v1.vector, &w1.vector); + } + + return GSL_SUCCESS; +} |