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

   1 /* linalg/householder.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004, 2007 Gerard Jungman, 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 #include <config.h>
  21 #include <stdlib.h>
  22 #include <gsl/gsl_math.h>
  23 #include <gsl/gsl_vector.h>
  24 #include <gsl/gsl_matrix.h>
  25 #include <gsl/gsl_blas.h>
  26
  27 #include <gsl/gsl_linalg.h>
  28
  29 double
  30 gsl_linalg_householder_transform (gsl_vector * v)
  31 {
  32   /* replace v[0:n-1] with a householder vector (v[0:n-1]) and
  33      coefficient tau that annihilate v[1:n-1] */
  34
  35   const size_t n = v->size ;
  36
  37   if (n == 1)
  38     {
  39       return 0.0; /* tau = 0 */
  40     }
  41   else
  42     {
  43       double alpha, beta, tau ;
  44
  45       gsl_vector_view x = gsl_vector_subvector (v, 1, n - 1) ;
  46
  47       double xnorm = gsl_blas_dnrm2 (&x.vector);
  48
  49       if (xnorm == 0)
  50         {
  51           return 0.0; /* tau = 0 */
  52         }
  53
  54       alpha = gsl_vector_get (v, 0) ;
  55       beta = - (alpha >= 0.0 ? +1.0 : -1.0) * hypot(alpha, xnorm) ;
  56       tau = (beta - alpha) / beta ;
  57
  58       gsl_blas_dscal (1.0 / (alpha - beta), &x.vector);
  59       gsl_vector_set (v, 0, beta) ;
  60
  61       return tau;
  62     }
  63 }
  64
  65 int
  66 gsl_linalg_householder_hm (double tau, const gsl_vector * v, gsl_matrix * A)
  67 {
  68   /* applies a householder transformation v,tau to matrix m */
  69
  70   if (tau == 0.0)
  71     {
  72       return GSL_SUCCESS;
  73     }
  74
  75 #ifdef USE_BLAS
  76   {
  77     gsl_vector_const_view v1 = gsl_vector_const_subvector (v, 1, v->size - 1);
  78     gsl_matrix_view A1 = gsl_matrix_submatrix (A, 1, 0, A->size1 - 1, A->size2);
  79     size_t j;
  80
  81     for (j = 0; j < A->size2; j++)
  82       {
  83         double wj = 0.0;
  84         gsl_vector_view A1j = gsl_matrix_column(&A1.matrix, j);
  85         gsl_blas_ddot (&A1j.vector, &v1.vector, &wj);
  86         wj += gsl_matrix_get(A,0,j);
  87
  88         {
  89           double A0j = gsl_matrix_get (A, 0, j);
  90           gsl_matrix_set (A, 0, j, A0j - tau *  wj);
  91         }
  92
  93         gsl_blas_daxpy (-tau * wj, &v1.vector, &A1j.vector);
  94       }
  95   }
  96 #else
  97   {
  98     size_t i, j;
  99
 100     for (j = 0; j < A->size2; j++)
 101       {
 102         /* Compute wj = Akj vk */
 103
 104         double wj = gsl_matrix_get(A,0,j);
 105
 106         for (i = 1; i < A->size1; i++)  /* note, computed for v(0) = 1 above */
 107           {
 108             wj += gsl_matrix_get(A,i,j) * gsl_vector_get(v,i);
 109           }
 110
 111         /* Aij = Aij - tau vi wj */
 112
 113         /* i = 0 */
 114         {
 115           double A0j = gsl_matrix_get (A, 0, j);
 116           gsl_matrix_set (A, 0, j, A0j - tau *  wj);
 117         }
 118
 119         /* i = 1 .. M-1 */
 120
 121         for (i = 1; i < A->size1; i++)
 122           {
 123             double Aij = gsl_matrix_get (A, i, j);
 124             double vi = gsl_vector_get (v, i);
 125             gsl_matrix_set (A, i, j, Aij - tau * vi * wj);
 126           }
 127       }
 128   }
 129 #endif
 130
 131   return GSL_SUCCESS;
 132 }
 133
 134 int
 135 gsl_linalg_householder_mh (double tau, const gsl_vector * v, gsl_matrix * A)
 136 {
 137   /* applies a householder transformation v,tau to matrix m from the
 138      right hand side in order to zero out rows */
 139
 140   if (tau == 0)
 141     return GSL_SUCCESS;
 142
 143   /* A = A - tau w v' */
 144
 145 #ifdef USE_BLAS
 146   {
 147     gsl_vector_const_view v1 = gsl_vector_const_subvector (v, 1, v->size - 1);
 148     gsl_matrix_view A1 = gsl_matrix_submatrix (A, 0, 1, A->size1, A->size2-1);
 149     size_t i;
 150
 151     for (i = 0; i < A->size1; i++)
 152       {
 153         double wi = 0.0;
 154         gsl_vector_view A1i = gsl_matrix_row(&A1.matrix, i);
 155         gsl_blas_ddot (&A1i.vector, &v1.vector, &wi);
 156         wi += gsl_matrix_get(A,i,0);
 157
 158         {
 159           double Ai0 = gsl_matrix_get (A, i, 0);
 160           gsl_matrix_set (A, i, 0, Ai0 - tau *  wi);
 161         }
 162
 163         gsl_blas_daxpy(-tau * wi, &v1.vector, &A1i.vector);
 164       }
 165   }
 166 #else
 167   {
 168     size_t i, j;
 169
 170     for (i = 0; i < A->size1; i++)
 171       {
 172         double wi = gsl_matrix_get(A,i,0);
 173
 174         for (j = 1; j < A->size2; j++)  /* note, computed for v(0) = 1 above */
 175           {
 176             wi += gsl_matrix_get(A,i,j) * gsl_vector_get(v,j);
 177           }
 178
 179         /* j = 0 */
 180
 181         {
 182           double Ai0 = gsl_matrix_get (A, i, 0);
 183           gsl_matrix_set (A, i, 0, Ai0 - tau *  wi);
 184         }
 185
 186         /* j = 1 .. N-1 */
 187
 188         for (j = 1; j < A->size2; j++)
 189           {
 190             double vj = gsl_vector_get (v, j);
 191             double Aij = gsl_matrix_get (A, i, j);
 192             gsl_matrix_set (A, i, j, Aij - tau * wi * vj);
 193           }
 194       }
 195   }
 196 #endif
 197
 198   return GSL_SUCCESS;
 199 }
 200
 201 int
 202 gsl_linalg_householder_hv (double tau, const gsl_vector * v, gsl_vector * w)
 203 {
 204   /* applies a householder transformation v to vector w */
 205   const size_t N = v->size;
 206
 207   if (tau == 0)
 208     return GSL_SUCCESS ;
 209
 210   {
 211     /* compute d = v'w */
 212
 213     double d0 = gsl_vector_get(w,0);
 214     double d1, d;
 215
 216     gsl_vector_const_view v1 = gsl_vector_const_subvector(v, 1, N-1);
 217     gsl_vector_view w1 = gsl_vector_subvector(w, 1, N-1);
 218
 219     gsl_blas_ddot (&v1.vector, &w1.vector, &d1);
 220
 221     d = d0 + d1;
 222
 223     /* compute w = w - tau (v) (v'w) */
 224
 225     {
 226       double w0 = gsl_vector_get (w,0);
 227       gsl_vector_set (w, 0, w0 - tau * d);
 228     }
 229
 230     gsl_blas_daxpy (-tau * d, &v1.vector, &w1.vector);
 231   }
 232
 233   return GSL_SUCCESS;
 234 }
 235
 236
 237 int
 238 gsl_linalg_householder_hm1 (double tau, gsl_matrix * A)
 239 {
 240   /* applies a householder transformation v,tau to a matrix being
 241      build up from the identity matrix, using the first column of A as
 242      a householder vector */
 243
 244   if (tau == 0)
 245     {
 246       size_t i,j;
 247
 248       gsl_matrix_set (A, 0, 0, 1.0);
 249
 250       for (j = 1; j < A->size2; j++)
 251         {
 252           gsl_matrix_set (A, 0, j, 0.0);
 253         }
 254
 255       for (i = 1; i < A->size1; i++)
 256         {
 257           gsl_matrix_set (A, i, 0, 0.0);
 258         }
 259
 260       return GSL_SUCCESS;
 261     }
 262
 263   /* w = A' v */
 264
 265 #ifdef USE_BLAS
 266   {
 267     gsl_matrix_view A1 = gsl_matrix_submatrix (A, 1, 0, A->size1 - 1, A->size2);
 268     gsl_vector_view v1 = gsl_matrix_column (&A1.matrix, 0);
 269     size_t j;
 270
 271     for (j = 1; j < A->size2; j++)
 272       {
 273         double wj = 0.0;   /* A0j * v0 */
 274
 275         gsl_vector_view A1j = gsl_matrix_column(&A1.matrix, j);
 276         gsl_blas_ddot (&A1j.vector, &v1.vector, &wj);
 277
 278         /* A = A - tau v w' */
 279
 280         gsl_matrix_set (A, 0, j, - tau *  wj);
 281
 282         gsl_blas_daxpy(-tau*wj, &v1.vector, &A1j.vector);
 283       }
 284
 285     gsl_blas_dscal(-tau, &v1.vector);
 286
 287     gsl_matrix_set (A, 0, 0, 1.0 - tau);
 288   }
 289 #else
 290   {
 291     size_t i, j;
 292
 293     for (j = 1; j < A->size2; j++)
 294       {
 295         double wj = 0.0;   /* A0j * v0 */
 296
 297         for (i = 1; i < A->size1; i++)
 298           {
 299             double vi = gsl_matrix_get(A, i, 0);
 300             wj += gsl_matrix_get(A,i,j) * vi;
 301           }
 302
 303         /* A = A - tau v w' */
 304
 305         gsl_matrix_set (A, 0, j, - tau *  wj);
 306
 307         for (i = 1; i < A->size1; i++)
 308           {
 309             double vi = gsl_matrix_get (A, i, 0);
 310             double Aij = gsl_matrix_get (A, i, j);
 311             gsl_matrix_set (A, i, j, Aij - tau * vi * wj);
 312           }
 313       }
 314
 315     for (i = 1; i < A->size1; i++)
 316       {
 317         double vi = gsl_matrix_get(A, i, 0);
 318         gsl_matrix_set(A, i, 0, -tau * vi);
 319       }
 320
 321     gsl_matrix_set (A, 0, 0, 1.0 - tau);
 322   }
 323 #endif
 324
 325   return GSL_SUCCESS;
 326 }