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

   1 /* linalg/sytd.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 symmetric matrix A into
  21  *
  22  * A = Q T Q'
  23  *
  24  * where Q is orthogonal and T is 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  * Q 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_1 Q_2 ... Q_(N-2)
  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+1,i), A(i+2,i), ... , A(N,i)]
  44  *
  45  * This storage scheme is the same as in LAPACK.  See LAPACK's
  46  * ssytd2.f for details.
  47  *
  48  * See Golub & Van Loan, "Matrix Computations" (3rd ed), Section 8.3
  49  *
  50  * Note: this description uses 1-based indices. The code below uses
  51  * 0-based indices
  52  */
  53
  54 #include <config.h>
  55 #include <stdlib.h>
  56 #include <gsl/gsl_math.h>
  57 #include <gsl/gsl_vector.h>
  58 #include <gsl/gsl_matrix.h>
  59 #include <gsl/gsl_blas.h>
  60
  61 #include <gsl/gsl_linalg.h>
  62
  63 int
  64 gsl_linalg_symmtd_decomp (gsl_matrix * A, gsl_vector * tau)
  65 {
  66   if (A->size1 != A->size2)
  67     {
  68       GSL_ERROR ("symmetric tridiagonal decomposition requires square matrix",
  69                  GSL_ENOTSQR);
  70     }
  71   else if (tau->size + 1 != A->size1)
  72     {
  73       GSL_ERROR ("size of tau must be (matrix size - 1)", GSL_EBADLEN);
  74     }
  75   else
  76     {
  77       const size_t N = A->size1;
  78       size_t i;
  79
  80       for (i = 0 ; i < N - 2; i++)
  81         {
  82           gsl_vector_view c = gsl_matrix_column (A, i);
  83           gsl_vector_view v = gsl_vector_subvector (&c.vector, i + 1, N - (i + 1));
  84           double tau_i = gsl_linalg_householder_transform (&v.vector);
  85
  86           /* Apply the transformation H^T A H to the remaining columns */
  87
  88           if (tau_i != 0.0)
  89             {
  90               gsl_matrix_view m = gsl_matrix_submatrix (A, i + 1, i + 1,
  91                                                         N - (i+1), N - (i+1));
  92               double ei = gsl_vector_get(&v.vector, 0);
  93               gsl_vector_view x = gsl_vector_subvector (tau, i, N-(i+1));
  94               gsl_vector_set (&v.vector, 0, 1.0);
  95
  96               /* x = tau * A * v */
  97               gsl_blas_dsymv (CblasLower, tau_i, &m.matrix, &v.vector, 0.0, &x.vector);
  98
  99               /* w = x - (1/2) tau * (x' * v) * v  */
 100               {
 101                 double xv, alpha;
 102                 gsl_blas_ddot(&x.vector, &v.vector, &xv);
 103                 alpha = - (tau_i / 2.0) * xv;
 104                 gsl_blas_daxpy(alpha, &v.vector, &x.vector);
 105               }
 106
 107               /* apply the transformation A = A - v w' - w v' */
 108               gsl_blas_dsyr2(CblasLower, -1.0, &v.vector, &x.vector, &m.matrix);
 109
 110               gsl_vector_set (&v.vector, 0, ei);
 111             }
 112
 113           gsl_vector_set (tau, i, tau_i);
 114         }
 115
 116       return GSL_SUCCESS;
 117     }
 118 }
 119
 120
 121 /*  Form the orthogonal matrix Q from the packed QR matrix */
 122
 123 int
 124 gsl_linalg_symmtd_unpack (const gsl_matrix * A,
 125                           const gsl_vector * tau,
 126                           gsl_matrix * Q,
 127                           gsl_vector * diag,
 128                           gsl_vector * sdiag)
 129 {
 130   if (A->size1 !=  A->size2)
 131     {
 132       GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
 133     }
 134   else if (tau->size + 1 != A->size1)
 135     {
 136       GSL_ERROR ("size of tau must be (matrix size - 1)", GSL_EBADLEN);
 137     }
 138   else if (Q->size1 != A->size1 || Q->size2 != A->size1)
 139     {
 140       GSL_ERROR ("size of Q must match size of A", GSL_EBADLEN);
 141     }
 142   else if (diag->size != A->size1)
 143     {
 144       GSL_ERROR ("size of diagonal must match size of A", GSL_EBADLEN);
 145     }
 146   else if (sdiag->size + 1 != A->size1)
 147     {
 148       GSL_ERROR ("size of subdiagonal must be (matrix size - 1)", GSL_EBADLEN);
 149     }
 150   else
 151     {
 152       const size_t N = A->size1;
 153
 154       size_t i;
 155
 156       /* Initialize Q to the identity */
 157
 158       gsl_matrix_set_identity (Q);
 159
 160       for (i = N - 2; i > 0 && i--;)
 161         {
 162           gsl_vector_const_view c = gsl_matrix_const_column (A, i);
 163           gsl_vector_const_view h = gsl_vector_const_subvector (&c.vector, i + 1, N - (i+1));
 164           double ti = gsl_vector_get (tau, i);
 165
 166           gsl_matrix_view m = gsl_matrix_submatrix (Q, i + 1, i + 1, N-(i+1), N-(i+1));
 167
 168           gsl_linalg_householder_hm (ti, &h.vector, &m.matrix);
 169         }
 170
 171       /* Copy diagonal into diag */
 172
 173       for (i = 0; i < N; i++)
 174         {
 175           double Aii = gsl_matrix_get (A, i, i);
 176           gsl_vector_set (diag, i, Aii);
 177         }
 178
 179       /* Copy subdiagonal into sd */
 180
 181       for (i = 0; i < N - 1; i++)
 182         {
 183           double Aji = gsl_matrix_get (A, i+1, i);
 184           gsl_vector_set (sdiag, i, Aji);
 185         }
 186
 187       return GSL_SUCCESS;
 188     }
 189 }
 190
 191 int
 192 gsl_linalg_symmtd_unpack_T (const gsl_matrix * A,
 193                             gsl_vector * diag,
 194                             gsl_vector * sdiag)
 195 {
 196   if (A->size1 !=  A->size2)
 197     {
 198       GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
 199     }
 200   else if (diag->size != A->size1)
 201     {
 202       GSL_ERROR ("size of diagonal must match size of A", GSL_EBADLEN);
 203     }
 204   else if (sdiag->size + 1 != A->size1)
 205     {
 206       GSL_ERROR ("size of subdiagonal must be (matrix size - 1)", GSL_EBADLEN);
 207     }
 208   else
 209     {
 210       const size_t N = A->size1;
 211
 212       size_t i;
 213
 214       /* Copy diagonal into diag */
 215
 216       for (i = 0; i < N; i++)
 217         {
 218           double Aii = gsl_matrix_get (A, i, i);
 219           gsl_vector_set (diag, i, Aii);
 220         }
 221
 222       /* Copy subdiagonal into sdiag */
 223
 224       for (i = 0; i < N - 1; i++)
 225         {
 226           double Aij = gsl_matrix_get (A, i+1, i);
 227           gsl_vector_set (sdiag, i, Aij);
 228         }
 229
 230       return GSL_SUCCESS;
 231     }
 232 }