diff options
Diffstat (limited to 'gsl-1.9/linalg/apply_givens.c')
| -rw-r--r-- | gsl-1.9/linalg/apply_givens.c | 125 |
1 files changed, 125 insertions, 0 deletions
diff --git a/gsl-1.9/linalg/apply_givens.c b/gsl-1.9/linalg/apply_givens.c new file mode 100644 index 0000000..2d4fe55 --- /dev/null +++ b/gsl-1.9/linalg/apply_givens.c @@ -0,0 +1,125 @@ +/* linalg/apply_givens.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001 Gerard Jungman, Brian Gough + * Copyright (C) 2004 Joerg Wensch, modifications for LQ. + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +inline static void +apply_givens_qr (size_t M, size_t N, gsl_matrix * Q, gsl_matrix * R, + size_t i, size_t j, double c, double s) +{ + size_t k; + + /* Apply rotation to matrix Q, Q' = Q G */ + +#if USE_BLAS + { + gsl_matrix_view Q0M = gsl_matrix_submatrix(Q,0,0,M,j+1); + gsl_vector_view Qi = gsl_matrix_column(&Q0M.matrix,i); + gsl_vector_view Qj = gsl_matrix_column(&Q0M.matrix,j); + gsl_blas_drot(&Qi.vector, &Qj.vector, c, -s); + } +#else + for (k = 0; k < M; k++) + { + double qki = gsl_matrix_get (Q, k, i); + double qkj = gsl_matrix_get (Q, k, j); + gsl_matrix_set (Q, k, i, qki * c - qkj * s); + gsl_matrix_set (Q, k, j, qki * s + qkj * c); + } +#endif + + /* Apply rotation to matrix R, R' = G^T R (note: upper triangular so + zero for column < row) */ + +#if USE_BLAS + { + k = GSL_MIN(i,j); + gsl_matrix_view R0 = gsl_matrix_submatrix(R, 0, k, j+1, N-k); + gsl_vector_view Ri = gsl_matrix_row(&R0.matrix,i); + gsl_vector_view Rj = gsl_matrix_row(&R0.matrix,j); + gsl_blas_drot(&Ri.vector, &Rj.vector, c, -s); + } +#else + for (k = GSL_MIN (i, j); k < N; k++) + { + double rik = gsl_matrix_get (R, i, k); + double rjk = gsl_matrix_get (R, j, k); + gsl_matrix_set (R, i, k, c * rik - s * rjk); + gsl_matrix_set (R, j, k, s * rik + c * rjk); + } +#endif +} + +inline static void +apply_givens_lq (size_t M, size_t N, gsl_matrix * Q, gsl_matrix * L, + size_t i, size_t j, double c, double s) +{ + size_t k; + + /* Apply rotation to matrix Q, Q' = G Q */ + +#if USE_BLAS + { + gsl_matrix_view Q0M = gsl_matrix_submatrix(Q,0,0,j+1,M); + gsl_vector_view Qi = gsl_matrix_row(&Q0M.matrix,i); + gsl_vector_view Qj = gsl_matrix_row(&Q0M.matrix,j); + gsl_blas_drot(&Qi.vector, &Qj.vector, c, -s); + } +#else + for (k = 0; k < M; k++) + { + double qik = gsl_matrix_get (Q, i, k); + double qjk = gsl_matrix_get (Q, j, k); + gsl_matrix_set (Q, i, k, qik * c - qjk * s); + gsl_matrix_set (Q, j, k, qik * s + qjk * c); + } +#endif + + /* Apply rotation to matrix L, L' = L G^T (note: lower triangular so + zero for column > row) */ + +#if USE_BLAS + { + k = GSL_MIN(i,j); + gsl_matrix_view L0 = gsl_matrix_submatrix(L, k, 0, N-k, j+1); + gsl_vector_view Li = gsl_matrix_column(&L0.matrix,i); + gsl_vector_view Lj = gsl_matrix_column(&L0.matrix,j); + gsl_blas_drot(&Li.vector, &Lj.vector, c, -s); + } +#else + for (k = GSL_MIN (i, j); k < N; k++) + { + double lki = gsl_matrix_get (L, k, i); + double lkj = gsl_matrix_get (L, k, j); + gsl_matrix_set (L, k, i, c * lki - s * lkj); + gsl_matrix_set (L, k, j, s * lki + c * lkj); + } +#endif +} + +inline static void +apply_givens_vec (gsl_vector * v, size_t i, size_t j, double c, double s) +{ + /* Apply rotation to vector v' = G^T v */ + + double vi = gsl_vector_get (v, i); + double vj = gsl_vector_get (v, j); + gsl_vector_set (v, i, c * vi - s * vj); + gsl_vector_set (v, j, s * vi + c * vj); +} + |