|
Packit |
67cb25 |
/* linalg/householder.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004, 2007, 2010 Gerard Jungman, Brian Gough
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_vector.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_matrix.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_blas.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <gsl/gsl_linalg.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_linalg_householder_transform()
|
|
Packit |
67cb25 |
Compute a householder transformation (tau,v) of a vector
|
|
Packit |
67cb25 |
x so that P x = [ I - tau*v*v' ] x annihilates x(1:n-1)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: v - on input, x vector
|
|
Packit |
67cb25 |
on output, householder vector v
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Notes:
|
|
Packit |
67cb25 |
1) on output, v is normalized so that v[0] = 1. The 1 is
|
|
Packit |
67cb25 |
not actually stored; instead v[0] = -sign(x[0])*||x|| so
|
|
Packit |
67cb25 |
that:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
P x = v[0] * e_1
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Therefore external routines should take care when applying
|
|
Packit |
67cb25 |
the projection matrix P to vectors, taking into account
|
|
Packit |
67cb25 |
that v[0] should be 1 when doing so.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_linalg_householder_transform (gsl_vector * v)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* replace v[0:n-1] with a householder vector (v[0:n-1]) and
|
|
Packit |
67cb25 |
coefficient tau that annihilate v[1:n-1] */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const size_t n = v->size ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (n == 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return 0.0; /* tau = 0 */
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double alpha, beta, tau ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_view x = gsl_vector_subvector (v, 1, n - 1) ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double xnorm = gsl_blas_dnrm2 (&x.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (xnorm == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return 0.0; /* tau = 0 */
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
alpha = gsl_vector_get (v, 0) ;
|
|
Packit |
67cb25 |
beta = - (alpha >= 0.0 ? +1.0 : -1.0) * hypot(alpha, xnorm) ;
|
|
Packit |
67cb25 |
tau = (beta - alpha) / beta ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double s = (alpha - beta);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (fabs(s) > GSL_DBL_MIN)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_blas_dscal (1.0 / s, &x.vector);
|
|
Packit |
67cb25 |
gsl_vector_set (v, 0, beta) ;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_blas_dscal (GSL_DBL_EPSILON / s, &x.vector);
|
|
Packit |
67cb25 |
gsl_blas_dscal (1.0 / GSL_DBL_EPSILON, &x.vector);
|
|
Packit |
67cb25 |
gsl_vector_set (v, 0, beta) ;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return tau;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_householder_hm (double tau, const gsl_vector * v, gsl_matrix * A)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* applies a householder transformation v,tau to matrix m */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (tau == 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#ifdef USE_BLAS
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_const_view v1 = gsl_vector_const_subvector (v, 1, v->size - 1);
|
|
Packit |
67cb25 |
gsl_matrix_view A1 = gsl_matrix_submatrix (A, 1, 0, A->size1 - 1, A->size2);
|
|
Packit |
67cb25 |
size_t j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < A->size2; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double wj = 0.0;
|
|
Packit |
67cb25 |
gsl_vector_view A1j = gsl_matrix_column(&A1.matrix, j);
|
|
Packit |
67cb25 |
gsl_blas_ddot (&A1j.vector, &v1.vector, &wj;;
|
|
Packit |
67cb25 |
wj += gsl_matrix_get(A,0,j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double A0j = gsl_matrix_get (A, 0, j);
|
|
Packit |
67cb25 |
gsl_matrix_set (A, 0, j, A0j - tau * wj);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_daxpy (-tau * wj, &v1.vector, &A1j.vector);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
#else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < A->size2; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Compute wj = Akj vk */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double wj = gsl_matrix_get(A,0,j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 1; i < A->size1; i++) /* note, computed for v(0) = 1 above */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
wj += gsl_matrix_get(A,i,j) * gsl_vector_get(v,i);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Aij = Aij - tau vi wj */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* i = 0 */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double A0j = gsl_matrix_get (A, 0, j);
|
|
Packit |
67cb25 |
gsl_matrix_set (A, 0, j, A0j - tau * wj);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* i = 1 .. M-1 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 1; i < A->size1; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double Aij = gsl_matrix_get (A, i, j);
|
|
Packit |
67cb25 |
double vi = gsl_vector_get (v, i);
|
|
Packit |
67cb25 |
gsl_matrix_set (A, i, j, Aij - tau * vi * wj);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
#endif
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_householder_mh (double tau, const gsl_vector * v, gsl_matrix * A)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* applies a householder transformation v,tau to matrix m from the
|
|
Packit |
67cb25 |
right hand side in order to zero out rows */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (tau == 0)
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* A = A - tau w v' */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#ifdef USE_BLAS
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_const_view v1 = gsl_vector_const_subvector (v, 1, v->size - 1);
|
|
Packit |
67cb25 |
gsl_matrix_view A1 = gsl_matrix_submatrix (A, 0, 1, A->size1, A->size2-1);
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < A->size1; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double wi = 0.0;
|
|
Packit |
67cb25 |
gsl_vector_view A1i = gsl_matrix_row(&A1.matrix, i);
|
|
Packit |
67cb25 |
gsl_blas_ddot (&A1i.vector, &v1.vector, &wi;;
|
|
Packit |
67cb25 |
wi += gsl_matrix_get(A,i,0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double Ai0 = gsl_matrix_get (A, i, 0);
|
|
Packit |
67cb25 |
gsl_matrix_set (A, i, 0, Ai0 - tau * wi);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_daxpy(-tau * wi, &v1.vector, &A1i.vector);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
#else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < A->size1; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double wi = gsl_matrix_get(A,i,0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 1; j < A->size2; j++) /* note, computed for v(0) = 1 above */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
wi += gsl_matrix_get(A,i,j) * gsl_vector_get(v,j);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* j = 0 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double Ai0 = gsl_matrix_get (A, i, 0);
|
|
Packit |
67cb25 |
gsl_matrix_set (A, i, 0, Ai0 - tau * wi);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* j = 1 .. N-1 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 1; j < A->size2; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double vj = gsl_vector_get (v, j);
|
|
Packit |
67cb25 |
double Aij = gsl_matrix_get (A, i, j);
|
|
Packit |
67cb25 |
gsl_matrix_set (A, i, j, Aij - tau * wi * vj);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
#endif
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_householder_hv (double tau, const gsl_vector * v, gsl_vector * w)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* applies a householder transformation v to vector w */
|
|
Packit |
67cb25 |
const size_t N = v->size;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (tau == 0)
|
|
Packit |
67cb25 |
return GSL_SUCCESS ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* compute d = v'w */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double w0 = gsl_vector_get(w,0);
|
|
Packit |
67cb25 |
double d1, d;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_const_view v1 = gsl_vector_const_subvector(v, 1, N-1);
|
|
Packit |
67cb25 |
gsl_vector_view w1 = gsl_vector_subvector(w, 1, N-1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute d1 = v(2:n)'w(2:n) */
|
|
Packit |
67cb25 |
gsl_blas_ddot (&v1.vector, &w1.vector, &d1;;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute d = v'w = w(1) + d1 since v(1) = 1 */
|
|
Packit |
67cb25 |
d = w0 + d1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute w = w - tau (v) (v'w) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set (w, 0, w0 - tau * d);
|
|
Packit |
67cb25 |
gsl_blas_daxpy (-tau * d, &v1.vector, &w1.vector);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_householder_hm1 (double tau, gsl_matrix * A)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* applies a householder transformation v,tau to a matrix being
|
|
Packit |
67cb25 |
build up from the identity matrix, using the first column of A as
|
|
Packit |
67cb25 |
a householder vector */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (tau == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i,j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_set (A, 0, 0, 1.0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 1; j < A->size2; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_set (A, 0, j, 0.0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 1; i < A->size1; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_set (A, i, 0, 0.0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* w = A' v */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#ifdef USE_BLAS
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_view A1 = gsl_matrix_submatrix (A, 1, 0, A->size1 - 1, A->size2);
|
|
Packit |
67cb25 |
gsl_vector_view v1 = gsl_matrix_column (&A1.matrix, 0);
|
|
Packit |
67cb25 |
size_t j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 1; j < A->size2; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double wj = 0.0; /* A0j * v0 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_view A1j = gsl_matrix_column(&A1.matrix, j);
|
|
Packit |
67cb25 |
gsl_blas_ddot (&A1j.vector, &v1.vector, &wj;;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* A = A - tau v w' */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_set (A, 0, j, - tau * wj);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_daxpy(-tau*wj, &v1.vector, &A1j.vector);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_dscal(-tau, &v1.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_set (A, 0, 0, 1.0 - tau);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
#else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 1; j < A->size2; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double wj = 0.0; /* A0j * v0 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 1; i < A->size1; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double vi = gsl_matrix_get(A, i, 0);
|
|
Packit |
67cb25 |
wj += gsl_matrix_get(A,i,j) * vi;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* A = A - tau v w' */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_set (A, 0, j, - tau * wj);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 1; i < A->size1; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double vi = gsl_matrix_get (A, i, 0);
|
|
Packit |
67cb25 |
double Aij = gsl_matrix_get (A, i, j);
|
|
Packit |
67cb25 |
gsl_matrix_set (A, i, j, Aij - tau * vi * wj);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 1; i < A->size1; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double vi = gsl_matrix_get(A, i, 0);
|
|
Packit |
67cb25 |
gsl_matrix_set(A, i, 0, -tau * vi);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_set (A, 0, 0, 1.0 - tau);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
#endif
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|