htswanalysis/MACS/lib/gsl/gsl-1.11/linalg/householdercomplex.c

   1 /* linalg/householdercomplex.c
   2  *
   3  * Copyright (C) 2001, 2007 Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* Computes a householder transformation matrix H such that
  21  *
  22  *       H' v = -/+ |v| e_1
  23  *
  24  * where e_1 is the first unit vector.  On exit the matrix H can be
  25  * computed from the return values (tau, v)
  26  *
  27  *       H = I - tau * w * w'
  28  *
  29  * where w = (1, v(2), ..., v(N)). The nonzero element of the result
  30  * vector -/+|v| e_1 is stored in v(1).
  31  *
  32  * Note that the matrix H' in the householder transformation is the
  33  * hermitian conjugate of H.  To compute H'v, pass the conjugate of
  34  * tau as the first argument to gsl_linalg_householder_hm() rather
  35  * than tau itself. See the LAPACK function CLARFG for details of this
  36  * convention.  */
  37
  38 #include <config.h>
  39 #include <stdlib.h>
  40 #include <gsl/gsl_math.h>
  41 #include <gsl/gsl_vector.h>
  42 #include <gsl/gsl_matrix.h>
  43 #include <gsl/gsl_blas.h>
  44 #include <gsl/gsl_complex_math.h>
  45
  46 #include <gsl/gsl_linalg.h>
  47
  48 gsl_complex
  49 gsl_linalg_complex_householder_transform (gsl_vector_complex * v)
  50 {
  51   /* replace v[0:n-1] with a householder vector (v[0:n-1]) and
  52      coefficient tau that annihilate v[1:n-1] */
  53
  54   const size_t n = v->size ;
  55
  56   if (n == 1)
  57     {
  58       gsl_complex alpha = gsl_vector_complex_get (v, 0) ;
  59       double absa = gsl_complex_abs (alpha);
  60       double beta_r = - (GSL_REAL(alpha) >= 0 ? +1 : -1) * absa ;
  61
  62       gsl_complex tau;
  63
  64       if (beta_r == 0.0)
  65         {
  66           GSL_REAL(tau) = 0.0;
  67           GSL_IMAG(tau) = 0.0;
  68         }
  69       else
  70         {
  71           GSL_REAL(tau) = (beta_r - GSL_REAL(alpha)) / beta_r ;
  72           GSL_IMAG(tau) = - GSL_IMAG(alpha) / beta_r ;
  73
  74           {
  75             gsl_complex beta = gsl_complex_rect (beta_r, 0.0);
  76             gsl_vector_complex_set (v, 0, beta) ;
  77           }
  78         }
  79
  80       return tau;
  81     }
  82   else
  83     {
  84       gsl_complex tau ;
  85       double beta_r;
  86
  87       gsl_vector_complex_view x = gsl_vector_complex_subvector (v, 1, n - 1) ;
  88       gsl_complex alpha = gsl_vector_complex_get (v, 0) ;
  89       double absa = gsl_complex_abs (alpha);
  90       double xnorm = gsl_blas_dznrm2 (&x.vector);
  91
  92       if (xnorm == 0 && GSL_IMAG(alpha) == 0)
  93         {
  94           gsl_complex zero = gsl_complex_rect(0.0, 0.0);
  95           return zero; /* tau = 0 */
  96         }
  97
  98       beta_r = - (GSL_REAL(alpha) >= 0 ? +1 : -1) * hypot(absa, xnorm) ;
  99
 100       GSL_REAL(tau) = (beta_r - GSL_REAL(alpha)) / beta_r ;
 101       GSL_IMAG(tau) = - GSL_IMAG(alpha) / beta_r ;
 102
 103       {
 104         gsl_complex amb = gsl_complex_sub_real(alpha, beta_r);
 105         gsl_complex s = gsl_complex_inverse(amb);
 106         gsl_blas_zscal (s, &x.vector);
 107       }
 108
 109       {
 110         gsl_complex beta = gsl_complex_rect (beta_r, 0.0);
 111         gsl_vector_complex_set (v, 0, beta) ;
 112       }
 113
 114       return tau;
 115     }
 116 }
 117
 118 int
 119 gsl_linalg_complex_householder_hm (gsl_complex tau, const gsl_vector_complex * v, gsl_matrix_complex * A)
 120 {
 121   /* applies a householder transformation v,tau to matrix m */
 122
 123   size_t i, j;
 124
 125   if (GSL_REAL(tau) == 0.0 && GSL_IMAG(tau) == 0.0)
 126     {
 127       return GSL_SUCCESS;
 128     }
 129
 130   /* w = (v' A)^T */
 131
 132   for (j = 0; j < A->size2; j++)
 133     {
 134       gsl_complex tauwj;
 135       gsl_complex wj = gsl_matrix_complex_get(A,0,j);
 136
 137       for (i = 1; i < A->size1; i++)  /* note, computed for v(0) = 1 above */
 138         {
 139           gsl_complex Aij = gsl_matrix_complex_get(A,i,j);
 140           gsl_complex vi = gsl_vector_complex_get(v,i);
 141           gsl_complex Av = gsl_complex_mul (Aij, gsl_complex_conjugate(vi));
 142           wj = gsl_complex_add (wj, Av);
 143         }
 144
 145       tauwj = gsl_complex_mul (tau, wj);
 146
 147       /* A = A - v w^T */
 148
 149       {
 150         gsl_complex A0j = gsl_matrix_complex_get (A, 0, j);
 151         gsl_complex Atw = gsl_complex_sub (A0j, tauwj);
 152         /* store A0j - tau  * wj */
 153         gsl_matrix_complex_set (A, 0, j, Atw);
 154       }
 155
 156       for (i = 1; i < A->size1; i++)
 157         {
 158           gsl_complex vi = gsl_vector_complex_get (v, i);
 159           gsl_complex tauvw = gsl_complex_mul(vi, tauwj);
 160           gsl_complex Aij = gsl_matrix_complex_get (A, i, j);
 161           gsl_complex Atwv = gsl_complex_sub (Aij, tauvw);
 162           /* store Aij - tau * vi * wj */
 163           gsl_matrix_complex_set (A, i, j, Atwv);
 164         }
 165     }
 166
 167   return GSL_SUCCESS;
 168 }
 169
 170 int
 171 gsl_linalg_complex_householder_mh (gsl_complex tau, const gsl_vector_complex * v, gsl_matrix_complex * A)
 172 {
 173   /* applies a householder transformation v,tau to matrix m on the right */
 174
 175   size_t i, j;
 176
 177   if (GSL_REAL(tau) == 0.0 && GSL_IMAG(tau) == 0.0)
 178     {
 179       return GSL_SUCCESS;
 180     }
 181
 182   /* A -> A - A*tau*v*v^h */
 183
 184   for (i = 0; i < A->size1; i++)
 185     {
 186       gsl_complex tauwi;
 187       gsl_complex Ai0 = gsl_matrix_complex_get (A, i, 0);
 188       gsl_complex wi = Ai0;
 189
 190       /* compute w = A v */
 191       for (j = 1; j < A->size2; j++)  /* note, computed for v(0) = 1 above */
 192         {
 193           gsl_complex Aij = gsl_matrix_complex_get(A, i, j);
 194           gsl_complex vj = gsl_vector_complex_get(v, j);
 195           gsl_complex Av = gsl_complex_mul (Aij, vj);
 196           wi = gsl_complex_add (wi, Av);
 197         }
 198
 199       tauwi = gsl_complex_mul (tau, wi);
 200
 201       /* A = A - w v^H */
 202
 203       {
 204         gsl_complex Atw = gsl_complex_sub (Ai0, tauwi);
 205         /* store Ai0 - tau  * wi */
 206         gsl_matrix_complex_set (A, i, 0, Atw);
 207       }
 208
 209       for (j = 1; j < A->size2; j++)
 210         {
 211           gsl_complex vj = gsl_vector_complex_get (v, j);
 212           gsl_complex tauwv = gsl_complex_mul(gsl_complex_conjugate(vj), tauwi);
 213           gsl_complex Aij = gsl_matrix_complex_get (A, i, j);
 214           gsl_complex Atwv = gsl_complex_sub (Aij, tauwv);
 215           /* store Aij - tau * wi * conj(vj) */
 216           gsl_matrix_complex_set (A, i, j, Atwv);
 217         }
 218     }
 219
 220   return GSL_SUCCESS;
 221 }
 222
 223 int
 224 gsl_linalg_complex_householder_hv (gsl_complex tau, const gsl_vector_complex * v, gsl_vector_complex *  w)
 225 {
 226   const size_t N = v->size;
 227
 228   if (GSL_REAL(tau) == 0.0 && GSL_IMAG(tau) == 0.0)
 229       return GSL_SUCCESS;
 230
 231   {
 232     /* compute z = v'w */
 233
 234     gsl_complex z0 = gsl_vector_complex_get(w,0);
 235     gsl_complex z1, z;
 236     gsl_complex tz, ntz;
 237
 238     gsl_vector_complex_const_view v1 = gsl_vector_complex_const_subvector(v, 1, N-1);
 239     gsl_vector_complex_view w1 = gsl_vector_complex_subvector(w, 1, N-1);
 240
 241     gsl_blas_zdotc(&v1.vector, &w1.vector, &z1);
 242
 243     z = gsl_complex_add (z0, z1);
 244
 245     tz = gsl_complex_mul(tau, z);
 246     ntz = gsl_complex_negative (tz);
 247
 248     /* compute w = w - tau * (v'w) * v   */
 249
 250     {
 251       gsl_complex w0 = gsl_vector_complex_get(w, 0);
 252       gsl_complex w0ntz = gsl_complex_add (w0, ntz);
 253       gsl_vector_complex_set (w, 0, w0ntz);
 254     }
 255
 256     gsl_blas_zaxpy(ntz, &v1.vector, &w1.vector);
 257   }
 258
 259   return GSL_SUCCESS;
 260 }