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

   1 /* linalg/qr.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, 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 /* Author:  G. Jungman */
  21
  22 #include <config.h>
  23 #include <stdlib.h>
  24 #include <string.h>
  25 #include <gsl/gsl_math.h>
  26 #include <gsl/gsl_vector.h>
  27 #include <gsl/gsl_matrix.h>
  28 #include <gsl/gsl_blas.h>
  29
  30 #include <gsl/gsl_linalg.h>
  31
  32 #define REAL double
  33
  34 #include "givens.c"
  35 #include "apply_givens.c"
  36
  37 /* Factorise a general M x N matrix A into
  38  *
  39  *   A = Q R
  40  *
  41  * where Q is orthogonal (M x M) and R is upper triangular (M x N).
  42  *
  43  * Q is stored as a packed set of Householder transformations in the
  44  * strict lower triangular part of the input matrix.
  45  *
  46  * R is stored in the diagonal and upper triangle of the input matrix.
  47  *
  48  * The full matrix for Q can be obtained as the product
  49  *
  50  *       Q = Q_k .. Q_2 Q_1
  51  *
  52  * where k = MIN(M,N) and
  53  *
  54  *       Q_i = (I - tau_i * v_i * v_i')
  55  *
  56  * and where v_i is a Householder vector
  57  *
  58  *       v_i = [1, m(i+1,i), m(i+2,i), ... , m(M,i)]
  59  *
  60  * This storage scheme is the same as in LAPACK.  */
  61
  62 int
  63 gsl_linalg_QR_decomp (gsl_matrix * A, gsl_vector * tau)
  64 {
  65   const size_t M = A->size1;
  66   const size_t N = A->size2;
  67
  68   if (tau->size != GSL_MIN (M, N))
  69     {
  70       GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
  71     }
  72   else
  73     {
  74       size_t i;
  75
  76       for (i = 0; i < GSL_MIN (M, N); i++)
  77         {
  78           /* Compute the Householder transformation to reduce the j-th
  79              column of the matrix to a multiple of the j-th unit vector */
  80
  81           gsl_vector_view c_full = gsl_matrix_column (A, i);
  82           gsl_vector_view c = gsl_vector_subvector (&(c_full.vector), i, M-i);
  83
  84           double tau_i = gsl_linalg_householder_transform (&(c.vector));
  85
  86           gsl_vector_set (tau, i, tau_i);
  87
  88           /* Apply the transformation to the remaining columns and
  89              update the norms */
  90
  91           if (i + 1 < N)
  92             {
  93               gsl_matrix_view m = gsl_matrix_submatrix (A, i, i + 1, M - i, N - (i + 1));
  94               gsl_linalg_householder_hm (tau_i, &(c.vector), &(m.matrix));
  95             }
  96         }
  97
  98       return GSL_SUCCESS;
  99     }
 100 }
 101
 102 /* Solves the system A x = b using the QR factorisation,
 103
 104  *  R x = Q^T b
 105  *
 106  * to obtain x. Based on SLATEC code.
 107  */
 108
 109 int
 110 gsl_linalg_QR_solve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_vector * b, gsl_vector * x)
 111 {
 112   if (QR->size1 != QR->size2)
 113     {
 114       GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
 115     }
 116   else if (QR->size1 != b->size)
 117     {
 118       GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
 119     }
 120   else if (QR->size2 != x->size)
 121     {
 122       GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
 123     }
 124   else
 125     {
 126       /* Copy x <- b */
 127
 128       gsl_vector_memcpy (x, b);
 129
 130       /* Solve for x */
 131
 132       gsl_linalg_QR_svx (QR, tau, x);
 133
 134       return GSL_SUCCESS;
 135     }
 136 }
 137
 138 /* Solves the system A x = b in place using the QR factorisation,
 139
 140  *  R x = Q^T b
 141  *
 142  * to obtain x. Based on SLATEC code.
 143  */
 144
 145 int
 146 gsl_linalg_QR_svx (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * x)
 147 {
 148
 149   if (QR->size1 != QR->size2)
 150     {
 151       GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
 152     }
 153   else if (QR->size1 != x->size)
 154     {
 155       GSL_ERROR ("matrix size must match x/rhs size", GSL_EBADLEN);
 156     }
 157   else
 158     {
 159       /* compute rhs = Q^T b */
 160
 161       gsl_linalg_QR_QTvec (QR, tau, x);
 162
 163       /* Solve R x = rhs, storing x in-place */
 164
 165       gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, QR, x);
 166
 167       return GSL_SUCCESS;
 168     }
 169 }
 170
 171
 172 /* Find the least squares solution to the overdetermined system
 173  *
 174  *   A x = b
 175  *
 176  * for M >= N using the QR factorization A = Q R.
 177  */
 178
 179 int
 180 gsl_linalg_QR_lssolve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_vector * b, gsl_vector * x, gsl_vector * residual)
 181 {
 182   const size_t M = QR->size1;
 183   const size_t N = QR->size2;
 184
 185   if (M < N)
 186     {
 187       GSL_ERROR ("QR matrix must have M>=N", GSL_EBADLEN);
 188     }
 189   else if (M != b->size)
 190     {
 191       GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
 192     }
 193   else if (N != x->size)
 194     {
 195       GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
 196     }
 197   else if (M != residual->size)
 198     {
 199       GSL_ERROR ("matrix size must match residual size", GSL_EBADLEN);
 200     }
 201   else
 202     {
 203       gsl_matrix_const_view R = gsl_matrix_const_submatrix (QR, 0, 0, N, N);
 204       gsl_vector_view c = gsl_vector_subvector(residual, 0, N);
 205
 206       gsl_vector_memcpy(residual, b);
 207
 208       /* compute rhs = Q^T b */
 209
 210       gsl_linalg_QR_QTvec (QR, tau, residual);
 211
 212       /* Solve R x = rhs */
 213
 214       gsl_vector_memcpy(x, &(c.vector));
 215
 216       gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, &(R.matrix), x);
 217
 218       /* Compute residual = b - A x = Q (Q^T b - R x) */
 219
 220       gsl_vector_set_zero(&(c.vector));
 221
 222       gsl_linalg_QR_Qvec(QR, tau, residual);
 223
 224       return GSL_SUCCESS;
 225     }
 226 }
 227
 228
 229 int
 230 gsl_linalg_QR_Rsolve (const gsl_matrix * QR, const gsl_vector * b, gsl_vector * x)
 231 {
 232   if (QR->size1 != QR->size2)
 233     {
 234       GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
 235     }
 236   else if (QR->size1 != b->size)
 237     {
 238       GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
 239     }
 240   else if (QR->size2 != x->size)
 241     {
 242       GSL_ERROR ("matrix size must match x size", GSL_EBADLEN);
 243     }
 244   else
 245     {
 246       /* Copy x <- b */
 247
 248       gsl_vector_memcpy (x, b);
 249
 250       /* Solve R x = b, storing x in-place */
 251
 252       gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, QR, x);
 253
 254       return GSL_SUCCESS;
 255     }
 256 }
 257
 258
 259 int
 260 gsl_linalg_QR_Rsvx (const gsl_matrix * QR, gsl_vector * x)
 261 {
 262   if (QR->size1 != QR->size2)
 263     {
 264       GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
 265     }
 266   else if (QR->size1 != x->size)
 267     {
 268       GSL_ERROR ("matrix size must match rhs size", GSL_EBADLEN);
 269     }
 270   else
 271     {
 272       /* Solve R x = b, storing x in-place */
 273
 274       gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, QR, x);
 275
 276       return GSL_SUCCESS;
 277     }
 278 }
 279
 280 int
 281 gsl_linalg_R_solve (const gsl_matrix * R, const gsl_vector * b, gsl_vector * x)
 282 {
 283   if (R->size1 != R->size2)
 284     {
 285       GSL_ERROR ("R matrix must be square", GSL_ENOTSQR);
 286     }
 287   else if (R->size1 != b->size)
 288     {
 289       GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
 290     }
 291   else if (R->size2 != x->size)
 292     {
 293       GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
 294     }
 295   else
 296     {
 297       /* Copy x <- b */
 298
 299       gsl_vector_memcpy (x, b);
 300
 301       /* Solve R x = b, storing x inplace in b */
 302
 303       gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, R, x);
 304
 305       return GSL_SUCCESS;
 306     }
 307 }
 308
 309 int
 310 gsl_linalg_R_svx (const gsl_matrix * R, gsl_vector * x)
 311 {
 312   if (R->size1 != R->size2)
 313     {
 314       GSL_ERROR ("R matrix must be square", GSL_ENOTSQR);
 315     }
 316   else if (R->size2 != x->size)
 317     {
 318       GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
 319     }
 320   else
 321     {
 322       /* Solve R x = b, storing x inplace in b */
 323
 324       gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, R, x);
 325
 326       return GSL_SUCCESS;
 327     }
 328 }
 329
 330
 331
 332 /* Form the product Q^T v  from a QR factorized matrix
 333  */
 334
 335 int
 336 gsl_linalg_QR_QTvec (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * v)
 337 {
 338   const size_t M = QR->size1;
 339   const size_t N = QR->size2;
 340
 341   if (tau->size != GSL_MIN (M, N))
 342     {
 343       GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
 344     }
 345   else if (v->size != M)
 346     {
 347       GSL_ERROR ("vector size must be N", GSL_EBADLEN);
 348     }
 349   else
 350     {
 351       size_t i;
 352
 353       /* compute Q^T v */
 354
 355       for (i = 0; i < GSL_MIN (M, N); i++)
 356         {
 357           gsl_vector_const_view c = gsl_matrix_const_column (QR, i);
 358           gsl_vector_const_view h = gsl_vector_const_subvector (&(c.vector), i, M - i);
 359           gsl_vector_view w = gsl_vector_subvector (v, i, M - i);
 360           double ti = gsl_vector_get (tau, i);
 361           gsl_linalg_householder_hv (ti, &(h.vector), &(w.vector));
 362         }
 363       return GSL_SUCCESS;
 364     }
 365 }
 366
 367
 368 int
 369 gsl_linalg_QR_Qvec (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * v)
 370 {
 371   const size_t M = QR->size1;
 372   const size_t N = QR->size2;
 373
 374   if (tau->size != GSL_MIN (M, N))
 375     {
 376       GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
 377     }
 378   else if (v->size != M)
 379     {
 380       GSL_ERROR ("vector size must be N", GSL_EBADLEN);
 381     }
 382   else
 383     {
 384       size_t i;
 385
 386       /* compute Q^T v */
 387
 388       for (i = GSL_MIN (M, N); i > 0 && i--;)
 389         {
 390           gsl_vector_const_view c = gsl_matrix_const_column (QR, i);
 391           gsl_vector_const_view h = gsl_vector_const_subvector (&(c.vector),
 392                                                                 i, M - i);
 393           gsl_vector_view w = gsl_vector_subvector (v, i, M - i);
 394           double ti = gsl_vector_get (tau, i);
 395           gsl_linalg_householder_hv (ti, &h.vector, &w.vector);
 396         }
 397       return GSL_SUCCESS;
 398     }
 399 }
 400
 401 /* Form the product Q^T A  from a QR factorized matrix */
 402
 403 int
 404 gsl_linalg_QR_QTmat (const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * A)
 405 {
 406   const size_t M = QR->size1;
 407   const size_t N = QR->size2;
 408
 409   if (tau->size != GSL_MIN (M, N))
 410     {
 411       GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
 412     }
 413   else if (A->size1 != M)
 414     {
 415       GSL_ERROR ("matrix must have M rows", GSL_EBADLEN);
 416     }
 417   else
 418     {
 419       size_t i;
 420
 421       /* compute Q^T A */
 422
 423       for (i = 0; i < GSL_MIN (M, N); i++)
 424         {
 425           gsl_vector_const_view c = gsl_matrix_const_column (QR, i);
 426           gsl_vector_const_view h = gsl_vector_const_subvector (&(c.vector), i, M - i);
 427           gsl_matrix_view m = gsl_matrix_submatrix(A, i, 0, M - i, A->size2);
 428           double ti = gsl_vector_get (tau, i);
 429           gsl_linalg_householder_hm (ti, &(h.vector), &(m.matrix));
 430         }
 431       return GSL_SUCCESS;
 432     }
 433 }
 434
 435
 436 /*  Form the orthogonal matrix Q from the packed QR matrix */
 437
 438 int
 439 gsl_linalg_QR_unpack (const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * Q, gsl_matrix * R)
 440 {
 441   const size_t M = QR->size1;
 442   const size_t N = QR->size2;
 443
 444   if (Q->size1 != M || Q->size2 != M)
 445     {
 446       GSL_ERROR ("Q matrix must be M x M", GSL_ENOTSQR);
 447     }
 448   else if (R->size1 != M || R->size2 != N)
 449     {
 450       GSL_ERROR ("R matrix must be M x N", GSL_ENOTSQR);
 451     }
 452   else if (tau->size != GSL_MIN (M, N))
 453     {
 454       GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
 455     }
 456   else
 457     {
 458       size_t i, j;
 459
 460       /* Initialize Q to the identity */
 461
 462       gsl_matrix_set_identity (Q);
 463
 464       for (i = GSL_MIN (M, N); i > 0 && i--;)
 465         {
 466           gsl_vector_const_view c = gsl_matrix_const_column (QR, i);
 467           gsl_vector_const_view h = gsl_vector_const_subvector (&c.vector,
 468                                                                 i, M - i);
 469           gsl_matrix_view m = gsl_matrix_submatrix (Q, i, i, M - i, M - i);
 470           double ti = gsl_vector_get (tau, i);
 471           gsl_linalg_householder_hm (ti, &h.vector, &m.matrix);
 472         }
 473
 474       /*  Form the right triangular matrix R from a packed QR matrix */
 475
 476       for (i = 0; i < M; i++)
 477         {
 478           for (j = 0; j < i && j < N; j++)
 479             gsl_matrix_set (R, i, j, 0.0);
 480
 481           for (j = i; j < N; j++)
 482             gsl_matrix_set (R, i, j, gsl_matrix_get (QR, i, j));
 483         }
 484
 485       return GSL_SUCCESS;
 486     }
 487 }
 488
 489
 490 /* Update a QR factorisation for A= Q R ,  A' = A + u v^T,
 491
 492  * Q' R' = QR + u v^T
 493  *       = Q (R + Q^T u v^T)
 494  *       = Q (R + w v^T)
 495  *
 496  * where w = Q^T u.
 497  *
 498  * Algorithm from Golub and Van Loan, "Matrix Computations", Section
 499  * 12.5 (Updating Matrix Factorizations, Rank-One Changes)
 500  */
 501
 502 int
 503 gsl_linalg_QR_update (gsl_matrix * Q, gsl_matrix * R,
 504                       gsl_vector * w, const gsl_vector * v)
 505 {
 506   const size_t M = R->size1;
 507   const size_t N = R->size2;
 508
 509   if (Q->size1 != M || Q->size2 != M)
 510     {
 511       GSL_ERROR ("Q matrix must be M x M if R is M x N", GSL_ENOTSQR);
 512     }
 513   else if (w->size != M)
 514     {
 515       GSL_ERROR ("w must be length M if R is M x N", GSL_EBADLEN);
 516     }
 517   else if (v->size != N)
 518     {
 519       GSL_ERROR ("v must be length N if R is M x N", GSL_EBADLEN);
 520     }
 521   else
 522     {
 523       size_t j, k;
 524       double w0;
 525
 526       /* Apply Given's rotations to reduce w to (|w|, 0, 0, ... , 0)
 527
 528          J_1^T .... J_(n-1)^T w = +/- |w| e_1
 529
 530          simultaneously applied to R,  H = J_1^T ... J^T_(n-1) R
 531          so that H is upper Hessenberg.  (12.5.2) */
 532
 533       for (k = M - 1; k > 0; k--)
 534         {
 535           double c, s;
 536           double wk = gsl_vector_get (w, k);
 537           double wkm1 = gsl_vector_get (w, k - 1);
 538
 539           create_givens (wkm1, wk, &c, &s);
 540           apply_givens_vec (w, k - 1, k, c, s);
 541           apply_givens_qr (M, N, Q, R, k - 1, k, c, s);
 542         }
 543
 544       w0 = gsl_vector_get (w, 0);
 545
 546       /* Add in w v^T  (Equation 12.5.3) */
 547
 548       for (j = 0; j < N; j++)
 549         {
 550           double r0j = gsl_matrix_get (R, 0, j);
 551           double vj = gsl_vector_get (v, j);
 552           gsl_matrix_set (R, 0, j, r0j + w0 * vj);
 553         }
 554
 555       /* Apply Givens transformations R' = G_(n-1)^T ... G_1^T H
 556          Equation 12.5.4 */
 557
 558       for (k = 1; k < GSL_MIN(M,N+1); k++)
 559         {
 560           double c, s;
 561           double diag = gsl_matrix_get (R, k - 1, k - 1);
 562           double offdiag = gsl_matrix_get (R, k, k - 1);
 563
 564           create_givens (diag, offdiag, &c, &s);
 565           apply_givens_qr (M, N, Q, R, k - 1, k, c, s);
 566
 567           gsl_matrix_set (R, k, k - 1, 0.0);    /* exact zero of G^T */
 568         }
 569
 570       return GSL_SUCCESS;
 571     }
 572 }
 573
 574 int
 575 gsl_linalg_QR_QRsolve (gsl_matrix * Q, gsl_matrix * R, const gsl_vector * b, gsl_vector * x)
 576 {
 577   const size_t M = R->size1;
 578   const size_t N = R->size2;
 579
 580   if (M != N)
 581     {
 582       return GSL_ENOTSQR;
 583     }
 584   else if (Q->size1 != M || b->size != M || x->size != M)
 585     {
 586       return GSL_EBADLEN;
 587     }
 588   else
 589     {
 590       /* compute sol = Q^T b */
 591
 592       gsl_blas_dgemv (CblasTrans, 1.0, Q, b, 0.0, x);
 593
 594       /* Solve R x = sol, storing x in-place */
 595
 596       gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, R, x);
 597
 598       return GSL_SUCCESS;
 599     }
 600 }