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

   1 /* linalg/hesstri.c
   2  *
   3  * Copyright (C) 2006, 2007 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 <stdlib.h>
  21 #include <math.h>
  22
  23 #include <config.h>
  24 #include <gsl/gsl_linalg.h>
  25 #include <gsl/gsl_matrix.h>
  26 #include <gsl/gsl_vector.h>
  27 #include <gsl/gsl_blas.h>
  28
  29 #include "givens.c"
  30
  31 /*
  32  * This module contains routines related to the Hessenberg-Triangular
  33  * reduction of two general real matrices
  34  *
  35  * See Golub & Van Loan, "Matrix Computations", 3rd ed, sec 7.7.4
  36  */
  37
  38 /*
  39 gsl_linalg_hesstri_decomp()
  40   Perform a reduction to generalized upper Hessenberg form.
  41 Given A and B, this function overwrites A with an upper Hessenberg
  42 matrix H = U^T A V and B with an upper triangular matrix R = U^T B V
  43 with U and V orthogonal.
  44
  45 See Golub & Van Loan, "Matrix Computations" (3rd ed) algorithm 7.7.1
  46
  47 Inputs: A    - real square matrix
  48         B    - real square matrix
  49         U    - (output) if non-null, U is stored here on output
  50         V    - (output) if non-null, V is stored here on output
  51         work - workspace (length n)
  52
  53 Return: success or error
  54 */
  55
  56 int
  57 gsl_linalg_hesstri_decomp(gsl_matrix * A, gsl_matrix * B, gsl_matrix * U,
  58                           gsl_matrix * V, gsl_vector * work)
  59 {
  60   const size_t N = A->size1;
  61
  62   if ((N != A->size2) || (N != B->size1) || (N != B->size2))
  63     {
  64       GSL_ERROR ("Hessenberg-triangular reduction requires square matrices",
  65                  GSL_ENOTSQR);
  66     }
  67   else if (N != work->size)
  68     {
  69       GSL_ERROR ("length of workspace must match matrix dimension",
  70                  GSL_EBADLEN);
  71     }
  72   else
  73     {
  74       double cs, sn;          /* rotation parameters */
  75       size_t i, j;            /* looping */
  76       gsl_vector_view xv, yv; /* temporary views */
  77
  78       /* B -> Q^T B = R (upper triangular) */
  79       gsl_linalg_QR_decomp(B, work);
  80
  81       /* A -> Q^T A */
  82       gsl_linalg_QR_QTmat(B, work, A);
  83
  84       /* initialize U and V if desired */
  85
  86       if (U)
  87         {
  88           gsl_linalg_QR_unpack(B, work, U, B);
  89         }
  90       else
  91         {
  92           /* zero out lower triangle of B */
  93           for (j = 0; j < N - 1; ++j)
  94             {
  95               for (i = j + 1; i < N; ++i)
  96                 gsl_matrix_set(B, i, j, 0.0);
  97             }
  98         }
  99
 100       if (V)
 101         gsl_matrix_set_identity(V);
 102
 103       if (N < 3)
 104         return GSL_SUCCESS; /* nothing more to do */
 105
 106       /* reduce A and B */
 107       for (j = 0; j < N - 2; ++j)
 108         {
 109           for (i = N - 1; i >= (j + 2); --i)
 110             {
 111               /* step 1: rotate rows i - 1, i to kill A(i,j) */
 112
 113               /*
 114                * compute G = [ CS SN ] so that G^t [ A(i-1,j) ] = [ * ]
 115                *             [-SN CS ]             [ A(i, j)  ]   [ 0 ]
 116                */
 117               create_givens(gsl_matrix_get(A, i - 1, j),
 118                             gsl_matrix_get(A, i, j),
 119                             &cs,
 120                             &sn);
 121               /* invert so drot() works correctly (G -> G^t) */
 122               sn = -sn;
 123
 124               /* compute G^t A(i-1:i, j:n) */
 125               xv = gsl_matrix_subrow(A, i - 1, j, N - j);
 126               yv = gsl_matrix_subrow(A, i, j, N - j);
 127               gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
 128
 129               /* compute G^t B(i-1:i, i-1:n) */
 130               xv = gsl_matrix_subrow(B, i - 1, i - 1, N - i + 1);
 131               yv = gsl_matrix_subrow(B, i, i - 1, N - i + 1);
 132               gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
 133
 134               if (U)
 135                 {
 136                   /* accumulate U: U -> U G */
 137                   xv = gsl_matrix_column(U, i - 1);
 138                   yv = gsl_matrix_column(U, i);
 139                   gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
 140                 }
 141
 142               /* step 2: rotate columns i, i - 1 to kill B(i, i - 1) */
 143
 144               create_givens(-gsl_matrix_get(B, i, i),
 145                             gsl_matrix_get(B, i, i - 1),
 146                             &cs,
 147                             &sn);
 148               /* invert so drot() works correctly (G -> G^t) */
 149               sn = -sn;
 150
 151               /* compute B(1:i, i-1:i) G */
 152               xv = gsl_matrix_subcolumn(B, i - 1, 0, i + 1);
 153               yv = gsl_matrix_subcolumn(B, i, 0, i + 1);
 154               gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
 155
 156               /* apply to A(1:n, i-1:i) */
 157               xv = gsl_matrix_column(A, i - 1);
 158               yv = gsl_matrix_column(A, i);
 159               gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
 160
 161               if (V)
 162                 {
 163                   /* accumulate V: V -> V G */
 164                   xv = gsl_matrix_column(V, i - 1);
 165                   yv = gsl_matrix_column(V, i);
 166                   gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
 167                 }
 168             }
 169         }
 170
 171       return GSL_SUCCESS;
 172     }
 173 } /* gsl_linalg_hesstri_decomp() */