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

   1 /* linalg/bidiag.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 matrix A into
  21  *
  22  * A = U B V^T
  23  *
  24  * where U and V are orthogonal and B is upper bidiagonal.
  25  *
  26  * On exit, B is stored in the diagonal and first superdiagonal of A.
  27  *
  28  * U is stored as a packed set of Householder transformations in the
  29  * lower triangular part of the input matrix below the diagonal.
  30  *
  31  * V is stored as a packed set of Householder transformations in the
  32  * upper triangular part of the input matrix above the first
  33  * superdiagonal.
  34  *
  35  * The full matrix for U can be obtained as the product
  36  *
  37  *       U = U_1 U_2 .. U_N
  38  *
  39  * where
  40  *
  41  *       U_i = (I - tau_i * u_i * u_i')
  42  *
  43  * and where u_i is a Householder vector
  44  *
  45  *       u_i = [0, .. , 0, 1, A(i+1,i), A(i+3,i), .. , A(M,i)]
  46  *
  47  * The full matrix for V can be obtained as the product
  48  *
  49  *       V = V_1 V_2 .. V_(N-2)
  50  *
  51  * where
  52  *
  53  *       V_i = (I - tau_i * v_i * v_i')
  54  *
  55  * and where v_i is a Householder vector
  56  *
  57  *       v_i = [0, .. , 0, 1, A(i,i+2), A(i,i+3), .. , A(i,N)]
  58  *
  59  * See Golub & Van Loan, "Matrix Computations" (3rd ed), Algorithm 5.4.2
  60  *
  61  * Note: this description uses 1-based indices. The code below uses
  62  * 0-based indices
  63  */
  64
  65 #include <config.h>
  66 #include <stdlib.h>
  67 #include <gsl/gsl_math.h>
  68 #include <gsl/gsl_vector.h>
  69 #include <gsl/gsl_matrix.h>
  70 #include <gsl/gsl_blas.h>
  71
  72 #include <gsl/gsl_linalg.h>
  73
  74 int
  75 gsl_linalg_bidiag_decomp (gsl_matrix * A, gsl_vector * tau_U, gsl_vector * tau_V)
  76 {
  77   if (A->size1 < A->size2)
  78     {
  79       GSL_ERROR ("bidiagonal decomposition requires M>=N", GSL_EBADLEN);
  80     }
  81   else if (tau_U->size  != A->size2)
  82     {
  83       GSL_ERROR ("size of tau_U must be N", GSL_EBADLEN);
  84     }
  85   else if (tau_V->size + 1 != A->size2)
  86     {
  87       GSL_ERROR ("size of tau_V must be (N - 1)", GSL_EBADLEN);
  88     }
  89   else
  90     {
  91       const size_t M = A->size1;
  92       const size_t N = A->size2;
  93       size_t i;
  94
  95       for (i = 0 ; i < N; i++)
  96         {
  97           /* Apply Householder transformation to current column */
  98
  99           {
 100             gsl_vector_view c = gsl_matrix_column (A, i);
 101             gsl_vector_view v = gsl_vector_subvector (&c.vector, i, M - i);
 102             double tau_i = gsl_linalg_householder_transform (&v.vector);
 103
 104             /* Apply the transformation to the remaining columns */
 105
 106             if (i + 1 < N)
 107               {
 108                 gsl_matrix_view m =
 109                   gsl_matrix_submatrix (A, i, i + 1, M - i, N - (i + 1));
 110                 gsl_linalg_householder_hm (tau_i, &v.vector, &m.matrix);
 111               }
 112
 113             gsl_vector_set (tau_U, i, tau_i);
 114
 115           }
 116
 117           /* Apply Householder transformation to current row */
 118
 119           if (i + 1 < N)
 120             {
 121               gsl_vector_view r = gsl_matrix_row (A, i);
 122               gsl_vector_view v = gsl_vector_subvector (&r.vector, i + 1, N - (i + 1));
 123               double tau_i = gsl_linalg_householder_transform (&v.vector);
 124
 125               /* Apply the transformation to the remaining rows */
 126
 127               if (i + 1 < M)
 128                 {
 129                   gsl_matrix_view m =
 130                     gsl_matrix_submatrix (A, i+1, i+1, M - (i+1), N - (i+1));
 131                   gsl_linalg_householder_mh (tau_i, &v.vector, &m.matrix);
 132                 }
 133
 134               gsl_vector_set (tau_V, i, tau_i);
 135             }
 136         }
 137     }
 138
 139   return GSL_SUCCESS;
 140 }
 141
 142 /* Form the orthogonal matrices U, V, diagonal d and superdiagonal sd
 143    from the packed bidiagonal matrix A */
 144
 145 int
 146 gsl_linalg_bidiag_unpack (const gsl_matrix * A,
 147                           const gsl_vector * tau_U,
 148                           gsl_matrix * U,
 149                           const gsl_vector * tau_V,
 150                           gsl_matrix * V,
 151                           gsl_vector * diag,
 152                           gsl_vector * superdiag)
 153 {
 154   const size_t M = A->size1;
 155   const size_t N = A->size2;
 156
 157   const size_t K = GSL_MIN(M, N);
 158
 159   if (M < N)
 160     {
 161       GSL_ERROR ("matrix A must have M >= N", GSL_EBADLEN);
 162     }
 163   else if (tau_U->size != K)
 164     {
 165       GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
 166     }
 167   else if (tau_V->size + 1 != K)
 168     {
 169       GSL_ERROR ("size of tau must be MIN(M,N) - 1", GSL_EBADLEN);
 170     }
 171   else if (U->size1 != M || U->size2 != N)
 172     {
 173       GSL_ERROR ("size of U must be M x N", GSL_EBADLEN);
 174     }
 175   else if (V->size1 != N || V->size2 != N)
 176     {
 177       GSL_ERROR ("size of V must be N x N", GSL_EBADLEN);
 178     }
 179   else if (diag->size != K)
 180     {
 181       GSL_ERROR ("size of diagonal must match size of A", GSL_EBADLEN);
 182     }
 183   else if (superdiag->size + 1 != K)
 184     {
 185       GSL_ERROR ("size of subdiagonal must be (diagonal size - 1)", GSL_EBADLEN);
 186     }
 187   else
 188     {
 189       size_t i, j;
 190
 191       /* Copy diagonal into diag */
 192
 193       for (i = 0; i < N; i++)
 194         {
 195           double Aii = gsl_matrix_get (A, i, i);
 196           gsl_vector_set (diag, i, Aii);
 197         }
 198
 199       /* Copy superdiagonal into superdiag */
 200
 201       for (i = 0; i < N - 1; i++)
 202         {
 203           double Aij = gsl_matrix_get (A, i, i+1);
 204           gsl_vector_set (superdiag, i, Aij);
 205         }
 206
 207       /* Initialize V to the identity */
 208
 209       gsl_matrix_set_identity (V);
 210
 211       for (i = N - 1; i > 0 && i--;)
 212         {
 213           /* Householder row transformation to accumulate V */
 214           gsl_vector_const_view r = gsl_matrix_const_row (A, i);
 215           gsl_vector_const_view h =
 216             gsl_vector_const_subvector (&r.vector, i + 1, N - (i+1));
 217
 218           double ti = gsl_vector_get (tau_V, i);
 219
 220           gsl_matrix_view m =
 221             gsl_matrix_submatrix (V, i + 1, i + 1, N-(i+1), N-(i+1));
 222
 223           gsl_linalg_householder_hm (ti, &h.vector, &m.matrix);
 224         }
 225
 226       /* Initialize U to the identity */
 227
 228       gsl_matrix_set_identity (U);
 229
 230       for (j = N; j > 0 && j--;)
 231         {
 232           /* Householder column transformation to accumulate U */
 233           gsl_vector_const_view c = gsl_matrix_const_column (A, j);
 234           gsl_vector_const_view h = gsl_vector_const_subvector (&c.vector, j, M - j);
 235           double tj = gsl_vector_get (tau_U, j);
 236
 237           gsl_matrix_view m =
 238             gsl_matrix_submatrix (U, j, j, M-j, N-j);
 239
 240           gsl_linalg_householder_hm (tj, &h.vector, &m.matrix);
 241         }
 242
 243       return GSL_SUCCESS;
 244     }
 245 }
 246
 247 int
 248 gsl_linalg_bidiag_unpack2 (gsl_matrix * A,
 249                            gsl_vector * tau_U,
 250                            gsl_vector * tau_V,
 251                            gsl_matrix * V)
 252 {
 253   const size_t M = A->size1;
 254   const size_t N = A->size2;
 255
 256   const size_t K = GSL_MIN(M, N);
 257
 258   if (M < N)
 259     {
 260       GSL_ERROR ("matrix A must have M >= N", GSL_EBADLEN);
 261     }
 262   else if (tau_U->size != K)
 263     {
 264       GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
 265     }
 266   else if (tau_V->size + 1 != K)
 267     {
 268       GSL_ERROR ("size of tau must be MIN(M,N) - 1", GSL_EBADLEN);
 269     }
 270   else if (V->size1 != N || V->size2 != N)
 271     {
 272       GSL_ERROR ("size of V must be N x N", GSL_EBADLEN);
 273     }
 274   else
 275     {
 276       size_t i, j;
 277
 278       /* Initialize V to the identity */
 279
 280       gsl_matrix_set_identity (V);
 281
 282       for (i = N - 1; i > 0 && i--;)
 283         {
 284           /* Householder row transformation to accumulate V */
 285           gsl_vector_const_view r = gsl_matrix_const_row (A, i);
 286           gsl_vector_const_view h =
 287             gsl_vector_const_subvector (&r.vector, i + 1, N - (i+1));
 288
 289           double ti = gsl_vector_get (tau_V, i);
 290
 291           gsl_matrix_view m =
 292             gsl_matrix_submatrix (V, i + 1, i + 1, N-(i+1), N-(i+1));
 293
 294           gsl_linalg_householder_hm (ti, &h.vector, &m.matrix);
 295         }
 296
 297       /* Copy superdiagonal into tau_v */
 298
 299       for (i = 0; i < N - 1; i++)
 300         {
 301           double Aij = gsl_matrix_get (A, i, i+1);
 302           gsl_vector_set (tau_V, i, Aij);
 303         }
 304
 305       /* Allow U to be unpacked into the same memory as A, copy
 306          diagonal into tau_U */
 307
 308       for (j = N; j > 0 && j--;)
 309         {
 310           /* Householder column transformation to accumulate U */
 311           double tj = gsl_vector_get (tau_U, j);
 312           double Ajj = gsl_matrix_get (A, j, j);
 313           gsl_matrix_view m = gsl_matrix_submatrix (A, j, j, M-j, N-j);
 314
 315           gsl_vector_set (tau_U, j, Ajj);
 316           gsl_linalg_householder_hm1 (tj, &m.matrix);
 317         }
 318
 319       return GSL_SUCCESS;
 320     }
 321 }
 322
 323
 324 int
 325 gsl_linalg_bidiag_unpack_B (const gsl_matrix * A,
 326                             gsl_vector * diag,
 327                             gsl_vector * superdiag)
 328 {
 329   const size_t M = A->size1;
 330   const size_t N = A->size2;
 331
 332   const size_t K = GSL_MIN(M, N);
 333
 334   if (diag->size != K)
 335     {
 336       GSL_ERROR ("size of diagonal must match size of A", GSL_EBADLEN);
 337     }
 338   else if (superdiag->size + 1 != K)
 339     {
 340       GSL_ERROR ("size of subdiagonal must be (matrix size - 1)", GSL_EBADLEN);
 341     }
 342   else
 343     {
 344       size_t i;
 345
 346       /* Copy diagonal into diag */
 347
 348       for (i = 0; i < K; i++)
 349         {
 350           double Aii = gsl_matrix_get (A, i, i);
 351           gsl_vector_set (diag, i, Aii);
 352         }
 353
 354       /* Copy superdiagonal into superdiag */
 355
 356       for (i = 0; i < K - 1; i++)
 357         {
 358           double Aij = gsl_matrix_get (A, i, i+1);
 359           gsl_vector_set (superdiag, i, Aij);
 360         }
 361
 362       return GSL_SUCCESS;
 363     }
 364 }