Hi,
The function gsl_linalg_complex_householder_hv(tau, v, w) [apply the
Householder transform defined by (tau, v) to the vector w] does not
check for the mildly pathological case where length of the involved
vectors is 1. In this case, a subvector()-call in the code cannot be
executed as it exceeds the index limits and the code breaks.
I encountered this problem when implementing the QR decomposition for
complex-valued matrices.
In the real-valued case, a vector of length 1 results in tau=0 when
calculating the HH transform using gsl_linalg_householder_transform()
and the corresponding function gsl_linalg_househoulder_hv() explicitly
checks for tau=0 and exits in this case.
In the complex-valued case, the HH transform of a vector of length 1
does *not* necessarily imply tau=0 which is why the case N=1 has to be
handled separately in gsl_linalg_complex_householder_hv().
I have attached a patch for linalg/householdercomplex.c which handles
the N=1 case.
To be patched against GSL 2.4; I have cloned the github repository a few
days ago. As that's my first patch, please let me know if I need to
provide anything else or if a different format would be better.
Christian
--- householdercomplex.c 2017-11-27 14:35:24.665197607 +0000
+++ householdercomplex.new.c 2017-11-27 14:35:40.105198153 +0000
@@ -1,6 +1,7 @@
/* linalg/householdercomplex.c
*
* Copyright (C) 2001, 2007 Brian Gough
+ * Copyright (C) 2017 Christian Krueger (bugfix in gsl_linalg_complex_householder_hv)
*
* 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
@@ -228,6 +229,19 @@
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 */
