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

   1 /* linalg/ptlq.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
   4  * Copyright (C) 2004 Joerg Wensch, modifications for LQ.
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU General Public License as published by
   8  * the Free Software Foundation; either version 3 of the License, or (at
   9  * your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful, but
  12  * WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License
  17  * along with this program; if not, write to the Free Software
  18  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  19  */
  20
  21 #include <config.h>
  22 #include <stdlib.h>
  23 #include <string.h>
  24 #include <gsl/gsl_blas.h>
  25 #include <gsl/gsl_math.h>
  26 #include <gsl/gsl_vector.h>
  27 #include <gsl/gsl_matrix.h>
  28 #include <gsl/gsl_permute_vector.h>
  29 #include <gsl/gsl_linalg.h>
  30
  31 #include "givens.c"
  32 #include "apply_givens.c"
  33
  34 /* The purpose of this package is to speed up QR-decomposition for
  35    large matrices.  Because QR-decomposition is column oriented, but
  36    GSL uses a row-oriented matrix format, there can considerable
  37    speedup obtained by computing the LQ-decomposition of the
  38    transposed matrix instead.  This package provides LQ-decomposition
  39    and related algorithms.  */
  40
  41 /* Factorise a general N x M matrix A into
  42  *
  43  *   P A = L Q
  44  *
  45  * where Q is orthogonal (M x M) and L is lower triangular (N x M).
  46  * When A is rank deficient, r = rank(A) < n, then the permutation is
  47  * used to ensure that the lower n - r columns of L are zero and the first
  48  * l rows of Q form an orthonormal basis for the rows of A.
  49  *
  50  * Q is stored as a packed set of Householder transformations in the
  51  * strict upper triangular part of the input matrix.
  52  *
  53  * L is stored in the diagonal and lower triangle of the input matrix.
  54  *
  55  * P: column j of P is column k of the identity matrix, where k =
  56  * permutation->data[j]
  57  *
  58  * The full matrix for Q can be obtained as the product
  59  *
  60  *       Q = Q_k .. Q_2 Q_1
  61  *
  62  * where k = MIN(M,N) and
  63  *
  64  *       Q_i = (I - tau_i * v_i * v_i')
  65  *
  66  * and where v_i is a Householder vector
  67  *
  68  *       v_i = [1, m(i,i+1), m(i,i+2), ... , m(i,M)]
  69  *
  70  * This storage scheme is the same as in LAPACK.  See LAPACK's
  71  * dgeqpf.f for details.
  72  *
  73  */
  74
  75 int
  76 gsl_linalg_PTLQ_decomp (gsl_matrix * A, gsl_vector * tau, gsl_permutation * p, int *signum, gsl_vector * norm)
  77 {
  78   const size_t N = A->size1;
  79   const size_t M = A->size2;
  80
  81   if (tau->size != GSL_MIN (M, N))
  82     {
  83       GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
  84     }
  85   else if (p->size != N)
  86     {
  87       GSL_ERROR ("permutation size must be N", GSL_EBADLEN);
  88     }
  89   else if (norm->size != N)
  90     {
  91       GSL_ERROR ("norm size must be N", GSL_EBADLEN);
  92     }
  93   else
  94     {
  95       size_t i;
  96
  97       *signum = 1;
  98
  99       gsl_permutation_init (p); /* set to identity */
 100
 101       /* Compute column norms and store in workspace */
 102
 103       for (i = 0; i < N; i++)
 104         {
 105           gsl_vector_view c = gsl_matrix_row (A, i);
 106           double x = gsl_blas_dnrm2 (&c.vector);
 107           gsl_vector_set (norm, i, x);
 108         }
 109
 110       for (i = 0; i < GSL_MIN (M, N); i++)
 111         {
 112           /* Bring the column of largest norm into the pivot position */
 113
 114           double max_norm = gsl_vector_get(norm, i);
 115           size_t j, kmax = i;
 116
 117           for (j = i + 1; j < N; j++)
 118             {
 119               double x = gsl_vector_get (norm, j);
 120
 121               if (x > max_norm)
 122                 {
 123                   max_norm = x;
 124                   kmax = j;
 125                 }
 126             }
 127
 128           if (kmax != i)
 129             {
 130               gsl_matrix_swap_rows (A, i, kmax);
 131               gsl_permutation_swap (p, i, kmax);
 132               gsl_vector_swap_elements(norm,i,kmax);
 133
 134               (*signum) = -(*signum);
 135             }
 136
 137           /* Compute the Householder transformation to reduce the j-th
 138              column of the matrix to a multiple of the j-th unit vector */
 139
 140           {
 141             gsl_vector_view c_full = gsl_matrix_row (A, i);
 142             gsl_vector_view c = gsl_vector_subvector (&c_full.vector,
 143                                                       i, M - i);
 144             double tau_i = gsl_linalg_householder_transform (&c.vector);
 145
 146             gsl_vector_set (tau, i, tau_i);
 147
 148             /* Apply the transformation to the remaining columns */
 149
 150             if (i + 1 < N)
 151               {
 152                 gsl_matrix_view m = gsl_matrix_submatrix (A, i +1, i, N - (i+1), M - i);
 153
 154                 gsl_linalg_householder_mh (tau_i, &c.vector, &m.matrix);
 155               }
 156           }
 157
 158           /* Update the norms of the remaining columns too */
 159
 160           if (i + 1 < M)
 161             {
 162               for (j = i + 1; j < N; j++)
 163                 {
 164                   double x = gsl_vector_get (norm, j);
 165
 166                   if (x > 0.0)
 167                     {
 168                       double y = 0;
 169                       double temp= gsl_matrix_get (A, j, i) / x;
 170
 171                       if (fabs (temp) >= 1)
 172                         y = 0.0;
 173                       else
 174                         y = x * sqrt (1 - temp * temp);
 175
 176                       /* recompute norm to prevent loss of accuracy */
 177
 178                       if (fabs (y / x) < sqrt (20.0) * GSL_SQRT_DBL_EPSILON)
 179                         {
 180                           gsl_vector_view c_full = gsl_matrix_row (A, j);
 181                           gsl_vector_view c =
 182                             gsl_vector_subvector(&c_full.vector,
 183                                                  i+1, M - (i+1));
 184                           y = gsl_blas_dnrm2 (&c.vector);
 185                         }
 186
 187                       gsl_vector_set (norm, j, y);
 188                     }
 189                 }
 190             }
 191         }
 192
 193       return GSL_SUCCESS;
 194     }
 195 }
 196
 197 int
 198 gsl_linalg_PTLQ_decomp2 (const gsl_matrix * A, gsl_matrix * q, gsl_matrix * r, gsl_vector * tau, gsl_permutation * p, int *signum, gsl_vector * norm)
 199 {
 200   const size_t N = A->size1;
 201   const size_t M = A->size2;
 202
 203   if (q->size1 != M || q->size2 !=M)
 204     {
 205       GSL_ERROR ("q must be M x M", GSL_EBADLEN);
 206     }
 207   else if (r->size1 != N || r->size2 !=M)
 208     {
 209       GSL_ERROR ("r must be N x M", GSL_EBADLEN);
 210     }
 211   else if (tau->size != GSL_MIN (M, N))
 212     {
 213       GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
 214     }
 215   else if (p->size != N)
 216     {
 217       GSL_ERROR ("permutation size must be N", GSL_EBADLEN);
 218     }
 219   else if (norm->size != N)
 220     {
 221       GSL_ERROR ("norm size must be N", GSL_EBADLEN);
 222     }
 223
 224   gsl_matrix_memcpy (r, A);
 225
 226   gsl_linalg_PTLQ_decomp (r, tau, p, signum, norm);
 227
 228   /* FIXME:  aliased arguments depends on behavior of unpack routine! */
 229
 230   gsl_linalg_LQ_unpack (r, tau, q, r);
 231
 232   return GSL_SUCCESS;
 233 }
 234
 235
 236 /* Solves the system x^T A = b^T using the P^T L Q  factorisation,
 237
 238    z^T L = b^T Q^T
 239
 240    x = P z;
 241
 242    to obtain x. Based on SLATEC code. */
 243
 244 int
 245 gsl_linalg_PTLQ_solve_T (const gsl_matrix * QR,
 246                          const gsl_vector * tau,
 247                          const gsl_permutation * p,
 248                          const gsl_vector * b,
 249                          gsl_vector * x)
 250 {
 251   if (QR->size1 != QR->size2)
 252     {
 253       GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
 254     }
 255   else if (QR->size2 != p->size)
 256     {
 257       GSL_ERROR ("matrix size must match permutation size", GSL_EBADLEN);
 258     }
 259   else if (QR->size2 != b->size)
 260     {
 261       GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
 262     }
 263   else if (QR->size1 != x->size)
 264     {
 265       GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
 266     }
 267   else
 268     {
 269       gsl_vector_memcpy (x, b);
 270
 271       gsl_linalg_PTLQ_svx_T (QR, tau, p, x);
 272
 273       return GSL_SUCCESS;
 274     }
 275 }
 276
 277 int
 278 gsl_linalg_PTLQ_svx_T (const gsl_matrix * LQ,
 279                        const gsl_vector * tau,
 280                        const gsl_permutation * p,
 281                        gsl_vector * x)
 282 {
 283   if (LQ->size1 != LQ->size2)
 284     {
 285       GSL_ERROR ("LQ matrix must be square", GSL_ENOTSQR);
 286     }
 287   else if (LQ->size2 != p->size)
 288     {
 289       GSL_ERROR ("matrix size must match permutation size", GSL_EBADLEN);
 290     }
 291   else if (LQ->size1 != x->size)
 292     {
 293       GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
 294     }
 295   else
 296     {
 297       /* compute sol = b^T Q^T */
 298
 299       gsl_linalg_LQ_vecQT (LQ, tau, x);
 300
 301       /* Solve  L^T x = sol, storing x inplace in sol */
 302
 303       gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, LQ, x);
 304
 305       gsl_permute_vector_inverse (p, x);
 306
 307       return GSL_SUCCESS;
 308     }
 309 }
 310
 311
 312 int
 313 gsl_linalg_PTLQ_LQsolve_T (const gsl_matrix * Q, const gsl_matrix * L,
 314                            const gsl_permutation * p,
 315                            const gsl_vector * b,
 316                            gsl_vector * x)
 317 {
 318   if (Q->size1 != Q->size2 || L->size1 != L->size2)
 319     {
 320       return GSL_ENOTSQR;
 321     }
 322   else if (Q->size1 != p->size || Q->size1 != L->size1
 323            || Q->size1 != b->size)
 324     {
 325       return GSL_EBADLEN;
 326     }
 327   else
 328     {
 329       /* compute b' = Q b */
 330
 331       gsl_blas_dgemv (CblasNoTrans, 1.0, Q, b, 0.0, x);
 332
 333       /* Solve L^T x = b', storing x inplace */
 334
 335       gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, L, x);
 336
 337       /* Apply permutation to solution in place */
 338
 339       gsl_permute_vector_inverse (p, x);
 340
 341       return GSL_SUCCESS;
 342     }
 343 }
 344
 345 int
 346 gsl_linalg_PTLQ_Lsolve_T (const gsl_matrix * LQ,
 347                         const gsl_permutation * p,
 348                         const gsl_vector * b,
 349                         gsl_vector * x)
 350 {
 351   if (LQ->size1 != LQ->size2)
 352     {
 353       GSL_ERROR ("LQ matrix must be square", GSL_ENOTSQR);
 354     }
 355   else if (LQ->size1 != b->size)
 356     {
 357       GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
 358     }
 359   else if (LQ->size2 != x->size)
 360     {
 361       GSL_ERROR ("matrix size must match x size", GSL_EBADLEN);
 362     }
 363   else if (p->size != x->size)
 364     {
 365       GSL_ERROR ("permutation size must match x size", GSL_EBADLEN);
 366     }
 367   else
 368     {
 369       /* Copy x <- b */
 370
 371       gsl_vector_memcpy (x, b);
 372
 373       /* Solve L^T x = b, storing x inplace */
 374
 375       gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, LQ, x);
 376
 377       gsl_permute_vector_inverse (p, x);
 378
 379       return GSL_SUCCESS;
 380     }
 381 }
 382
 383
 384 int
 385 gsl_linalg_PTLQ_Lsvx_T (const gsl_matrix * LQ,
 386                         const gsl_permutation * p,
 387                         gsl_vector * x)
 388 {
 389   if (LQ->size1 != LQ->size2)
 390     {
 391       GSL_ERROR ("LQ matrix must be square", GSL_ENOTSQR);
 392     }
 393   else if (LQ->size2 != x->size)
 394     {
 395       GSL_ERROR ("matrix size must match x size", GSL_EBADLEN);
 396     }
 397   else if (p->size != x->size)
 398     {
 399       GSL_ERROR ("permutation size must match x size", GSL_EBADLEN);
 400     }
 401   else
 402     {
 403       /* Solve L^T x = b, storing x inplace */
 404
 405       gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, LQ, x);
 406
 407       gsl_permute_vector_inverse (p, x);
 408
 409       return GSL_SUCCESS;
 410     }
 411 }
 412
 413
 414
 415 /* Update a P^T L Q factorisation for P A= L Q ,  A' =  A +  v u^T,
 416                                                  PA' = PA + Pv u^T
 417
 418  * P^T L' Q' = P^T LQ + v u^T
 419  *       = P^T (L + (P v) u^T Q^T) Q
 420  *       = P^T (L + (P v) w^T) Q
 421  *
 422  * where w = Q^T u.
 423  *
 424  * Algorithm from Golub and Van Loan, "Matrix Computations", Section
 425  * 12.5 (Updating Matrix Factorizations, Rank-One Changes)
 426  */
 427
 428 int
 429 gsl_linalg_PTLQ_update (gsl_matrix * Q, gsl_matrix * L,
 430                         const gsl_permutation * p,
 431                         const gsl_vector * v, gsl_vector * w)
 432 {
 433   if (Q->size1 != Q->size2 || L->size1 != L->size2)
 434     {
 435       return GSL_ENOTSQR;
 436     }
 437   else if (L->size1 != Q->size2 || v->size != Q->size2 || w->size != Q->size2)
 438     {
 439       return GSL_EBADLEN;
 440     }
 441   else
 442     {
 443       size_t j, k;
 444       const size_t N = Q->size1;
 445       const size_t M = Q->size2;
 446       double w0;
 447
 448       /* Apply Given's rotations to reduce w to (|w|, 0, 0, ... , 0)
 449
 450          J_1^T .... J_(n-1)^T w = +/- |w| e_1
 451
 452          simultaneously applied to L,  H = J_1^T ... J^T_(n-1) L
 453          so that H is upper Hessenberg.  (12.5.2) */
 454
 455       for (k = M - 1; k > 0; k--)
 456         {
 457           double c, s;
 458           double wk = gsl_vector_get (w, k);
 459           double wkm1 = gsl_vector_get (w, k - 1);
 460
 461           create_givens (wkm1, wk, &c, &s);
 462           apply_givens_vec (w, k - 1, k, c, s);
 463           apply_givens_lq (M, N, Q, L, k - 1, k, c, s);
 464         }
 465
 466       w0 = gsl_vector_get (w, 0);
 467
 468       /* Add in v w^T  (Equation 12.5.3) */
 469
 470       for (j = 0; j < N; j++)
 471         {
 472           double lj0 = gsl_matrix_get (L, j, 0);
 473           size_t p_j = gsl_permutation_get (p, j);
 474           double vj = gsl_vector_get (v, p_j);
 475           gsl_matrix_set (L, j, 0, lj0 + w0 * vj);
 476         }
 477
 478       /* Apply Givens transformations L' = G_(n-1)^T ... G_1^T H
 479          Equation 12.5.4 */
 480
 481       for (k = 1; k < N; k++)
 482         {
 483           double c, s;
 484           double diag = gsl_matrix_get (L, k - 1, k - 1);
 485           double offdiag = gsl_matrix_get (L, k - 1, k );
 486
 487           create_givens (diag, offdiag, &c, &s);
 488           apply_givens_lq (M, N, Q, L, k - 1, k, c, s);
 489         }
 490
 491       return GSL_SUCCESS;
 492     }
 493 }