|
Packit |
67cb25 |
/* linalg/qr.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 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 |
/* Author: G. Jungman */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <string.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_linalg.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 "apply_givens.c"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Factorise a general M x N matrix A into
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* A = Q R
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* where Q is orthogonal (M x M) and R is upper triangular (M x N).
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Q is stored as a packed set of Householder transformations in the
|
|
Packit |
67cb25 |
* strict lower triangular part of the input matrix.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* R is stored in the diagonal and upper triangle of the input matrix.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* The full matrix for Q can be obtained as the product
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Q = Q_k .. Q_2 Q_1
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* where k = MIN(M,N) and
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Q_i = (I - tau_i * v_i * v_i')
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* and where v_i is a Householder vector
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* v_i = [1, m(i+1,i), m(i+2,i), ... , m(M,i)]
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This storage scheme is the same as in LAPACK. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_QR_decomp (gsl_matrix * A, gsl_vector * tau)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t M = A->size1;
|
|
Packit |
67cb25 |
const size_t N = A->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (tau->size != GSL_MIN (M, N))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < GSL_MIN (M, N); i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Compute the Householder transformation to reduce the j-th
|
|
Packit |
67cb25 |
column of the matrix to a multiple of the j-th unit vector */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_view c_full = gsl_matrix_column (A, i);
|
|
Packit |
67cb25 |
gsl_vector_view c = gsl_vector_subvector (&(c_full.vector), i, M-i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double tau_i = gsl_linalg_householder_transform (&(c.vector));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set (tau, i, tau_i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Apply the transformation to the remaining columns and
|
|
Packit |
67cb25 |
update the norms */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (i + 1 < N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_view m = gsl_matrix_submatrix (A, i, i + 1, M - i, N - (i + 1));
|
|
Packit |
67cb25 |
gsl_linalg_householder_hm (tau_i, &(c.vector), &(m.matrix));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Solves the system A x = b using the QR factorisation,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* R x = Q^T b
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* to obtain x. Based on SLATEC code.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_QR_solve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_vector * b, gsl_vector * x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (QR->size1 != QR->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (QR->size1 != b->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (QR->size2 != x->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Copy x <- b */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_memcpy (x, b);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Solve for x */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_linalg_QR_svx (QR, tau, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Solves the system A x = b in place using the QR factorisation,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* R x = Q^T b
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* to obtain x. Based on SLATEC code.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_QR_svx (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (QR->size1 != QR->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (QR->size1 != x->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match x/rhs size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* compute rhs = Q^T b */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_linalg_QR_QTvec (QR, tau, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Solve R x = rhs, storing x in-place */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, QR, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Find the least squares solution to the overdetermined system
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* A x = b
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* for M >= N using the QR factorization A = Q R.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_QR_lssolve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_vector * b, gsl_vector * x, gsl_vector * residual)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t M = QR->size1;
|
|
Packit |
67cb25 |
const size_t N = QR->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (M < N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("QR matrix must have M>=N", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (M != b->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (N != x->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (M != residual->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match residual size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_const_view R = gsl_matrix_const_submatrix (QR, 0, 0, N, N);
|
|
Packit |
67cb25 |
gsl_vector_view c = gsl_vector_subvector(residual, 0, N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_memcpy(residual, b);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute rhs = Q^T b */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_linalg_QR_QTvec (QR, tau, residual);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Solve R x = rhs */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_memcpy(x, &(c.vector));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, &(R.matrix), x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute residual = b - A x = Q (Q^T b - R x) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set_zero(&(c.vector));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_linalg_QR_Qvec(QR, tau, residual);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_QR_Rsolve (const gsl_matrix * QR, const gsl_vector * b, gsl_vector * x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (QR->size1 != QR->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (QR->size1 != b->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (QR->size2 != x->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match x size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Copy x <- b */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_memcpy (x, b);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Solve R x = b, storing x in-place */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, QR, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_QR_Rsvx (const gsl_matrix * QR, gsl_vector * x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (QR->size1 != QR->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (QR->size1 != x->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match rhs size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Solve R x = b, storing x in-place */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, QR, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_R_solve (const gsl_matrix * R, const gsl_vector * b, gsl_vector * x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (R->size1 != R->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("R matrix must be square", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (R->size1 != b->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (R->size2 != x->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Copy x <- b */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_memcpy (x, b);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Solve R x = b, storing x inplace in b */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, R, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_R_svx (const gsl_matrix * R, gsl_vector * x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (R->size1 != R->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("R matrix must be square", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (R->size2 != x->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Solve R x = b, storing x inplace in b */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, R, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Form the product Q^T v from a QR factorized matrix
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_QR_QTvec (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * v)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t M = QR->size1;
|
|
Packit |
67cb25 |
const size_t N = QR->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (tau->size != GSL_MIN (M, N))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (v->size != M)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("vector size must be M", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute Q^T v */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < GSL_MIN (M, N); i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_const_view c = gsl_matrix_const_column (QR, i);
|
|
Packit |
67cb25 |
gsl_vector_const_view h = gsl_vector_const_subvector (&(c.vector), i, M - i);
|
|
Packit |
67cb25 |
gsl_vector_view w = gsl_vector_subvector (v, i, M - i);
|
|
Packit |
67cb25 |
double ti = gsl_vector_get (tau, i);
|
|
Packit |
67cb25 |
gsl_linalg_householder_hv (ti, &(h.vector), &(w.vector));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_QR_Qvec (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * v)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t M = QR->size1;
|
|
Packit |
67cb25 |
const size_t N = QR->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (tau->size != GSL_MIN (M, N))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (v->size != M)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("vector size must be M", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute Q v */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = GSL_MIN (M, N); i-- > 0;)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_const_view c = gsl_matrix_const_column (QR, i);
|
|
Packit |
67cb25 |
gsl_vector_const_view h = gsl_vector_const_subvector (&(c.vector),
|
|
Packit |
67cb25 |
i, M - i);
|
|
Packit |
67cb25 |
gsl_vector_view w = gsl_vector_subvector (v, i, M - i);
|
|
Packit |
67cb25 |
double ti = gsl_vector_get (tau, i);
|
|
Packit |
67cb25 |
gsl_linalg_householder_hv (ti, &h.vector, &w.vector);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Form the product Q^T A from a QR factorized matrix */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_QR_QTmat (const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * A)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t M = QR->size1;
|
|
Packit |
67cb25 |
const size_t N = QR->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (tau->size != GSL_MIN (M, N))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (A->size1 != M)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix must have M rows", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute Q^T A */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < GSL_MIN (M, N); i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_const_view c = gsl_matrix_const_column (QR, i);
|
|
Packit |
67cb25 |
gsl_vector_const_view h = gsl_vector_const_subvector (&(c.vector), i, M - i);
|
|
Packit |
67cb25 |
gsl_matrix_view m = gsl_matrix_submatrix(A, i, 0, M - i, A->size2);
|
|
Packit |
67cb25 |
double ti = gsl_vector_get (tau, i);
|
|
Packit |
67cb25 |
gsl_linalg_householder_hm (ti, &(h.vector), &(m.matrix));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Form the product A Q from a QR factorized matrix */
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_QR_matQ (const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * A)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t M = QR->size1;
|
|
Packit |
67cb25 |
const size_t N = QR->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (tau->size != GSL_MIN (M, N))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (A->size2 != M)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix must have M columns", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute A Q */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < GSL_MIN (M, N); i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_const_view c = gsl_matrix_const_column (QR, i);
|
|
Packit |
67cb25 |
gsl_vector_const_view h = gsl_vector_const_subvector (&(c.vector), i, M - i);
|
|
Packit |
67cb25 |
gsl_matrix_view m = gsl_matrix_submatrix(A, 0, i, A->size1, M - i);
|
|
Packit |
67cb25 |
double ti = gsl_vector_get (tau, i);
|
|
Packit |
67cb25 |
gsl_linalg_householder_mh (ti, &(h.vector), &(m.matrix));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Form the orthogonal matrix Q from the packed QR matrix */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_QR_unpack (const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * Q, gsl_matrix * R)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t M = QR->size1;
|
|
Packit |
67cb25 |
const size_t N = QR->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (Q->size1 != M || Q->size2 != M)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("Q matrix must be M x M", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (R->size1 != M || R->size2 != N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("R matrix must be M x N", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (tau->size != GSL_MIN (M, N))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Initialize Q to the identity */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_set_identity (Q);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = GSL_MIN (M, N); i-- > 0;)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_const_view c = gsl_matrix_const_column (QR, i);
|
|
Packit |
67cb25 |
gsl_vector_const_view h = gsl_vector_const_subvector (&c.vector,
|
|
Packit |
67cb25 |
i, M - i);
|
|
Packit |
67cb25 |
gsl_matrix_view m = gsl_matrix_submatrix (Q, i, i, M - i, M - i);
|
|
Packit |
67cb25 |
double ti = gsl_vector_get (tau, i);
|
|
Packit |
67cb25 |
gsl_linalg_householder_hm (ti, &h.vector, &m.matrix);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Form the right triangular matrix R from a packed QR matrix */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < M; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (j = 0; j < i && j < N; j++)
|
|
Packit |
67cb25 |
gsl_matrix_set (R, i, j, 0.0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = i; j < N; j++)
|
|
Packit |
67cb25 |
gsl_matrix_set (R, i, j, gsl_matrix_get (QR, i, j));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Update a QR factorisation for A= Q R , A' = A + u v^T,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* Q' R' = QR + u v^T
|
|
Packit |
67cb25 |
* = Q (R + Q^T u v^T)
|
|
Packit |
67cb25 |
* = Q (R + w v^T)
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* where w = Q^T u.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Algorithm from Golub and Van Loan, "Matrix Computations", Section
|
|
Packit |
67cb25 |
* 12.5 (Updating Matrix Factorizations, Rank-One Changes)
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_QR_update (gsl_matrix * Q, gsl_matrix * R,
|
|
Packit |
67cb25 |
gsl_vector * w, const gsl_vector * v)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t M = R->size1;
|
|
Packit |
67cb25 |
const size_t N = R->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (Q->size1 != M || Q->size2 != M)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("Q matrix must be M x M if R is M x N", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (w->size != M)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("w must be length M if R is M x N", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (v->size != N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("v must be length N if R is M x N", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t j, k;
|
|
Packit |
67cb25 |
double w0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Apply Given's rotations to reduce w to (|w|, 0, 0, ... , 0)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
J_1^T .... J_(n-1)^T w = +/- |w| e_1
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
simultaneously applied to R, H = J_1^T ... J^T_(n-1) R
|
|
Packit |
67cb25 |
so that H is upper Hessenberg. (12.5.2) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (k = M - 1; k > 0; k--) /* loop from k = M-1 to 1 */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double c, s;
|
|
Packit |
67cb25 |
double wk = gsl_vector_get (w, k);
|
|
Packit |
67cb25 |
double wkm1 = gsl_vector_get (w, k - 1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_linalg_givens (wkm1, wk, &c, &s);
|
|
Packit |
67cb25 |
gsl_linalg_givens_gv (w, k - 1, k, c, s);
|
|
Packit |
67cb25 |
apply_givens_qr (M, N, Q, R, k - 1, k, c, s);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
w0 = gsl_vector_get (w, 0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Add in w v^T (Equation 12.5.3) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < N; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double r0j = gsl_matrix_get (R, 0, j);
|
|
Packit |
67cb25 |
double vj = gsl_vector_get (v, j);
|
|
Packit |
67cb25 |
gsl_matrix_set (R, 0, j, r0j + w0 * vj);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Apply Givens transformations R' = G_(n-1)^T ... G_1^T H
|
|
Packit |
67cb25 |
Equation 12.5.4 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (k = 1; k < GSL_MIN(M,N+1); k++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double c, s;
|
|
Packit |
67cb25 |
double diag = gsl_matrix_get (R, k - 1, k - 1);
|
|
Packit |
67cb25 |
double offdiag = gsl_matrix_get (R, k, k - 1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_linalg_givens (diag, offdiag, &c, &s);
|
|
Packit |
67cb25 |
apply_givens_qr (M, N, Q, R, k - 1, k, c, s);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_set (R, k, k - 1, 0.0); /* exact zero of G^T */
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_QR_QRsolve (gsl_matrix * Q, gsl_matrix * R, const gsl_vector * b, gsl_vector * x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t M = R->size1;
|
|
Packit |
67cb25 |
const size_t N = R->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (M != N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return GSL_ENOTSQR;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (Q->size1 != M || b->size != M || x->size != M)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return GSL_EBADLEN;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* compute sol = Q^T b */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_dgemv (CblasTrans, 1.0, Q, b, 0.0, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Solve R x = sol, storing x in-place */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, R, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|