|
Packit |
67cb25 |
/* linalg/ptlq.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
|
|
Packit |
67cb25 |
* Copyright (C) 2004 Joerg Wensch, modifications for LQ.
|
|
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 <string.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_blas.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_permute_vector.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_linalg.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "apply_givens.c"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* The purpose of this package is to speed up QR-decomposition for
|
|
Packit |
67cb25 |
large matrices. Because QR-decomposition is column oriented, but
|
|
Packit |
67cb25 |
GSL uses a row-oriented matrix format, there can considerable
|
|
Packit |
67cb25 |
speedup obtained by computing the LQ-decomposition of the
|
|
Packit |
67cb25 |
transposed matrix instead. This package provides LQ-decomposition
|
|
Packit |
67cb25 |
and related algorithms. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Factorise a general N x M matrix A into
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* P A = L Q
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* where Q is orthogonal (M x M) and L is lower triangular (N x M).
|
|
Packit |
67cb25 |
* When A is rank deficient, r = rank(A) < n, then the permutation is
|
|
Packit |
67cb25 |
* used to ensure that the lower n - r columns of L are zero and the first
|
|
Packit |
67cb25 |
* l rows of Q form an orthonormal basis for the rows of A.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Q is stored as a packed set of Householder transformations in the
|
|
Packit |
67cb25 |
* strict upper triangular part of the input matrix.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* L is stored in the diagonal and lower triangle of the input matrix.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* P: column j of P is column k of the identity matrix, where k =
|
|
Packit |
67cb25 |
* permutation->data[j]
|
|
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,i+1), m(i,i+2), ... , m(i,M)]
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This storage scheme is the same as in LAPACK. See LAPACK's
|
|
Packit |
67cb25 |
* dgeqpf.f for details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_PTLQ_decomp (gsl_matrix * A, gsl_vector * tau, gsl_permutation * p, int *signum, gsl_vector * norm)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = A->size1;
|
|
Packit |
67cb25 |
const size_t M = 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 if (p->size != N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("permutation size must be N", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (norm->size != N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("norm size must be N", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*signum = 1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_permutation_init (p); /* set to identity */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute column norms and store in workspace */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_view c = gsl_matrix_row (A, i);
|
|
Packit |
67cb25 |
double x = gsl_blas_dnrm2 (&c.vector);
|
|
Packit |
67cb25 |
gsl_vector_set (norm, i, x);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < GSL_MIN (M, N); i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Bring the column of largest norm into the pivot position */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double max_norm = gsl_vector_get(norm, i);
|
|
Packit |
67cb25 |
size_t j, kmax = i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = i + 1; j < N; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double x = gsl_vector_get (norm, j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (x > max_norm)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
max_norm = x;
|
|
Packit |
67cb25 |
kmax = j;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (kmax != i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_swap_rows (A, i, kmax);
|
|
Packit |
67cb25 |
gsl_permutation_swap (p, i, kmax);
|
|
Packit |
67cb25 |
gsl_vector_swap_elements(norm,i,kmax);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
(*signum) = -(*signum);
|
|
Packit |
67cb25 |
}
|
|
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 |
{
|
|
Packit |
67cb25 |
gsl_vector_view c_full = gsl_matrix_row (A, i);
|
|
Packit |
67cb25 |
gsl_vector_view c = gsl_vector_subvector (&c_full.vector,
|
|
Packit |
67cb25 |
i, M - i);
|
|
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 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (i + 1 < N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_view m = gsl_matrix_submatrix (A, i +1, i, N - (i+1), M - i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_linalg_householder_mh (tau_i, &c.vector, &m.matrix);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Update the norms of the remaining columns too */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (i + 1 < M)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (j = i + 1; j < N; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double x = gsl_vector_get (norm, j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (x > 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double y = 0;
|
|
Packit |
67cb25 |
double temp= gsl_matrix_get (A, j, i) / x;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (fabs (temp) >= 1)
|
|
Packit |
67cb25 |
y = 0.0;
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
y = x * sqrt (1 - temp * temp);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* recompute norm to prevent loss of accuracy */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (fabs (y / x) < sqrt (20.0) * GSL_SQRT_DBL_EPSILON)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_view c_full = gsl_matrix_row (A, j);
|
|
Packit |
67cb25 |
gsl_vector_view c =
|
|
Packit |
67cb25 |
gsl_vector_subvector(&c_full.vector,
|
|
Packit |
67cb25 |
i+1, M - (i+1));
|
|
Packit |
67cb25 |
y = gsl_blas_dnrm2 (&c.vector);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set (norm, j, y);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_PTLQ_decomp2 (const gsl_matrix * A, gsl_matrix * q, gsl_matrix * r, gsl_vector * tau, gsl_permutation * p, int *signum, gsl_vector * norm)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = A->size1;
|
|
Packit |
67cb25 |
const size_t M = A->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (q->size1 != M || q->size2 !=M)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("q must be M x M", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (r->size1 != N || r->size2 !=M)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("r must be N x M", GSL_EBADLEN);
|
|
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 if (p->size != N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("permutation size must be N", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (norm->size != N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("norm size must be N", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_memcpy (r, A);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_linalg_PTLQ_decomp (r, tau, p, signum, norm);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* FIXME: aliased arguments depends on behavior of unpack routine! */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_linalg_LQ_unpack (r, tau, q, r);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Solves the system x^T A = b^T using the P^T L Q factorisation,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
z^T L = b^T Q^T
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
x = P z;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
to obtain x. Based on SLATEC code. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_PTLQ_solve_T (const gsl_matrix * QR,
|
|
Packit |
67cb25 |
const gsl_vector * tau,
|
|
Packit |
67cb25 |
const gsl_permutation * p,
|
|
Packit |
67cb25 |
const gsl_vector * b,
|
|
Packit |
67cb25 |
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->size2 != p->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match permutation size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (QR->size2 != b->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (QR->size1 != 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 |
gsl_vector_memcpy (x, b);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_linalg_PTLQ_svx_T (QR, tau, p, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_PTLQ_svx_T (const gsl_matrix * LQ,
|
|
Packit |
67cb25 |
const gsl_vector * tau,
|
|
Packit |
67cb25 |
const gsl_permutation * p,
|
|
Packit |
67cb25 |
gsl_vector * x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (LQ->size1 != LQ->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("LQ matrix must be square", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (LQ->size2 != p->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match permutation size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (LQ->size1 != 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 |
/* compute sol = b^T Q^T */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_linalg_LQ_vecQT (LQ, tau, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Solve L^T x = sol, storing x inplace in sol */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, LQ, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_permute_vector_inverse (p, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_PTLQ_LQsolve_T (const gsl_matrix * Q, const gsl_matrix * L,
|
|
Packit |
67cb25 |
const gsl_permutation * p,
|
|
Packit |
67cb25 |
const gsl_vector * b,
|
|
Packit |
67cb25 |
gsl_vector * x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (Q->size1 != Q->size2 || L->size1 != L->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return GSL_ENOTSQR;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (Q->size1 != p->size || Q->size1 != L->size1
|
|
Packit |
67cb25 |
|| Q->size1 != b->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return GSL_EBADLEN;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* compute b' = Q b */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_dgemv (CblasNoTrans, 1.0, Q, b, 0.0, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Solve L^T x = b', storing x inplace */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, L, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Apply permutation to solution in place */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_permute_vector_inverse (p, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_PTLQ_Lsolve_T (const gsl_matrix * LQ,
|
|
Packit |
67cb25 |
const gsl_permutation * p,
|
|
Packit |
67cb25 |
const gsl_vector * b,
|
|
Packit |
67cb25 |
gsl_vector * x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (LQ->size1 != LQ->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("LQ matrix must be square", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (LQ->size1 != b->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (LQ->size2 != x->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match x size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (p->size != x->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("permutation 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 L^T x = b, storing x inplace */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, LQ, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_permute_vector_inverse (p, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_PTLQ_Lsvx_T (const gsl_matrix * LQ,
|
|
Packit |
67cb25 |
const gsl_permutation * p,
|
|
Packit |
67cb25 |
gsl_vector * x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (LQ->size1 != LQ->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("LQ matrix must be square", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (LQ->size2 != x->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match x size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (p->size != x->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("permutation size must match x size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Solve L^T x = b, storing x inplace */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, LQ, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_permute_vector_inverse (p, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Update a P^T L Q factorisation for P A= L Q , A' = A + v u^T,
|
|
Packit |
67cb25 |
PA' = PA + Pv u^T
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* P^T L' Q' = P^T LQ + v u^T
|
|
Packit |
67cb25 |
* = P^T (L + (P v) u^T Q^T) Q
|
|
Packit |
67cb25 |
* = P^T (L + (P v) w^T) Q
|
|
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_PTLQ_update (gsl_matrix * Q, gsl_matrix * L,
|
|
Packit |
67cb25 |
const gsl_permutation * p,
|
|
Packit |
67cb25 |
const gsl_vector * v, gsl_vector * w)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (Q->size1 != Q->size2 || L->size1 != L->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return GSL_ENOTSQR;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (L->size1 != Q->size2 || v->size != Q->size2 || w->size != Q->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return GSL_EBADLEN;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t j, k;
|
|
Packit |
67cb25 |
const size_t N = Q->size1;
|
|
Packit |
67cb25 |
const size_t M = Q->size2;
|
|
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 L, H = J_1^T ... J^T_(n-1) L
|
|
Packit |
67cb25 |
so that H is upper Hessenberg. (12.5.2) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (k = M - 1; k > 0; k--)
|
|
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_lq (M, N, Q, L, k - 1, k, c, s);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
w0 = gsl_vector_get (w, 0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Add in v w^T (Equation 12.5.3) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < N; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double lj0 = gsl_matrix_get (L, j, 0);
|
|
Packit |
67cb25 |
size_t p_j = gsl_permutation_get (p, j);
|
|
Packit |
67cb25 |
double vj = gsl_vector_get (v, p_j);
|
|
Packit |
67cb25 |
gsl_matrix_set (L, j, 0, lj0 + w0 * vj);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Apply Givens transformations L' = G_(n-1)^T ... G_1^T H
|
|
Packit |
67cb25 |
Equation 12.5.4 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (k = 1; k < N; k++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double c, s;
|
|
Packit |
67cb25 |
double diag = gsl_matrix_get (L, k - 1, k - 1);
|
|
Packit |
67cb25 |
double offdiag = gsl_matrix_get (L, k - 1, k );
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_linalg_givens (diag, offdiag, &c, &s);
|
|
Packit |
67cb25 |
apply_givens_lq (M, N, Q, L, k - 1, k, c, s);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|