summaryrefslogtreecommitdiff
path: root/gsl-1.9/linalg/householder.c
diff options
context:
space:
mode:
Diffstat (limited to 'gsl-1.9/linalg/householder.c')
-rw-r--r--gsl-1.9/linalg/householder.c326
1 files changed, 326 insertions, 0 deletions
diff --git a/gsl-1.9/linalg/householder.c b/gsl-1.9/linalg/householder.c
new file mode 100644
index 0000000..d9ab952
--- /dev/null
+++ b/gsl-1.9/linalg/householder.c
@@ -0,0 +1,326 @@
+/* linalg/householder.c
+ *
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 Gerard Jungman, Brian Gough
+ *
+ * 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.
+ */
+
+#include <config.h>
+#include <stdlib.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_blas.h>
+
+#include <gsl/gsl_linalg.h>
+
+double
+gsl_linalg_householder_transform (gsl_vector * v)
+{
+ /* replace v[0:n-1] with a householder vector (v[0:n-1]) and
+ coefficient tau that annihilate v[1:n-1] */
+
+ const size_t n = v->size ;
+
+ if (n == 1)
+ {
+ return 0.0; /* tau = 0 */
+ }
+ else
+ {
+ double alpha, beta, tau ;
+
+ gsl_vector_view x = gsl_vector_subvector (v, 1, n - 1) ;
+
+ double xnorm = gsl_blas_dnrm2 (&x.vector);
+
+ if (xnorm == 0)
+ {
+ return 0.0; /* tau = 0 */
+ }
+
+ alpha = gsl_vector_get (v, 0) ;
+ beta = - (alpha >= 0.0 ? +1.0 : -1.0) * hypot(alpha, xnorm) ;
+ tau = (beta - alpha) / beta ;
+
+ gsl_blas_dscal (1.0 / (alpha - beta), &x.vector);
+ gsl_vector_set (v, 0, beta) ;
+
+ return tau;
+ }
+}
+
+int
+gsl_linalg_householder_hm (double tau, const gsl_vector * v, gsl_matrix * A)
+{
+ /* applies a householder transformation v,tau to matrix m */
+
+ if (tau == 0.0)
+ {
+ return GSL_SUCCESS;
+ }
+
+#ifdef USE_BLAS
+ {
+ gsl_vector_const_view v1 = gsl_vector_const_subvector (v, 1, v->size - 1);
+ gsl_matrix_view A1 = gsl_matrix_submatrix (A, 1, 0, A->size1 - 1, A->size2);
+ size_t j;
+
+ for (j = 0; j < A->size2; j++)
+ {
+ double wj = 0.0;
+ gsl_vector_view A1j = gsl_matrix_column(&A1.matrix, j);
+ gsl_blas_ddot (&A1j.vector, &v1.vector, &wj);
+ wj += gsl_matrix_get(A,0,j);
+
+ {
+ double A0j = gsl_matrix_get (A, 0, j);
+ gsl_matrix_set (A, 0, j, A0j - tau * wj);
+ }
+
+ gsl_blas_daxpy (-tau * wj, &v1.vector, &A1j.vector);
+ }
+ }
+#else
+ {
+ size_t i, j;
+
+ for (j = 0; j < A->size2; j++)
+ {
+ /* Compute wj = Akj vk */
+
+ double wj = gsl_matrix_get(A,0,j);
+
+ for (i = 1; i < A->size1; i++) /* note, computed for v(0) = 1 above */
+ {
+ wj += gsl_matrix_get(A,i,j) * gsl_vector_get(v,i);
+ }
+
+ /* Aij = Aij - tau vi wj */
+
+ /* i = 0 */
+ {
+ double A0j = gsl_matrix_get (A, 0, j);
+ gsl_matrix_set (A, 0, j, A0j - tau * wj);
+ }
+
+ /* i = 1 .. M-1 */
+
+ for (i = 1; i < A->size1; i++)
+ {
+ double Aij = gsl_matrix_get (A, i, j);
+ double vi = gsl_vector_get (v, i);
+ gsl_matrix_set (A, i, j, Aij - tau * vi * wj);
+ }
+ }
+ }
+#endif
+
+ return GSL_SUCCESS;
+}
+
+int
+gsl_linalg_householder_mh (double tau, const gsl_vector * v, gsl_matrix * A)
+{
+ /* applies a householder transformation v,tau to matrix m from the
+ right hand side in order to zero out rows */
+
+ if (tau == 0)
+ return GSL_SUCCESS;
+
+ /* A = A - tau w v' */
+
+#ifdef USE_BLAS
+ {
+ gsl_vector_const_view v1 = gsl_vector_const_subvector (v, 1, v->size - 1);
+ gsl_matrix_view A1 = gsl_matrix_submatrix (A, 0, 1, A->size1, A->size2-1);
+ size_t i;
+
+ for (i = 0; i < A->size1; i++)
+ {
+ double wi = 0.0;
+ gsl_vector_view A1i = gsl_matrix_row(&A1.matrix, i);
+ gsl_blas_ddot (&A1i.vector, &v1.vector, &wi);
+ wi += gsl_matrix_get(A,i,0);
+
+ {
+ double Ai0 = gsl_matrix_get (A, i, 0);
+ gsl_matrix_set (A, i, 0, Ai0 - tau * wi);
+ }
+
+ gsl_blas_daxpy(-tau * wi, &v1.vector, &A1i.vector);
+ }
+ }
+#else
+ {
+ size_t i, j;
+
+ for (i = 0; i < A->size1; i++)
+ {
+ double wi = gsl_matrix_get(A,i,0);
+
+ for (j = 1; j < A->size2; j++) /* note, computed for v(0) = 1 above */
+ {
+ wi += gsl_matrix_get(A,i,j) * gsl_vector_get(v,j);
+ }
+
+ /* j = 0 */
+
+ {
+ double Ai0 = gsl_matrix_get (A, i, 0);
+ gsl_matrix_set (A, i, 0, Ai0 - tau * wi);
+ }
+
+ /* j = 1 .. N-1 */
+
+ for (j = 1; j < A->size2; j++)
+ {
+ double vj = gsl_vector_get (v, j);
+ double Aij = gsl_matrix_get (A, i, j);
+ gsl_matrix_set (A, i, j, Aij - tau * wi * vj);
+ }
+ }
+ }
+#endif
+
+ return GSL_SUCCESS;
+}
+
+int
+gsl_linalg_householder_hv (double tau, const gsl_vector * v, gsl_vector * w)
+{
+ /* applies a householder transformation v to vector w */
+ const size_t N = v->size;
+
+ if (tau == 0)
+ return GSL_SUCCESS ;
+
+ {
+ /* compute d = v'w */
+
+ double d0 = gsl_vector_get(w,0);
+ double d1, d;
+
+ gsl_vector_const_view v1 = gsl_vector_const_subvector(v, 1, N-1);
+ gsl_vector_view w1 = gsl_vector_subvector(w, 1, N-1);
+
+ gsl_blas_ddot (&v1.vector, &w1.vector, &d1);
+
+ d = d0 + d1;
+
+ /* compute w = w - tau (v) (v'w) */
+
+ {
+ double w0 = gsl_vector_get (w,0);
+ gsl_vector_set (w, 0, w0 - tau * d);
+ }
+
+ gsl_blas_daxpy (-tau * d, &v1.vector, &w1.vector);
+ }
+
+ return GSL_SUCCESS;
+}
+
+
+int
+gsl_linalg_householder_hm1 (double tau, gsl_matrix * A)
+{
+ /* applies a householder transformation v,tau to a matrix being
+ build up from the identity matrix, using the first column of A as
+ a householder vector */
+
+ if (tau == 0)
+ {
+ size_t i,j;
+
+ gsl_matrix_set (A, 0, 0, 1.0);
+
+ for (j = 1; j < A->size2; j++)
+ {
+ gsl_matrix_set (A, 0, j, 0.0);
+ }
+
+ for (i = 1; i < A->size1; i++)
+ {
+ gsl_matrix_set (A, i, 0, 0.0);
+ }
+
+ return GSL_SUCCESS;
+ }
+
+ /* w = A' v */
+
+#ifdef USE_BLAS
+ {
+ gsl_matrix_view A1 = gsl_matrix_submatrix (A, 1, 0, A->size1 - 1, A->size2);
+ gsl_vector_view v1 = gsl_matrix_column (&A1.matrix, 0);
+ size_t j;
+
+ for (j = 1; j < A->size2; j++)
+ {
+ double wj = 0.0; /* A0j * v0 */
+
+ gsl_vector_view A1j = gsl_matrix_column(&A1.matrix, j);
+ gsl_blas_ddot (&A1j.vector, &v1.vector, &wj);
+
+ /* A = A - tau v w' */
+
+ gsl_matrix_set (A, 0, j, - tau * wj);
+
+ gsl_blas_daxpy(-tau*wj, &v1.vector, &A1j.vector);
+ }
+
+ gsl_blas_dscal(-tau, &v1.vector);
+
+ gsl_matrix_set (A, 0, 0, 1.0 - tau);
+ }
+#else
+ {
+ size_t i, j;
+
+ for (j = 1; j < A->size2; j++)
+ {
+ double wj = 0.0; /* A0j * v0 */
+
+ for (i = 1; i < A->size1; i++)
+ {
+ double vi = gsl_matrix_get(A, i, 0);
+ wj += gsl_matrix_get(A,i,j) * vi;
+ }
+
+ /* A = A - tau v w' */
+
+ gsl_matrix_set (A, 0, j, - tau * wj);
+
+ for (i = 1; i < A->size1; i++)
+ {
+ double vi = gsl_matrix_get (A, i, 0);
+ double Aij = gsl_matrix_get (A, i, j);
+ gsl_matrix_set (A, i, j, Aij - tau * vi * wj);
+ }
+ }
+
+ for (i = 1; i < A->size1; i++)
+ {
+ double vi = gsl_matrix_get(A, i, 0);
+ gsl_matrix_set(A, i, 0, -tau * vi);
+ }
+
+ gsl_matrix_set (A, 0, 0, 1.0 - tau);
+ }
+#endif
+
+ return GSL_SUCCESS;
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-01 02:35:45 +0000