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

   1 /* linalg/apply_givens.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 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 inline static void
  22 apply_givens_qr (size_t M, size_t N, gsl_matrix * Q, gsl_matrix * R,
  23                  size_t i, size_t j, double c, double s)
  24 {
  25   size_t k;
  26
  27   /* Apply rotation to matrix Q,  Q' = Q G */
  28
  29 #if USE_BLAS
  30   {
  31     gsl_matrix_view Q0M = gsl_matrix_submatrix(Q,0,0,M,j+1);
  32     gsl_vector_view Qi = gsl_matrix_column(&Q0M.matrix,i);
  33     gsl_vector_view Qj = gsl_matrix_column(&Q0M.matrix,j);
  34     gsl_blas_drot(&Qi.vector, &Qj.vector, c, -s);
  35   }
  36 #else
  37   for (k = 0; k < M; k++)
  38     {
  39       double qki = gsl_matrix_get (Q, k, i);
  40       double qkj = gsl_matrix_get (Q, k, j);
  41       gsl_matrix_set (Q, k, i, qki * c - qkj * s);
  42       gsl_matrix_set (Q, k, j, qki * s + qkj * c);
  43     }
  44 #endif
  45
  46   /* Apply rotation to matrix R, R' = G^T R (note: upper triangular so
  47      zero for column < row) */
  48
  49 #if USE_BLAS
  50   {
  51     k = GSL_MIN(i,j);
  52     gsl_matrix_view R0 = gsl_matrix_submatrix(R, 0, k, j+1, N-k);
  53     gsl_vector_view Ri = gsl_matrix_row(&R0.matrix,i);
  54     gsl_vector_view Rj = gsl_matrix_row(&R0.matrix,j);
  55     gsl_blas_drot(&Ri.vector, &Rj.vector, c, -s);
  56   }
  57 #else
  58   for (k = GSL_MIN (i, j); k < N; k++)
  59     {
  60       double rik = gsl_matrix_get (R, i, k);
  61       double rjk = gsl_matrix_get (R, j, k);
  62       gsl_matrix_set (R, i, k, c * rik - s * rjk);
  63       gsl_matrix_set (R, j, k, s * rik + c * rjk);
  64     }
  65 #endif
  66 }
  67
  68 inline static void
  69 apply_givens_lq (size_t M, size_t N, gsl_matrix * Q, gsl_matrix * L,
  70                  size_t i, size_t j, double c, double s)
  71 {
  72   size_t k;
  73
  74   /* Apply rotation to matrix Q,  Q' = G Q */
  75
  76 #if USE_BLAS
  77   {
  78     gsl_matrix_view Q0M = gsl_matrix_submatrix(Q,0,0,j+1,M);
  79     gsl_vector_view Qi = gsl_matrix_row(&Q0M.matrix,i);
  80     gsl_vector_view Qj = gsl_matrix_row(&Q0M.matrix,j);
  81     gsl_blas_drot(&Qi.vector, &Qj.vector, c, -s);
  82   }
  83 #else
  84   for (k = 0; k < M; k++)
  85     {
  86       double qik = gsl_matrix_get (Q, i, k);
  87       double qjk = gsl_matrix_get (Q, j, k);
  88       gsl_matrix_set (Q, i, k, qik * c - qjk * s);
  89       gsl_matrix_set (Q, j, k, qik * s + qjk * c);
  90     }
  91 #endif
  92
  93   /* Apply rotation to matrix L, L' = L G^T (note: lower triangular so
  94      zero for column > row) */
  95
  96 #if USE_BLAS
  97   {
  98     k = GSL_MIN(i,j);
  99     gsl_matrix_view L0 = gsl_matrix_submatrix(L, k, 0, N-k, j+1);
 100     gsl_vector_view Li = gsl_matrix_column(&L0.matrix,i);
 101     gsl_vector_view Lj = gsl_matrix_column(&L0.matrix,j);
 102     gsl_blas_drot(&Li.vector, &Lj.vector, c, -s);
 103   }
 104 #else
 105   for (k = GSL_MIN (i, j); k < N; k++)
 106     {
 107       double lki = gsl_matrix_get (L, k, i);
 108       double lkj = gsl_matrix_get (L, k, j);
 109       gsl_matrix_set (L, k, i, c * lki - s * lkj);
 110       gsl_matrix_set (L, k, j, s * lki + c * lkj);
 111     }
 112 #endif
 113 }
 114
 115 inline static void
 116 apply_givens_vec (gsl_vector * v, size_t i, size_t j, double c, double s)
 117 {
 118   /* Apply rotation to vector v' = G^T v */
 119
 120   double vi = gsl_vector_get (v, i);
 121   double vj = gsl_vector_get (v, j);
 122   gsl_vector_set (v, i, c * vi - s * vj);
 123   gsl_vector_set (v, j, s * vi + c * vj);
 124 }
 125