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

   1 /* linalg/hessenberg.c
   2  *
   3  * Copyright (C) 2006 Patrick Alken
   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 <gsl/gsl_linalg.h>
  22 #include <gsl/gsl_matrix.h>
  23 #include <gsl/gsl_vector.h>
  24
  25 /*
  26 gsl_linalg_hessenberg_decomp()
  27   Compute the Householder reduction to Hessenberg form of a
  28 square N-by-N matrix A.
  29
  30 H = U^t A U
  31
  32 See Golub & Van Loan, "Matrix Computations" (3rd ed), algorithm
  33 7.4.2
  34
  35 Inputs: A   - matrix to reduce
  36         tau - where to store scalar factors in Householder
  37               matrices; this vector must be of length N,
  38               where N is the order of A
  39
  40 Return: GSL_SUCCESS unless error occurs
  41
  42 Notes: on output, the upper triangular portion of A (including
  43 the diagaonal and subdiagonal) contains the Hessenberg matrix.
  44 The lower triangular portion (below the subdiagonal) contains
  45 the Householder vectors which can be used to construct
  46 the similarity transform matrix U.
  47
  48 The matrix U is
  49
  50 U = U(1) U(2) ... U(n - 2)
  51
  52 where
  53
  54 U(i) = I - tau(i) * v(i) * v(i)^t
  55
  56 and the vector v(i) is stored in column i of the matrix A
  57 underneath the subdiagonal. So the first element of v(i)
  58 is stored in row i + 2, column i, the second element at
  59 row i + 3, column i, and so on.
  60
  61 Also note that for the purposes of computing U(i),
  62 v(1:i) = 0, v(i + 1) = 1, and v(i+2:n) is what is stored in
  63 column i of A beneath the subdiagonal.
  64 */
  65
  66 int
  67 gsl_linalg_hessenberg_decomp(gsl_matrix *A, gsl_vector *tau)
  68 {
  69   const size_t N = A->size1;
  70
  71   if (N != A->size2)
  72     {
  73       GSL_ERROR ("Hessenberg reduction requires square matrix",
  74                  GSL_ENOTSQR);
  75     }
  76   else if (N != tau->size)
  77     {
  78       GSL_ERROR ("tau vector must match matrix size", GSL_EBADLEN);
  79     }
  80   else if (N < 3)
  81     {
  82       /* nothing to do */
  83       return GSL_SUCCESS;
  84     }
  85   else
  86     {
  87       size_t i;           /* looping */
  88       gsl_vector_view c,  /* matrix column */
  89                       hv; /* householder vector */
  90       gsl_matrix_view m;
  91       double tau_i;       /* beta in algorithm 7.4.2 */
  92
  93       for (i = 0; i < N - 2; ++i)
  94         {
  95           /*
  96            * make a copy of A(i + 1:n, i) and store it in the section
  97            * of 'tau' that we haven't stored coefficients in yet
  98            */
  99
 100           c = gsl_matrix_subcolumn(A, i, i + 1, N - i - 1);
 101
 102           hv = gsl_vector_subvector(tau, i + 1, N - (i + 1));
 103           gsl_vector_memcpy(&hv.vector, &c.vector);
 104
 105           /* compute householder transformation of A(i+1:n,i) */
 106           tau_i = gsl_linalg_householder_transform(&hv.vector);
 107
 108           /* apply left householder matrix (I - tau_i v v') to A */
 109           m = gsl_matrix_submatrix(A, i + 1, i, N - (i + 1), N - i);
 110           gsl_linalg_householder_hm(tau_i, &hv.vector, &m.matrix);
 111
 112           /* apply right householder matrix (I - tau_i v v') to A */
 113           m = gsl_matrix_submatrix(A, 0, i + 1, N, N - (i + 1));
 114           gsl_linalg_householder_mh(tau_i, &hv.vector, &m.matrix);
 115
 116           /* save Householder coefficient */
 117           gsl_vector_set(tau, i, tau_i);
 118
 119           /*
 120            * store Householder vector below the subdiagonal in column
 121            * i of the matrix. hv(1) does not need to be stored since
 122            * it is always 1.
 123            */
 124           c = gsl_vector_subvector(&c.vector, 1, c.vector.size - 1);
 125           hv = gsl_vector_subvector(&hv.vector, 1, hv.vector.size - 1);
 126           gsl_vector_memcpy(&c.vector, &hv.vector);
 127         }
 128
 129       return GSL_SUCCESS;
 130     }
 131 } /* gsl_linalg_hessenberg_decomp() */
 132
 133 /*
 134 gsl_linalg_hessenberg_unpack()
 135   Construct the matrix U which transforms a matrix A into
 136 its upper Hessenberg form:
 137
 138 H = U^t A U
 139
 140 by unpacking the information stored in H from gsl_linalg_hessenberg().
 141
 142 U is a product of Householder matrices:
 143
 144 U = U(1) U(2) ... U(n - 2)
 145
 146 where
 147
 148 U(i) = I - tau(i) * v(i) * v(i)^t
 149
 150 The v(i) are stored in the lower triangular part of H by
 151 gsl_linalg_hessenberg(). The tau(i) are stored in the vector tau.
 152
 153 Inputs: H       - Hessenberg matrix computed from
 154                   gsl_linalg_hessenberg()
 155         tau     - tau vector computed from gsl_linalg_hessenberg()
 156         U       - (output) where to store similarity matrix
 157
 158 Return: success or error
 159 */
 160
 161 int
 162 gsl_linalg_hessenberg_unpack(gsl_matrix * H, gsl_vector * tau,
 163                              gsl_matrix * U)
 164 {
 165   int s;
 166
 167   gsl_matrix_set_identity(U);
 168
 169   s = gsl_linalg_hessenberg_unpack_accum(H, tau, U);
 170
 171   return s;
 172 } /* gsl_linalg_hessenberg_unpack() */
 173
 174 /*
 175 gsl_linalg_hessenberg_unpack_accum()
 176   This routine is the same as gsl_linalg_hessenberg_unpack(), except
 177 instead of storing the similarity matrix in U, it accumulates it,
 178 so that
 179
 180 U -> U * [ U(1) U(2) ... U(n - 2) ]
 181
 182 instead of:
 183
 184 U -> U(1) U(2) ... U(n - 2)
 185
 186 Inputs: H       - Hessenberg matrix computed from
 187                   gsl_linalg_hessenberg()
 188         tau     - tau vector computed from gsl_linalg_hessenberg()
 189         V       - (input/output) where to accumulate similarity matrix
 190
 191 Return: success or error
 192
 193 Notes: 1) On input, V needs to be initialized. The Householder matrices
 194           are accumulated into V, so on output,
 195
 196             V_out = V_in * U(1) * U(2) * ... * U(n - 2)
 197
 198           so if you just want the product of the Householder matrices,
 199           initialize V to the identity matrix before calling this
 200           function.
 201
 202        2) V does not have to be square, but must have the same
 203           number of columns as the order of H
 204 */
 205
 206 int
 207 gsl_linalg_hessenberg_unpack_accum(gsl_matrix * H, gsl_vector * tau,
 208                                    gsl_matrix * V)
 209 {
 210   const size_t N = H->size1;
 211
 212   if (N != H->size2)
 213     {
 214       GSL_ERROR ("Hessenberg reduction requires square matrix",
 215                  GSL_ENOTSQR);
 216     }
 217   else if (N != tau->size)
 218     {
 219       GSL_ERROR ("tau vector must match matrix size", GSL_EBADLEN);
 220     }
 221   else if (N != V->size2)
 222     {
 223       GSL_ERROR ("V matrix has wrong dimension", GSL_EBADLEN);
 224     }
 225   else
 226     {
 227       size_t j;           /* looping */
 228       double tau_j;       /* householder coefficient */
 229       gsl_vector_view c,  /* matrix column */
 230                       hv; /* householder vector */
 231       gsl_matrix_view m;
 232
 233       if (N < 3)
 234         {
 235           /* nothing to do */
 236           return GSL_SUCCESS;
 237         }
 238
 239       for (j = 0; j < (N - 2); ++j)
 240         {
 241           c = gsl_matrix_column(H, j);
 242
 243           tau_j = gsl_vector_get(tau, j);
 244
 245           /*
 246            * get a view to the householder vector in column j, but
 247            * make sure hv(2) starts at the element below the
 248            * subdiagonal, since hv(1) was never stored and is always
 249            * 1
 250            */
 251           hv = gsl_vector_subvector(&c.vector, j + 1, N - (j + 1));
 252
 253           /*
 254            * Only operate on part of the matrix since the first
 255            * j + 1 entries of the real householder vector are 0
 256            *
 257            * V -> V * U(j)
 258            *
 259            * Note here that V->size1 is not necessarily equal to N
 260            */
 261           m = gsl_matrix_submatrix(V, 0, j + 1, V->size1, N - (j + 1));
 262
 263           /* apply right Householder matrix to V */
 264           gsl_linalg_householder_mh(tau_j, &hv.vector, &m.matrix);
 265         }
 266
 267       return GSL_SUCCESS;
 268     }
 269 } /* gsl_linalg_hessenberg_unpack_accum() */
 270
 271 /*
 272 gsl_linalg_hessenberg_set_zero()
 273   Zero out the lower triangular portion of the Hessenberg matrix H.
 274 This is useful when Householder vectors may be stored in the lower
 275 part of H, but eigenvalue solvers need some scratch space with zeros.
 276 */
 277
 278 int
 279 gsl_linalg_hessenberg_set_zero(gsl_matrix * H)
 280 {
 281   const size_t N = H->size1;
 282   size_t i, j;
 283
 284   if (N < 3)
 285     return GSL_SUCCESS;
 286
 287   for (j = 0; j < N - 2; ++j)
 288     {
 289       for (i = j + 2; i < N; ++i)
 290         {
 291           gsl_matrix_set(H, i, j, 0.0);
 292         }
 293     }
 294
 295   return GSL_SUCCESS;
 296 } /* gsl_linalg_hessenberg_set_zero() */
 297
 298 /*
 299 gsl_linalg_hessenberg_submatrix()
 300
 301   This routine does the same thing as gsl_linalg_hessenberg(),
 302 except that it operates on a submatrix of a larger matrix, but
 303 updates the larger matrix with the Householder transformations.
 304
 305 For example, suppose
 306
 307 M = [ M_{11} | M_{12} | M_{13} ]
 308     [   0    |   A    | M_{23} ]
 309     [   0    |   0    | M_{33} ]
 310
 311 where M_{11} and M_{33} are already in Hessenberg form, and we
 312 just want to reduce A to Hessenberg form. Applying the transformations
 313 to A alone will cause the larger matrix M to lose its similarity
 314 information. So this routine updates M_{12} and M_{23} as A gets
 315 reduced.
 316
 317 Inputs: M   - total matrix
 318         A   - (sub)matrix to reduce
 319         top - row index of top of A in M
 320         tau - where to store scalar factors in Householder
 321               matrices; this vector must be of length N,
 322               where N is the order of A
 323
 324 Return: GSL_SUCCESS unless error occurs
 325
 326 Notes: on output, the upper triangular portion of A (including
 327 the diagaonal and subdiagonal) contains the Hessenberg matrix.
 328 The lower triangular portion (below the subdiagonal) contains
 329 the Householder vectors which can be used to construct
 330 the similarity transform matrix U.
 331
 332 The matrix U is
 333
 334 U = U(1) U(2) ... U(n - 2)
 335
 336 where
 337
 338 U(i) = I - tau(i) * v(i) * v(i)^t
 339
 340 and the vector v(i) is stored in column i of the matrix A
 341 underneath the subdiagonal. So the first element of v(i)
 342 is stored in row i + 2, column i, the second element at
 343 row i + 3, column i, and so on.
 344
 345 Also note that for the purposes of computing U(i),
 346 v(1:i) = 0, v(i + 1) = 1, and v(i+2:n) is what is stored in
 347 column i of A beneath the subdiagonal.
 348 */
 349
 350 int
 351 gsl_linalg_hessenberg_submatrix(gsl_matrix *M, gsl_matrix *A, size_t top,
 352                                 gsl_vector *tau)
 353 {
 354   const size_t N = A->size1;
 355   const size_t N_M = M->size1;
 356
 357   if (N != A->size2)
 358     {
 359       GSL_ERROR ("Hessenberg reduction requires square matrix",
 360                  GSL_ENOTSQR);
 361     }
 362   else if (N != tau->size)
 363     {
 364       GSL_ERROR ("tau vector must match matrix size", GSL_EBADLEN);
 365     }
 366   else if (N < 3)
 367     {
 368       /* nothing to do */
 369       return GSL_SUCCESS;
 370     }
 371   else
 372     {
 373       size_t i;           /* looping */
 374       gsl_vector_view c,  /* matrix column */
 375                       hv; /* householder vector */
 376       gsl_matrix_view m;
 377       double tau_i;       /* beta in algorithm 7.4.2 */
 378
 379       for (i = 0; i < N - 2; ++i)
 380         {
 381           /*
 382            * make a copy of A(i + 1:n, i) and store it in the section
 383            * of 'tau' that we haven't stored coefficients in yet
 384            */
 385
 386           c = gsl_matrix_subcolumn(A, i, i + 1, N - i - 1);
 387
 388           hv = gsl_vector_subvector(tau, i + 1, N - (i + 1));
 389           gsl_vector_memcpy(&hv.vector, &c.vector);
 390
 391           /* compute householder transformation of A(i+1:n,i) */
 392           tau_i = gsl_linalg_householder_transform(&hv.vector);
 393
 394           /*
 395            * apply left householder matrix (I - tau_i v v') to
 396            * [ A | M_{23} ]
 397            */
 398           m = gsl_matrix_submatrix(M,
 399                                    top + i + 1,
 400                                    top + i,
 401                                    N - (i + 1),
 402                                    N_M - top - i);
 403           gsl_linalg_householder_hm(tau_i, &hv.vector, &m.matrix);
 404
 405           /*
 406            * apply right householder matrix (I - tau_i v v') to
 407            *
 408            * [ M_{12} ]
 409            * [   A    ]
 410            */
 411           m = gsl_matrix_submatrix(M,
 412                                    0,
 413                                    top + i + 1,
 414                                    top + N,
 415                                    N - (i + 1));
 416           gsl_linalg_householder_mh(tau_i, &hv.vector, &m.matrix);
 417
 418           /* save Householder coefficient */
 419           gsl_vector_set(tau, i, tau_i);
 420
 421           /*
 422            * store Householder vector below the subdiagonal in column
 423            * i of the matrix. hv(1) does not need to be stored since
 424            * it is always 1.
 425            */
 426           c = gsl_vector_subvector(&c.vector, 1, c.vector.size - 1);
 427           hv = gsl_vector_subvector(&hv.vector, 1, hv.vector.size - 1);
 428           gsl_vector_memcpy(&c.vector, &hv.vector);
 429         }
 430
 431       return GSL_SUCCESS;
 432     }
 433 } /* gsl_linalg_hessenberg_submatrix() */
 434
 435 /* To support gsl-1.9 interface: DEPRECATED */
 436 int
 437 gsl_linalg_hessenberg(gsl_matrix *A, gsl_vector *tau)
 438 {
 439   return gsl_linalg_hessenberg_decomp(A, tau);
 440 }