From b329da7d951d0c63080b730d859830c203cda04e Mon Sep 17 00:00:00 2001 From: Christian J Krueger Date: Wed, 3 Jun 2020 12:27:00 +0200 Subject: [PATCH 2/2] add complex QR decomposition * File qr.c copied into qrc.c and adjusted all functions to allow for complex values (except gsl_linalg_QR_update() as this also requires complex-valued Givens rotations) * Implement test cases for complex valued QR decomposition based on the tests for the real-valued version. --- linalg/Makefile.am | 2 +- linalg/gsl_linalg.h | 70 ++++++ linalg/qrc.c | 644 ++++++++++++++++++++++++++++++++++++++++++++++++++++ linalg/test.c | 351 ++++++++++++++++++++++++++++ 4 files changed, 1066 insertions(+), 1 deletion(-) create mode 100644 linalg/qrc.c diff --git a/linalg/Makefile.am b/linalg/Makefile.am index c7fc58ce..cdc630d9 100644 --- a/linalg/Makefile.am +++ b/linalg/Makefile.am @@ -4,7 +4,7 @@ pkginclude_HEADERS = gsl_linalg.h AM_CPPFLAGS = -I$(top_srcdir) -libgsllinalg_la_SOURCES = cod.c condest.c invtri.c invtri_complex.c multiply.c exponential.c tridiag.c tridiag.h lu.c luc.c hh.c ql.c qr.c qrpt.c qr_tr.c rqr.c lq.c ptlq.c svd.c householder.c householdercomplex.c hessenberg.c hesstri.c cholesky.c choleskyc.c mcholesky.c pcholesky.c cholesky_band.c ldlt.c ldlt_band.c symmtd.c hermtd.c bidiag.c balance.c balancemat.c inline.c trimult.c trimult_complex.c +libgsllinalg_la_SOURCES = cod.c condest.c invtri.c invtri_complex.c multiply.c exponential.c tridiag.c tridiag.h lu.c luc.c hh.c ql.c qr.c qrc.c qrpt.c qr_tr.c rqr.c lq.c ptlq.c svd.c householder.c householdercomplex.c hessenberg.c hesstri.c cholesky.c choleskyc.c mcholesky.c pcholesky.c cholesky_band.c ldlt.c ldlt_band.c symmtd.c hermtd.c bidiag.c balance.c balancemat.c inline.c trimult.c trimult_complex.c noinst_HEADERS = apply_givens.c cholesky_common.c recurse.h svdstep.c tridiag.h test_cholesky.c test_choleskyc.c test_cod.c test_common.c test_ldlt.c test_lu.c test_luc.c test_lq.c test_ql.c test_qr.c test_tri.c diff --git a/linalg/gsl_linalg.h b/linalg/gsl_linalg.h index 6267942d..4a1e0e45 100644 --- a/linalg/gsl_linalg.h +++ b/linalg/gsl_linalg.h @@ -298,6 +298,76 @@ int gsl_linalg_R_svx (const gsl_matrix * R, gsl_vector * x); int gsl_linalg_QR_rcond(const gsl_matrix * QR, double * rcond, gsl_vector * work); +/* complex QR decomposition */ + +int gsl_linalg_complex_QR_decomp (gsl_matrix_complex * A, + gsl_vector_complex * tau); + + +int gsl_linalg_complex_QR_solve (const gsl_matrix_complex * QR, + const gsl_vector_complex * tau, + const gsl_vector_complex * b, + gsl_vector_complex * x); + + +int gsl_linalg_complex_QR_svx (const gsl_matrix_complex * QR, + const gsl_vector_complex * tau, + gsl_vector_complex * x); + +int gsl_linalg_complex_QR_lssolve (const gsl_matrix_complex * QR, + const gsl_vector_complex * tau, + const gsl_vector_complex * b, + gsl_vector_complex * x, + gsl_vector_complex * residual); + +int gsl_linalg_complex_QR_QRsolve (gsl_matrix_complex * Q, + gsl_matrix_complex * R, + const gsl_vector_complex * b, + gsl_vector_complex * x); + +int gsl_linalg_complex_QR_Rsolve (const gsl_matrix_complex * QR, + const gsl_vector_complex * b, + gsl_vector_complex * x); + +int gsl_linalg_complex_QR_Rsvx (const gsl_matrix_complex * QR, + gsl_vector_complex * x); + +/* +int gsl_linalg_QR_update (gsl_matrix * Q, + gsl_matrix * R, + gsl_vector * w, + const gsl_vector * v); +*/ + +int gsl_linalg_complex_QR_QTvec (const gsl_matrix_complex * QR, + const gsl_vector_complex * tau, + gsl_vector_complex * v); + +int gsl_linalg_complex_QR_Qvec (const gsl_matrix_complex * QR, + const gsl_vector_complex * tau, + gsl_vector_complex * v); + +int gsl_linalg_complex_QR_QTmat (const gsl_matrix_complex * QR, + const gsl_vector_complex * tau, + gsl_matrix_complex * A); + +int gsl_linalg_complex_QR_matQ (const gsl_matrix_complex * QR, + const gsl_vector_complex * tau, + gsl_matrix_complex * A); + +int gsl_linalg_complex_QR_unpack (const gsl_matrix_complex * QR, + const gsl_vector_complex * tau, + gsl_matrix_complex * Q, + gsl_matrix_complex * R); + +int gsl_linalg_complex_R_solve (const gsl_matrix_complex * R, + const gsl_vector_complex * b, + gsl_vector_complex * x); + +int gsl_linalg_complex_R_svx (const gsl_matrix_complex * R, + gsl_vector_complex * x); + + /* Q R P^T decomposition */ int gsl_linalg_QRPT_decomp (gsl_matrix * A, diff --git a/linalg/qrc.c b/linalg/qrc.c new file mode 100644 index 00000000..2d2e2649 --- /dev/null +++ b/linalg/qrc.c @@ -0,0 +1,644 @@ +/* linalg/qrc.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough + * Copyright (C) 2017 Christian Krueger + * + * 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 3 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. + */ + +/* Author: G. Jungman, modified by C. Krueger */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "apply_givens.c" + +/* + * This code is essentially the same as in qr.c but adapted for complex + * valued matrices. Adapting the function gsl_linalg_QR_update() is + * a bit more involved and has not been done yet. + */ + +/* Factorise a general complex-valued M x N matrix A into + * + * A = Q R + * + * where Q is unitary (M x M) and R is upper triangular (M x N). + * + * Q is stored as a packed set of Householder transformations in the + * strict lower triangular part of the input matrix. + * + * R is stored in the diagonal and upper triangle of the input matrix. + * + * The full matrix for Q can be obtained as the product + * + * Q = Q_k .. Q_2 Q_1 + * + * where k = MIN(M,N) and + * + * Q_i = (I - tau_i * v_i * v_i') + * + * and where v_i is a Householder vector + * + * v_i = [1, m(i+1,i), m(i+2,i), ... , m(M,i)] + * + * This storage scheme is the same as in LAPACK. */ + +int +gsl_linalg_complex_QR_decomp (gsl_matrix_complex * A, gsl_vector_complex * tau) +{ + const size_t M = A->size1; + const size_t N = A->size2; + + if (tau->size != GSL_MIN (M, N)) + { + GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN); + } + else + { + size_t i; + + for (i = 0; i < GSL_MIN (M, N); i++) + { + /* Compute the Householder transformation to reduce the j-th + column of the matrix to a multiple of the j-th unit vector */ + + gsl_vector_complex_view c_full = gsl_matrix_complex_column (A, i); + gsl_vector_complex_view c = gsl_vector_complex_subvector (&(c_full.vector), i, M-i); + + gsl_complex tau_i = gsl_linalg_complex_householder_transform (&(c.vector)); + + gsl_vector_complex_set (tau, i, tau_i); + + /* Apply the transformation to the remaining columns and + update the norms */ + + if (i + 1 < N) + { + gsl_matrix_complex_view m = gsl_matrix_complex_submatrix (A, i, i + 1, M - i, N - (i + 1)); + gsl_complex tau_i_conj = gsl_complex_conjugate(tau_i); + gsl_linalg_complex_householder_hm (tau_i_conj, &(c.vector), &(m.matrix)); + } + } + + return GSL_SUCCESS; + } +} + +/* Solves the system A x = b using the QR factorisation, + + * R x = Q^* b + * + * to obtain x. Based on SLATEC code. + */ + +int +gsl_linalg_complex_QR_solve (const gsl_matrix_complex * QR, const gsl_vector_complex * tau, const gsl_vector_complex * b, gsl_vector_complex * x) +{ + if (QR->size1 != QR->size2) + { + GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR); + } + else if (QR->size1 != b->size) + { + GSL_ERROR ("matrix size must match b size", GSL_EBADLEN); + } + else if (QR->size2 != x->size) + { + GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN); + } + else + { + /* Copy x <- b */ + + gsl_vector_complex_memcpy (x, b); + + /* Solve for x */ + + gsl_linalg_complex_QR_svx (QR, tau, x); + + return GSL_SUCCESS; + } +} + +/* Solves the system A x = b in place using the QR factorisation, + + * R x = Q^* b + * + * to obtain x. Based on SLATEC code. + */ + +int +gsl_linalg_complex_QR_svx (const gsl_matrix_complex * QR, const gsl_vector_complex * tau, gsl_vector_complex * x) +{ + + if (QR->size1 != QR->size2) + { + GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR); + } + else if (QR->size1 != x->size) + { + GSL_ERROR ("matrix size must match x/rhs size", GSL_EBADLEN); + } + else + { + /* compute rhs = Q^* b */ + + gsl_linalg_complex_QR_QTvec (QR, tau, x); + + /* Solve R x = rhs, storing x in-place */ + + gsl_blas_ztrsv (CblasUpper, CblasNoTrans, CblasNonUnit, QR, x); + + return GSL_SUCCESS; + } +} + + +/* Find the least squares solution to the overdetermined system + * + * A x = b + * + * for M >= N using the QR factorization A = Q R. + */ + +int +gsl_linalg_complex_QR_lssolve (const gsl_matrix_complex * QR, const gsl_vector_complex * tau, const gsl_vector_complex * b, gsl_vector_complex * x, gsl_vector_complex * residual) +{ + const size_t M = QR->size1; + const size_t N = QR->size2; + + if (M < N) + { + GSL_ERROR ("QR matrix must have M>=N", GSL_EBADLEN); + } + else if (M != b->size) + { + GSL_ERROR ("matrix size must match b size", GSL_EBADLEN); + } + else if (N != x->size) + { + GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN); + } + else if (M != residual->size) + { + GSL_ERROR ("matrix size must match residual size", GSL_EBADLEN); + } + else + { + gsl_matrix_complex_const_view R = gsl_matrix_complex_const_submatrix (QR, 0, 0, N, N); + gsl_vector_complex_view c = gsl_vector_complex_subvector(residual, 0, N); + + gsl_vector_complex_memcpy(residual, b); + + /* compute rhs = Q^* b */ + + gsl_linalg_complex_QR_QTvec (QR, tau, residual); + + /* Solve R x = rhs */ + + gsl_vector_complex_memcpy(x, &(c.vector)); + + gsl_blas_ztrsv (CblasUpper, CblasNoTrans, CblasNonUnit, &(R.matrix), x); + + /* Compute residual = b - A x = Q (Q^* b - R x) */ + + gsl_vector_complex_set_zero(&(c.vector)); + + gsl_linalg_complex_QR_Qvec(QR, tau, residual); + + return GSL_SUCCESS; + } +} + + +int +gsl_linalg_complex_QR_Rsolve (const gsl_matrix_complex * QR, const gsl_vector_complex * b, gsl_vector_complex * x) +{ + if (QR->size1 != QR->size2) + { + GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR); + } + else if (QR->size1 != b->size) + { + GSL_ERROR ("matrix size must match b size", GSL_EBADLEN); + } + else if (QR->size2 != x->size) + { + GSL_ERROR ("matrix size must match x size", GSL_EBADLEN); + } + else + { + /* Copy x <- b */ + + gsl_vector_complex_memcpy (x, b); + + /* Solve R x = b, storing x in-place */ + + gsl_blas_ztrsv (CblasUpper, CblasNoTrans, CblasNonUnit, QR, x); + + return GSL_SUCCESS; + } +} + + +int +gsl_linalg_complex_QR_Rsvx (const gsl_matrix_complex * QR, gsl_vector_complex * x) +{ + if (QR->size1 != QR->size2) + { + GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR); + } + else if (QR->size1 != x->size) + { + GSL_ERROR ("matrix size must match rhs size", GSL_EBADLEN); + } + else + { + /* Solve R x = b, storing x in-place */ + + gsl_blas_ztrsv (CblasUpper, CblasNoTrans, CblasNonUnit, QR, x); + + return GSL_SUCCESS; + } +} + +int +gsl_linalg_complex_R_solve (const gsl_matrix_complex * R, const gsl_vector_complex * b, gsl_vector_complex * x) +{ + if (R->size1 != R->size2) + { + GSL_ERROR ("R matrix must be square", GSL_ENOTSQR); + } + else if (R->size1 != b->size) + { + GSL_ERROR ("matrix size must match b size", GSL_EBADLEN); + } + else if (R->size2 != x->size) + { + GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN); + } + else + { + /* Copy x <- b */ + + gsl_vector_complex_memcpy (x, b); + + /* Solve R x = b, storing x inplace in b */ + + gsl_blas_ztrsv (CblasUpper, CblasNoTrans, CblasNonUnit, R, x); + + return GSL_SUCCESS; + } +} + +int +gsl_linalg_complex_R_svx (const gsl_matrix_complex * R, gsl_vector_complex * x) +{ + if (R->size1 != R->size2) + { + GSL_ERROR ("R matrix must be square", GSL_ENOTSQR); + } + else if (R->size2 != x->size) + { + GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN); + } + else + { + /* Solve R x = b, storing x inplace in b */ + + gsl_blas_ztrsv (CblasUpper, CblasNoTrans, CblasNonUnit, R, x); + + return GSL_SUCCESS; + } +} + + +/* Form the product Q^* v from a QR factorized matrix + */ + +int +gsl_linalg_complex_QR_QTvec (const gsl_matrix_complex * QR, const gsl_vector_complex * tau, gsl_vector_complex * v) +{ + const size_t M = QR->size1; + const size_t N = QR->size2; + + if (tau->size != GSL_MIN (M, N)) + { + GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN); + } + else if (v->size != M) + { + GSL_ERROR ("vector size must be M", GSL_EBADLEN); + } + else + { + size_t i; + + /* compute Q^* v */ + + for (i = 0; i < GSL_MIN (M, N); i++) + { + gsl_vector_complex_const_view c = gsl_matrix_complex_const_column (QR, i); + gsl_vector_complex_const_view h = gsl_vector_complex_const_subvector (&(c.vector), i, M - i); + gsl_vector_complex_view w = gsl_vector_complex_subvector (v, i, M - i); + gsl_complex ti = gsl_vector_complex_get (tau, i); + gsl_complex ti_conj = gsl_complex_conjugate(ti); + gsl_linalg_complex_householder_hv (ti_conj, &(h.vector), &(w.vector)); + } + return GSL_SUCCESS; + } +} + + +int +gsl_linalg_complex_QR_Qvec (const gsl_matrix_complex * QR, const gsl_vector_complex * tau, gsl_vector_complex * v) +{ + const size_t M = QR->size1; + const size_t N = QR->size2; + + if (tau->size != GSL_MIN (M, N)) + { + GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN); + } + else if (v->size != M) + { + GSL_ERROR ("vector size must be M", GSL_EBADLEN); + } + else + { + size_t i; + + /* compute Q v */ + + for (i = GSL_MIN (M, N); i-- > 0;) + { + gsl_vector_complex_const_view c = gsl_matrix_complex_const_column (QR, i); + gsl_vector_complex_const_view h = gsl_vector_complex_const_subvector (&(c.vector), + i, M - i); + gsl_vector_complex_view w = gsl_vector_complex_subvector (v, i, M - i); + gsl_complex ti = gsl_vector_complex_get (tau, i); + /* we do not need the conjugate of ti here */ + gsl_linalg_complex_householder_hv (ti, &h.vector, &w.vector); + } + return GSL_SUCCESS; + } +} + +/* Form the product Q^* A from a QR factorized matrix */ + +int +gsl_linalg_complex_QR_QTmat (const gsl_matrix_complex * QR, const gsl_vector_complex * tau, gsl_matrix_complex * A) +{ + const size_t M = QR->size1; + const size_t N = QR->size2; + + if (tau->size != GSL_MIN (M, N)) + { + GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN); + } + else if (A->size1 != M) + { + GSL_ERROR ("matrix must have M rows", GSL_EBADLEN); + } + else + { + size_t i; + + /* compute Q^* A */ + + for (i = 0; i < GSL_MIN (M, N); i++) + { + gsl_vector_complex_const_view c = gsl_matrix_complex_const_column (QR, i); + gsl_vector_complex_const_view h = gsl_vector_complex_const_subvector (&(c.vector), i, M - i); + gsl_matrix_complex_view m = gsl_matrix_complex_submatrix(A, i, 0, M - i, A->size2); + gsl_complex ti = gsl_vector_complex_get (tau, i); + gsl_complex ti_conj = gsl_complex_conjugate(ti); + gsl_linalg_complex_householder_hm (ti_conj, &(h.vector), &(m.matrix)); + } + return GSL_SUCCESS; + } +} + +/* Form the product A Q from a QR factorized matrix */ +int +gsl_linalg_complex_QR_matQ (const gsl_matrix_complex * QR, const gsl_vector_complex * tau, gsl_matrix_complex * A) +{ + const size_t M = QR->size1; + const size_t N = QR->size2; + + if (tau->size != GSL_MIN (M, N)) + { + GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN); + } + else if (A->size2 != M) + { + GSL_ERROR ("matrix must have M columns", GSL_EBADLEN); + } + else + { + size_t i; + + /* compute A Q */ + + for (i = 0; i < GSL_MIN (M, N); i++) + { + gsl_vector_complex_const_view c = gsl_matrix_complex_const_column (QR, i); + gsl_vector_complex_const_view h = gsl_vector_complex_const_subvector (&(c.vector), i, M - i); + gsl_matrix_complex_view m = gsl_matrix_complex_submatrix(A, 0, i, A->size1, M - i); + gsl_complex ti = gsl_vector_complex_get (tau, i); + /* we do not need the conjugate of ti here */ + gsl_linalg_complex_householder_mh (ti, &(h.vector), &(m.matrix)); + } + return GSL_SUCCESS; + } +} + +/* Form the unitary matrix Q from the packed QR matrix */ + +int +gsl_linalg_complex_QR_unpack (const gsl_matrix_complex * QR, const gsl_vector_complex * tau, gsl_matrix_complex * Q, gsl_matrix_complex * R) +{ + const size_t M = QR->size1; + const size_t N = QR->size2; + + if (Q->size1 != M || Q->size2 != M) + { + GSL_ERROR ("Q matrix must be M x M", GSL_ENOTSQR); + } + else if (R->size1 != M || R->size2 != N) + { + GSL_ERROR ("R matrix must be M x N", GSL_ENOTSQR); + } + else if (tau->size != GSL_MIN (M, N)) + { + GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN); + } + else + { + size_t i, j; + + /* Initialize Q to the identity */ + + gsl_matrix_complex_set_identity (Q); + + for (i = GSL_MIN (M, N); i-- > 0;) + { + gsl_vector_complex_const_view c = gsl_matrix_complex_const_column (QR, i); + gsl_vector_complex_const_view h = gsl_vector_complex_const_subvector (&c.vector, i, M - i); + gsl_matrix_complex_view m = gsl_matrix_complex_submatrix (Q, i, i, M - i, M - i); + gsl_complex ti = gsl_vector_complex_get (tau, i); + /* we do not need the conjugate of ti here */ + gsl_linalg_complex_householder_hm (ti, &h.vector, &m.matrix); + } + + /* Form the right triangular matrix R from a packed QR matrix */ + + for (i = 0; i < M; i++) + { + for (j = 0; j < i && j < N; j++) + gsl_matrix_complex_set (R, i, j, GSL_COMPLEX_ZERO); + + for (j = i; j < N; j++) + gsl_matrix_complex_set (R, i, j, gsl_matrix_complex_get (QR, i, j)); + } + + return GSL_SUCCESS; + } +} +// +// +// TODO: This function still needs adapting to complex numbers. +// However, this also requires the functions supplied by GSL +// for the Givens rotations (gsl_linalg_givens_*) to be adapted. +// +// /* Update a QR factorisation for A= Q R , A' = A + u v^T, +// +// * Q' R' = QR + u v^T +// * = Q (R + Q^T u v^T) +// * = Q (R + w v^T) +// * +// * where w = Q^T u. +// * +// * Algorithm from Golub and Van Loan, "Matrix Computations", Section +// * 12.5 (Updating Matrix Factorizations, Rank-One Changes) +// */ +// +// int +// gsl_linalg_QR_update (gsl_matrix * Q, gsl_matrix * R, +// gsl_vector * w, const gsl_vector * v) +// { +// const size_t M = R->size1; +// const size_t N = R->size2; +// +// if (Q->size1 != M || Q->size2 != M) +// { +// GSL_ERROR ("Q matrix must be M x M if R is M x N", GSL_ENOTSQR); +// } +// else if (w->size != M) +// { +// GSL_ERROR ("w must be length M if R is M x N", GSL_EBADLEN); +// } +// else if (v->size != N) +// { +// GSL_ERROR ("v must be length N if R is M x N", GSL_EBADLEN); +// } +// else +// { +// size_t j, k; +// double w0; +// +// /* Apply Given's rotations to reduce w to (|w|, 0, 0, ... , 0) +// +// J_1^T .... J_(n-1)^T w = +/- |w| e_1 +// +// simultaneously applied to R, H = J_1^T ... J^T_(n-1) R +// so that H is upper Hessenberg. (12.5.2) */ +// +// for (k = M - 1; k > 0; k--) /* loop from k = M-1 to 1 */ +// { +// double c, s; +// double wk = gsl_vector_get (w, k); +// double wkm1 = gsl_vector_get (w, k - 1); +// +// gsl_linalg_givens (wkm1, wk, &c, &s); +// gsl_linalg_givens_gv (w, k - 1, k, c, s); +// apply_givens_qr (M, N, Q, R, k - 1, k, c, s); +// } +// +// w0 = gsl_vector_get (w, 0); +// +// /* Add in w v^T (Equation 12.5.3) */ +// +// for (j = 0; j < N; j++) +// { +// double r0j = gsl_matrix_get (R, 0, j); +// double vj = gsl_vector_get (v, j); +// gsl_matrix_set (R, 0, j, r0j + w0 * vj); +// } +// +// /* Apply Givens transformations R' = G_(n-1)^T ... G_1^T H +// Equation 12.5.4 */ +// +// for (k = 1; k < GSL_MIN(M,N+1); k++) +// { +// double c, s; +// double diag = gsl_matrix_get (R, k - 1, k - 1); +// double offdiag = gsl_matrix_get (R, k, k - 1); +// +// gsl_linalg_givens (diag, offdiag, &c, &s); +// apply_givens_qr (M, N, Q, R, k - 1, k, c, s); +// +// gsl_matrix_set (R, k, k - 1, 0.0); /* exact zero of G^T */ +// } +// +// return GSL_SUCCESS; +// } +// } + +int +gsl_linalg_complex_QR_QRsolve (gsl_matrix_complex * Q, gsl_matrix_complex * R, const gsl_vector_complex * b, gsl_vector_complex * x) +{ + const size_t M = R->size1; + const size_t N = R->size2; + + if (M != N) + { + return GSL_ENOTSQR; + } + else if (Q->size1 != M || b->size != M || x->size != M) + { + return GSL_EBADLEN; + } + else + { + /* compute sol = Q^* b */ + + gsl_blas_zgemv (CblasConjTrans, GSL_COMPLEX_ONE, Q, b, GSL_COMPLEX_ZERO, x); + + /* Solve R x = sol, storing x in-place */ + + gsl_blas_ztrsv (CblasUpper, CblasNoTrans, CblasNonUnit, R, x); + + return GSL_SUCCESS; + } +} diff --git a/linalg/test.c b/linalg/test.c index a5273085..e30b9f8c 100644 --- a/linalg/test.c +++ b/linalg/test.c @@ -2,6 +2,7 @@ * * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004, 2005, 2006, 2007, 2010 Gerard Jungman, Brian Gough * Copyright (C) Huan Wu (test_choleskyc_invert and test_choleskyc_invert_dim) + * Copyright (C) 2017 Christian Krueger (test_QRc_*) * * 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 @@ -158,6 +159,19 @@ create_general_matrix(unsigned long size1, unsigned long size2) return m; } +gsl_matrix_complex * +create_general_matrix_complex(unsigned long size1, unsigned long size2) +{ + unsigned long i, j; + gsl_matrix_complex * m = gsl_matrix_complex_alloc(size1, size2); + for(i=0; isize1; + + gsl_vector_complex * rhs = gsl_vector_complex_alloc(dim); + gsl_matrix_complex * qr = gsl_matrix_complex_alloc(dim,dim); + gsl_vector_complex * d = gsl_vector_complex_alloc(dim); + gsl_vector_complex * x = gsl_vector_complex_alloc(dim); + + gsl_matrix_complex_memcpy(qr,m); + for(i=0; isize1; + + gsl_vector_complex * rhs = gsl_vector_complex_alloc(dim); + gsl_matrix_complex * qr = gsl_matrix_complex_alloc(dim,dim); + gsl_matrix_complex * q = gsl_matrix_complex_alloc(dim,dim); + gsl_matrix_complex * r = gsl_matrix_complex_alloc(dim,dim); + gsl_vector_complex * d = gsl_vector_complex_alloc(dim); + gsl_vector_complex * x = gsl_vector_complex_alloc(dim); + + gsl_matrix_complex_memcpy(qr,m); + for(i=0; isize1, N = m->size2; + + gsl_vector_complex * rhs = gsl_vector_complex_alloc(M); + gsl_matrix_complex * qr = gsl_matrix_complex_alloc(M,N); + gsl_vector_complex * d = gsl_vector_complex_alloc(N); + gsl_vector_complex * x = gsl_vector_complex_alloc(N); + gsl_vector_complex * r = gsl_vector_complex_alloc(M); + gsl_vector_complex * res = gsl_vector_complex_alloc(M); + + gsl_matrix_complex_memcpy(qr,m); + for(i=0; isize1, N = m->size2; + + gsl_matrix_complex * qr = gsl_matrix_complex_alloc(M,N); + gsl_matrix_complex * a = gsl_matrix_complex_alloc(M,N); + gsl_matrix_complex * q = gsl_matrix_complex_alloc(M,M); + gsl_matrix_complex * r = gsl_matrix_complex_alloc(M,N); + gsl_vector_complex * d = gsl_vector_complex_alloc(GSL_MIN(M,N)); + + gsl_matrix_complex_memcpy(qr,m); + + s += gsl_linalg_complex_QR_decomp(qr, d); + s += gsl_linalg_complex_QR_unpack(qr, d, q, r); + + /* compute a = q r */ + gsl_blas_zgemm (CblasNoTrans, CblasNoTrans, GSL_COMPLEX_ONE, q, r, GSL_COMPLEX_ZERO, a); + + for(i=0; i