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

   1 /* linalg/hermtd.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 /* Factorise a hermitian matrix A into
  21  *
  22  * A = U T U'
  23  *
  24  * where U is unitary and T is real symmetric tridiagonal.  Only the
  25  * diagonal and lower triangular part of A is referenced and modified.
  26  *
  27  * On exit, T is stored in the diagonal and first subdiagonal of
  28  * A. Since T is symmetric the upper diagonal is not stored.
  29  *
  30  * U is stored as a packed set of Householder transformations in the
  31  * lower triangular part of the input matrix below the first subdiagonal.
  32  *
  33  * The full matrix for Q can be obtained as the product
  34  *
  35  *       Q = Q_N ... Q_2 Q_1
  36  *
  37  * where
  38  *
  39  *       Q_i = (I - tau_i * v_i * v_i')
  40  *
  41  * and where v_i is a Householder vector
  42  *
  43  *       v_i = [0, ..., 0, 1, A(i+2,i), A(i+3,i), ... , A(N,i)]
  44  *
  45  * This storage scheme is the same as in LAPACK.  See LAPACK's
  46  * chetd2.f for details.
  47  *
  48  * See Golub & Van Loan, "Matrix Computations" (3rd ed), Section 8.3 */
  49
  50 #include <config.h>
  51 #include <stdlib.h>
  52 #include <gsl/gsl_math.h>
  53 #include <gsl/gsl_vector.h>
  54 #include <gsl/gsl_matrix.h>
  55 #include <gsl/gsl_blas.h>
  56 #include <gsl/gsl_complex_math.h>
  57
  58 #include <gsl/gsl_linalg.h>
  59
  60 int
  61 gsl_linalg_hermtd_decomp (gsl_matrix_complex * A, gsl_vector_complex * tau)
  62 {
  63   if (A->size1 != A->size2)
  64     {
  65       GSL_ERROR ("hermitian tridiagonal decomposition requires square matrix",
  66                  GSL_ENOTSQR);
  67     }
  68   else if (tau->size + 1 != A->size1)
  69     {
  70       GSL_ERROR ("size of tau must be (matrix size - 1)", GSL_EBADLEN);
  71     }
  72   else
  73     {
  74       const size_t N = A->size1;
  75       size_t i;
  76
  77       const gsl_complex zero = gsl_complex_rect (0.0, 0.0);
  78       const gsl_complex one = gsl_complex_rect (1.0, 0.0);
  79       const gsl_complex neg_one = gsl_complex_rect (-1.0, 0.0);
  80
  81       for (i = 0 ; i < N - 1; i++)
  82         {
  83           gsl_vector_complex_view c = gsl_matrix_complex_column (A, i);
  84           gsl_vector_complex_view v = gsl_vector_complex_subvector (&c.vector, i + 1, N - (i + 1));
  85           gsl_complex tau_i = gsl_linalg_complex_householder_transform (&v.vector);
  86
  87           /* Apply the transformation H^T A H to the remaining columns */
  88
  89           if ((i + 1) < (N - 1)
  90               && !(GSL_REAL(tau_i) == 0.0 && GSL_IMAG(tau_i) == 0.0))
  91             {
  92               gsl_matrix_complex_view m =
  93                 gsl_matrix_complex_submatrix (A, i + 1, i + 1,
  94                                               N - (i+1), N - (i+1));
  95               gsl_complex ei = gsl_vector_complex_get(&v.vector, 0);
  96               gsl_vector_complex_view x = gsl_vector_complex_subvector (tau, i, N-(i+1));
  97               gsl_vector_complex_set (&v.vector, 0, one);
  98
  99               /* x = tau * A * v */
 100               gsl_blas_zhemv (CblasLower, tau_i, &m.matrix, &v.vector, zero, &x.vector);
 101
 102               /* w = x - (1/2) tau * (x' * v) * v  */
 103               {
 104                 gsl_complex xv, txv, alpha;
 105                 gsl_blas_zdotc(&x.vector, &v.vector, &xv);
 106                 txv = gsl_complex_mul(tau_i, xv);
 107                 alpha = gsl_complex_mul_real(txv, -0.5);
 108                 gsl_blas_zaxpy(alpha, &v.vector, &x.vector);
 109               }
 110
 111               /* apply the transformation A = A - v w' - w v' */
 112               gsl_blas_zher2(CblasLower, neg_one, &v.vector, &x.vector, &m.matrix);
 113
 114               gsl_vector_complex_set (&v.vector, 0, ei);
 115             }
 116
 117           gsl_vector_complex_set (tau, i, tau_i);
 118         }
 119
 120       return GSL_SUCCESS;
 121     }
 122 }
 123
 124
 125 /*  Form the orthogonal matrix Q from the packed QR matrix */
 126
 127 int
 128 gsl_linalg_hermtd_unpack (const gsl_matrix_complex * A,
 129                           const gsl_vector_complex * tau,
 130                           gsl_matrix_complex * Q,
 131                           gsl_vector * diag,
 132                           gsl_vector * sdiag)
 133 {
 134   if (A->size1 !=  A->size2)
 135     {
 136       GSL_ERROR ("matrix A must be sqaure", GSL_ENOTSQR);
 137     }
 138   else if (tau->size + 1 != A->size1)
 139     {
 140       GSL_ERROR ("size of tau must be (matrix size - 1)", GSL_EBADLEN);
 141     }
 142   else if (Q->size1 != A->size1 || Q->size2 != A->size1)
 143     {
 144       GSL_ERROR ("size of Q must match size of A", GSL_EBADLEN);
 145     }
 146   else if (diag->size != A->size1)
 147     {
 148       GSL_ERROR ("size of diagonal must match size of A", GSL_EBADLEN);
 149     }
 150   else if (sdiag->size + 1 != A->size1)
 151     {
 152       GSL_ERROR ("size of subdiagonal must be (matrix size - 1)", GSL_EBADLEN);
 153     }
 154   else
 155     {
 156       const size_t N = A->size1;
 157
 158       size_t i;
 159
 160       /* Initialize Q to the identity */
 161
 162       gsl_matrix_complex_set_identity (Q);
 163
 164       for (i = N - 1; i > 0 && i--;)
 165         {
 166           gsl_complex ti = gsl_vector_complex_get (tau, i);
 167
 168           gsl_vector_complex_const_view c = gsl_matrix_complex_const_column (A, i);
 169
 170           gsl_vector_complex_const_view h =
 171             gsl_vector_complex_const_subvector (&c.vector, i + 1, N - (i+1));
 172
 173           gsl_matrix_complex_view m =
 174             gsl_matrix_complex_submatrix (Q, i + 1, i + 1, N-(i+1), N-(i+1));
 175
 176           gsl_linalg_complex_householder_hm (ti, &h.vector, &m.matrix);
 177         }
 178
 179       /* Copy diagonal into diag */
 180
 181       for (i = 0; i < N; i++)
 182         {
 183           gsl_complex Aii = gsl_matrix_complex_get (A, i, i);
 184           gsl_vector_set (diag, i, GSL_REAL(Aii));
 185         }
 186
 187       /* Copy subdiagonal into sdiag */
 188
 189       for (i = 0; i < N - 1; i++)
 190         {
 191           gsl_complex Aji = gsl_matrix_complex_get (A, i+1, i);
 192           gsl_vector_set (sdiag, i, GSL_REAL(Aji));
 193         }
 194
 195       return GSL_SUCCESS;
 196     }
 197 }
 198
 199 int
 200 gsl_linalg_hermtd_unpack_T (const gsl_matrix_complex * A,
 201                             gsl_vector * diag,
 202                             gsl_vector * sdiag)
 203 {
 204   if (A->size1 !=  A->size2)
 205     {
 206       GSL_ERROR ("matrix A must be sqaure", GSL_ENOTSQR);
 207     }
 208   else if (diag->size != A->size1)
 209     {
 210       GSL_ERROR ("size of diagonal must match size of A", GSL_EBADLEN);
 211     }
 212   else if (sdiag->size + 1 != A->size1)
 213     {
 214       GSL_ERROR ("size of subdiagonal must be (matrix size - 1)", GSL_EBADLEN);
 215     }
 216   else
 217     {
 218       const size_t N = A->size1;
 219
 220       size_t i;
 221
 222       /* Copy diagonal into diag */
 223
 224       for (i = 0; i < N; i++)
 225         {
 226           gsl_complex Aii = gsl_matrix_complex_get (A, i, i);
 227           gsl_vector_set (diag, i, GSL_REAL(Aii));
 228         }
 229
 230       /* Copy subdiagonal into sd */
 231
 232       for (i = 0; i < N - 1; i++)
 233         {
 234           gsl_complex Aji = gsl_matrix_complex_get (A, i+1, i);
 235           gsl_vector_set (sdiag, i, GSL_REAL(Aji));
 236         }
 237
 238       return GSL_SUCCESS;
 239     }
 240 }