Hi Patrick,
I don't want to insist on getting my patches introduced. But as I have,
meanwhile, learned a bit more about git and such, I thought I would make
my patch more easily accessible. I've merged my changes onto the current
git commit (as of today -- bb7f532c) and attach two clean patches (one
for the Householder fix and one for the complex QR) that could be applied.
Cheers,
Christian
On 11/27/17 10:34 PM, Patrick Alken wrote:
> Hi Christian,
>
> This should be very useful. I will take a look at your code as soon
> as I can, though it may be a while as I am very busy with other things
> at the moment. I will add this to the bug tracker so it isn't lost.
>
> Thanks,
> Patrick
>
> On 11/27/2017 08:34 AM, Christian Krueger wrote:
>> Hi,
>>
>> I was looking for the QR decomposition for complex-valued matrices. As
>> GSL does not provide this yet, I have copied the code for the
>> real-valued QR and modified it to handle complex-valued matrices.
>>
>> As I've seen in the archive that this feature has been requested a few
>> times already, I am happy to share it with everyone who wants to use it.
>> My github repository can be found at:
>> https://github.com/mangroveck/gsl-qr-complex-devel.git
>>
>> The file linalg/qrc.c is essentially a copy of the file linalg/qr.c and
>> provides a bunch of gsl_linalg_complex_QR_* functions. I've had to add
>> "_complex" to the function calls and sometimes take the conjugate of tau
>> when calling the householder-functions but the whole process was fairly
>> straight-forward and the decomposition works.
>>
>> * My repository also includes my patch for the
>> gsl_linalg_complex_householder_hv() function to handle the N=1 case.
>> (I've just reported the bug on the bug-list).
>> * I've added various tests to the file linalg/test.c to test the new
>> gsl_linalg_complex_QR_* functions (in the same style as the tests for
>> the real-valued QR decomp).
>> * I've also adjusted the header file linalg/gsl_linalg.h and the file
>> Makefile.am
>>
>> The only exception is the function gsl_linalg_complex_QR_update() which
>> I haven't adjusted yet as it requires Givens rotations and they don't
>> take complex values yet (within GSL). All other functionality from qr.c
>> is now available for complex values, too.
>>
>> If it helps, I'm happy for this code to be included in GSL.
>>
>> Christian
>>
>>
>
>
>
>
>
From 30d95fcca76e0c350b85fd5ee0f5a6f33b435bc3 Mon Sep 17 00:00:00 2001
From: Christian J Krueger <chri.k...@gmail.com>
Date: Wed, 3 Jun 2020 12:25:32 +0200
Subject: [PATCH 1/2] fix pathological case in householdercomplex.c
The pathological case of 1 dimension in gsl_linalg_complex_householder_hv() is
treated individually.
---
linalg/householdercomplex.c | 13 +++++++++++++
1 file changed, 13 insertions(+)
diff --git a/linalg/householdercomplex.c b/linalg/householdercomplex.c
index 3953c481..0cb832ba 100644
--- a/linalg/householdercomplex.c
+++ b/linalg/householdercomplex.c
@@ -123,6 +123,19 @@ gsl_linalg_complex_householder_hv (gsl_complex tau, const
gsl_vector_complex * v
if (GSL_REAL(tau) == 0.0 && GSL_IMAG(tau) == 0.0)
return GSL_SUCCESS;
+ /* treat the N=1 case separately */
+ if (N == 1)
+ {
+ gsl_complex w0 = gsl_vector_complex_get(w, 0);
+ gsl_complex tz = gsl_complex_mul(tau, w0);
+ gsl_complex ntz = gsl_complex_negative(tz);
+ gsl_complex w0ntz = gsl_complex_add(w0, ntz);
+
+ gsl_vector_complex_set(w, 0, w0ntz);
+
+ return GSL_SUCCESS;
+ }
+
{
/* compute z = v'w */
--
2.11.0
From b329da7d951d0c63080b730d859830c203cda04e Mon Sep 17 00:00:00 2001
From: Christian J Krueger <chri.k...@gmail.com>
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 <config.h>
+#include <stdlib.h>
+#include <string.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_complex_math.h>
+
+#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; i<size1; i++) {
+ for(j=0; j<size2; j++) {
+ gsl_matrix_complex_set(m, i, j,
gsl_complex_rect(1.0/(i+j+1.0),1.0/(1+size1+size2-i-j)));
+ }
+ }
+ return m;
+}
+
gsl_matrix *
create_singular_matrix(unsigned long size1, unsigned long size2)
{
@@ -190,6 +204,20 @@ create_vandermonde_matrix(unsigned long size)
return m;
}
+gsl_matrix_complex *
+create_vandermonde_matrix_complex(unsigned long size)
+{
+ unsigned long i, j;
+ gsl_matrix_complex * m = gsl_matrix_complex_alloc(size, size);
+ for(i=0; i<size; i++) {
+ for(j=0; j<size; j++) {
+ gsl_complex a = gsl_complex_rect(i + 1.0, i*0.5);
+ gsl_matrix_complex_set(m, i, j, gsl_complex_pow_real(a, size - j - 1.0));
+ }
+ }
+ return m;
+}
+
gsl_matrix *
create_moler_matrix(unsigned long size)
{
@@ -330,6 +358,10 @@ gsl_matrix * m35;
gsl_matrix * m53;
gsl_matrix * m97;
+gsl_matrix_complex * m35c;
+gsl_matrix_complex * m53c;
+gsl_matrix_complex * m97c;
+
gsl_matrix * s35;
gsl_matrix * s53;
@@ -356,6 +388,10 @@ gsl_matrix * bigsparse;
double m53_lssolution[] = {52.5992295702070, -337.7263113752073,
351.8823436427604};
+double m53c_lssolution[] = {-2.42665224693892, -6.55069697384445,
+ -5.47280378950620, 35.11460791136691,
+ 13.84982610972922,-36.12898835153614};
+
double hilb2_solution[] = {-8.0, 18.0} ;
double hilb3_solution[] = {27.0, -192.0, 210.0};
double hilb4_solution[] = {-64.0, 900.0, -2520.0, 1820.0};
@@ -384,6 +420,18 @@ double vander12_solution[] = {0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 1.0, 0.0};
+gsl_matrix_complex * vander2c;
+gsl_matrix_complex * vander3c;
+gsl_matrix_complex * vander4c;
+gsl_matrix_complex * vander12c;
+
+double vander2c_solution[] = {1.0, 0.0, 0.0, 0.0};
+double vander3c_solution[] = {0.0, 0.0, 1.0, 0.0, 0.0, 0.0};
+double vander4c_solution[] = {0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0};
+double vander12c_solution[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0};
+
gsl_matrix * moler10;
#include "test_common.c"
@@ -643,6 +691,296 @@ int test_QR_lssolve(void)
}
int
+test_QRc_solve_dim(const gsl_matrix_complex * m, const double * actual, double
eps)
+{
+ int s = 0;
+ unsigned long i, dim = m->size1;
+
+ 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; i<dim; i++) {
+ gsl_complex a = gsl_complex_rect(i+1.0, 0.5*i);
+ gsl_vector_complex_set(rhs, i, a);
+ }
+ s += gsl_linalg_complex_QR_decomp(qr, d);
+ s += gsl_linalg_complex_QR_solve(qr, d, rhs, x);
+ for(i=0; i<dim; i++) {
+ int foo_r = check(GSL_REAL(gsl_vector_complex_get(x, i)), actual[2*i],
eps);
+ int foo_i = check(GSL_IMAG(gsl_vector_complex_get(x, i)), actual[2*i+1],
eps);
+ if(foo_r || foo_i) {
+ printf("%3lu[%lu]: %22.18g %22.18g\n", dim, i,
GSL_REAL(gsl_vector_complex_get(x, i)), actual[2*i]);
+ printf("%3lu[%lu]: %22.18g %22.18g\n", dim, i,
GSL_IMAG(gsl_vector_complex_get(x, i)), actual[2*i+1]);
+ }
+ s += foo_r + foo_i;
+ }
+
+ gsl_vector_complex_free(x);
+ gsl_vector_complex_free(d);
+ gsl_matrix_complex_free(qr);
+ gsl_vector_complex_free(rhs);
+
+ return s;
+}
+
+int test_QRc_solve(void)
+{
+ int f;
+ int s = 0;
+
+ f = test_QRc_solve_dim(vander2c, vander2c_solution, 16.0 * GSL_DBL_EPSILON);
+ gsl_test(f, " QRc_solve vander(2)");
+ s += f;
+
+ f = test_QRc_solve_dim(vander3c, vander3c_solution, 64.0 * GSL_DBL_EPSILON);
+ gsl_test(f, " QRc_solve vander(3)");
+ s += f;
+
+ f = test_QRc_solve_dim(vander4c, vander4c_solution, 1024.0 *
GSL_DBL_EPSILON);
+ gsl_test(f, " QRc_solve vander(4)");
+ s += f;
+
+ f = test_QRc_solve_dim(vander12c, vander12c_solution, 0.05);
+ gsl_test(f, " QRc_solve vander(12)");
+ s += f;
+
+ return s;
+}
+
+int
+test_QRc_QRsolve_dim(const gsl_matrix_complex * m, const double * actual,
double eps)
+{
+ int s = 0;
+ unsigned long i, dim = m->size1;
+
+ 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; i<dim; i++) {
+ gsl_complex a = gsl_complex_rect(i+1.0, 0.5*i);
+ gsl_vector_complex_set(rhs, i, a);
+ }
+ s += gsl_linalg_complex_QR_decomp(qr, d);
+ s += gsl_linalg_complex_QR_unpack(qr, d, q, r);
+ s += gsl_linalg_complex_QR_QRsolve(q, r, rhs, x);
+ for(i=0; i<dim; i++) {
+ int foo_r = check(GSL_REAL(gsl_vector_complex_get(x, i)), actual[2*i],
eps);
+ int foo_i = check(GSL_IMAG(gsl_vector_complex_get(x, i)), actual[2*i+1],
eps);
+ if(foo_r || foo_i) {
+ printf("%3lu[%lu]: %22.18g %22.18g\n", dim, i,
GSL_REAL(gsl_vector_complex_get(x, i)), actual[2*i]);
+ printf("%3lu[%lu]: %22.18g %22.18g\n", dim, i,
GSL_IMAG(gsl_vector_complex_get(x, i)), actual[2*i+1]);
+ }
+ s += foo_r || foo_i;
+ }
+
+ gsl_vector_complex_free(x);
+ gsl_vector_complex_free(d);
+ gsl_matrix_complex_free(qr);
+ gsl_matrix_complex_free(q);
+ gsl_matrix_complex_free(r);
+ gsl_vector_complex_free(rhs);
+ return s;
+}
+
+int test_QRc_QRsolve(void)
+{
+ int f;
+ int s = 0;
+
+ f = test_QRc_QRsolve_dim(vander2c, vander2c_solution, 8.0 * GSL_DBL_EPSILON);
+ gsl_test(f, " QRc_QRsolve vander(2)");
+ s += f;
+
+ f = test_QRc_QRsolve_dim(vander3c, vander3c_solution, 64.0 *
GSL_DBL_EPSILON);
+ gsl_test(f, " QRc_QRsolve vander(3)");
+ s += f;
+
+ f = test_QRc_QRsolve_dim(vander4c, vander4c_solution, 1024.0 *
GSL_DBL_EPSILON);
+ gsl_test(f, " QRc_QRsolve vander(4)");
+ s += f;
+
+ f = test_QRc_QRsolve_dim(vander12c, vander12c_solution, 0.05);
+ gsl_test(f, " QRc_QRsolve vander(12)");
+ s += f;
+
+ return s;
+}
+
+int
+test_QRc_lssolve_dim(const gsl_matrix_complex * m, const double * actual,
double eps)
+{
+ int s = 0;
+ unsigned long i, M = m->size1, 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; i<M; i++) {
+ gsl_complex a = gsl_complex_rect(i+1.0, 0.5*i);
+ gsl_vector_complex_set(rhs, i, a);
+ }
+ s += gsl_linalg_complex_QR_decomp(qr, d);
+ s += gsl_linalg_complex_QR_lssolve(qr, d, rhs, x, res);
+
+ for(i=0; i<N; i++) {
+ int foo_r = check(GSL_REAL(gsl_vector_complex_get(x, i)), actual[2*i],
eps);
+ int foo_i = check(GSL_IMAG(gsl_vector_complex_get(x, i)), actual[2*i+1],
eps);
+ if(foo_r || foo_i) {
+ printf("(%3lu,%3lu)[%lu]: %22.18g %22.18g\n", M, N, i,
GSL_REAL(gsl_vector_complex_get(x, i)), actual[2*i]);
+ printf("(%3lu,%3lu)[%lu]: %22.18g %22.18g\n", M, N, i,
GSL_IMAG(gsl_vector_complex_get(x, i)), actual[2*i+1]);
+ }
+ s += foo_r + foo_i;
+ }
+
+ /* compute residual r = b - m x */
+ if (M == N) {
+ gsl_vector_complex_set_zero(r);
+ } else {
+ gsl_vector_complex_memcpy(r, rhs);
+ gsl_blas_zgemv(CblasNoTrans, gsl_complex_rect(-1.0, 0.0), m, x,
GSL_COMPLEX_ONE, r);
+ };
+
+ for(i=0; i<N; i++) {
+ int foo_r = check(GSL_REAL(gsl_vector_complex_get(res, i)),
GSL_REAL(gsl_vector_complex_get(r,i)), sqrt(eps));
+ int foo_i = check(GSL_IMAG(gsl_vector_complex_get(res, i)),
GSL_IMAG(gsl_vector_complex_get(r,i)), sqrt(eps));
+ if(foo_r || foo_i) {
+ printf("(%3lu,%3lu)[%lu]: %22.18g %22.18g\n", M, N, i,
GSL_REAL(gsl_vector_complex_get(res, i)),
GSL_REAL(gsl_vector_complex_get(r,i)));
+ printf("(%3lu,%3lu)[%lu]: %22.18g %22.18g\n", M, N, i,
GSL_IMAG(gsl_vector_complex_get(res, i)),
GSL_IMAG(gsl_vector_complex_get(r,i)));
+ }
+ s += foo_r + foo_i;
+ }
+
+ gsl_vector_complex_free(r);
+ gsl_vector_complex_free(res);
+ gsl_vector_complex_free(x);
+ gsl_vector_complex_free(d);
+ gsl_matrix_complex_free(qr);
+ gsl_vector_complex_free(rhs);
+
+ return s;
+}
+
+int test_QRc_lssolve(void)
+{
+ int f;
+ int s = 0;
+
+ f = test_QRc_lssolve_dim(m53c, m53c_lssolution, 2 * 64.0 * GSL_DBL_EPSILON);
+ gsl_test(f, " QRc_lssolve m(5,3)");
+ s += f;
+
+ f = test_QRc_lssolve_dim(vander2c, vander2c_solution, 16.0 *
GSL_DBL_EPSILON);
+ gsl_test(f, " QRc_lssolve vander(2)");
+ s += f;
+
+ f = test_QRc_lssolve_dim(vander3c, vander3c_solution, 64.0 *
GSL_DBL_EPSILON);
+ gsl_test(f, " QRc_lssolve vander(3)");
+ s += f;
+
+ f = test_QRc_lssolve_dim(vander4c, vander4c_solution, 1024.0 *
GSL_DBL_EPSILON);
+ gsl_test(f, " QRc_lssolve vander(4)");
+ s += f;
+
+ f = test_QRc_lssolve_dim(vander12c, vander12c_solution, 0.05);
+ gsl_test(f, " QRc_lssolve vander(12)");
+ s += f;
+
+ return s;
+}
+
+int
+test_QRc_decomp_dim(const gsl_matrix_complex * m, double eps)
+{
+ int s = 0;
+ unsigned long i,j, M = m->size1, 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<M; i++) {
+ for(j=0; j<N; j++) {
+ gsl_complex aij = gsl_matrix_complex_get(a, i, j);
+ gsl_complex mij = gsl_matrix_complex_get(m, i, j);
+ int foo_r = check(GSL_REAL(aij), GSL_REAL(mij), eps);
+ int foo_i = check(GSL_IMAG(aij), GSL_IMAG(mij), eps);
+ if(foo_r || foo_i) {
+ printf("(%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n", M, N, i,j,
GSL_REAL(aij), GSL_REAL(mij));
+ printf("(%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n", M, N, i,j,
GSL_IMAG(aij), GSL_IMAG(mij));
+ }
+ s += foo_r + foo_i;
+ }
+ }
+
+ gsl_vector_complex_free(d);
+ gsl_matrix_complex_free(qr);
+ gsl_matrix_complex_free(a);
+ gsl_matrix_complex_free(q);
+ gsl_matrix_complex_free(r);
+
+ return s;
+}
+
+int test_QRc_decomp(void)
+{
+ int f;
+ int s = 0;
+
+ f = test_QRc_decomp_dim(m35c, 2 * 8.0 * GSL_DBL_EPSILON);
+ gsl_test(f, " QR_decomp m(3,5)");
+ s += f;
+
+ f = test_QRc_decomp_dim(m53c, 2 * 64.0 * GSL_DBL_EPSILON);
+ gsl_test(f, " QR_decomp m(5,3)");
+ s += f;
+
+ f = test_QRc_decomp_dim(m97c, 2 * 64.0 * GSL_DBL_EPSILON);
+ gsl_test(f, " QR_decomp m(9,7)");
+ s += f;
+
+ f = test_QRc_decomp_dim(vander2c, 8.0 * GSL_DBL_EPSILON);
+ gsl_test(f, " QRc_decomp vander(2)");
+ s += f;
+
+ f = test_QRc_decomp_dim(vander3c, 64.0 * GSL_DBL_EPSILON);
+ gsl_test(f, " QRc_decomp vander(3)");
+ s += f;
+
+ f = test_QRc_decomp_dim(vander4c, 1024.0 * GSL_DBL_EPSILON);
+ gsl_test(f, " QRc_decomp vander(4)");
+ s += f;
+
+ f = test_QRc_decomp_dim(vander12c, 0.0005); /* FIXME: bad accuracy */
+ gsl_test(f, " QRc_decomp vander(12)");
+ s += f;
+
+ return s;
+}
+
+int
test_QRPT_lssolve_dim(const gsl_matrix * m, const double * actual, double eps)
{
int s = 0;
@@ -3196,6 +3534,10 @@ main(void)
m53 = create_general_matrix(5,3);
m97 = create_general_matrix(9,7);
+ m35c = create_general_matrix_complex(3,5);
+ m53c = create_general_matrix_complex(5,3);
+ m97c = create_general_matrix_complex(9,7);
+
s35 = create_singular_matrix(3,5);
s53 = create_singular_matrix(5,3);
@@ -3209,6 +3551,11 @@ main(void)
vander4 = create_vandermonde_matrix(4);
vander12 = create_vandermonde_matrix(12);
+ vander2c = create_vandermonde_matrix_complex(2);
+ vander3c = create_vandermonde_matrix_complex(3);
+ vander4c = create_vandermonde_matrix_complex(4);
+ vander12c = create_vandermonde_matrix_complex(12);
+
moler10 = create_moler_matrix(10);
c7 = create_complex_matrix(7);
@@ -3256,7 +3603,9 @@ main(void)
gsl_test(test_LUc_solve(r), "Complex LU Solve");
gsl_test(test_LUc_invert(r), "Complex LU Inverse");
gsl_test(test_QR_decomp(), "QR Decomposition");
+ gsl_test(test_QRc_decomp(), "Complex QR Decomposition");
gsl_test(test_QR_solve(), "QR Solve");
+ gsl_test(test_QRc_solve(), "Complex QR Solve");
gsl_test(test_LQ_solve(), "LQ Solve");
gsl_test(test_PTLQ_solve(), "PTLQ Solve");
@@ -3278,7 +3627,9 @@ main(void)
gsl_test(test_PTLQ_solve(), "PTLQ Solve");
gsl_test(test_QR_QRsolve(), "QR QR Solve");
+ gsl_test(test_QRc_QRsolve(), "Complex QR QR Solve");
gsl_test(test_QR_lssolve(), "QR LS Solve");
+ gsl_test(test_QRc_lssolve(), "Complex QR LS Solve");
gsl_test(test_QR_update(), "QR Rank-1 Update");
gsl_test(test_QRPT_decomp(), "QRPT Decomposition");
gsl_test(test_QRPT_lssolve(), "QRPT LS Solve");
--
2.11.0
