bugGNU Scientific Library - Bugs: bug #52540, Complex QR decomposition

 
 

You are not allowed to post comments on this tracker with your current authentication level.

bug #52540: Complex QR decomposition

Submitter:  Patrick Alken <psa>
Submitted:  Mon 27 Nov 2017 09:39:58 PM UTC
   
 
Category:  None Severity:  3 - Normal
Operating System:  Status:  Fixed
Assigned to:  None Open/Closed:  Closed
Release: 

Mon 08 Jun 2020 03:31:37 AM UTC, comment #2: 

Complex QR has been added to GSL

Patrick Alken <psa>
Group administrator
Mon 27 Nov 2017 09:40:45 PM UTC, comment #1: 

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

(file #42506)

Patrick Alken <psa>
Group administrator
Mon 27 Nov 2017 09:39:58 PM UTC, original submission:  

from chri =dot= krue =at= gmail =dot= com

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

Patrick Alken <psa>
Group administrator

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attached Files
file #42506:  patch_householdercomplex.c added by psa (948B - text/x-csrc)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by psa (Submitted the item)
  •  

    Follow 3 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2020-06-08 psa StatusNone Fixed
    2020-06-08 psa Open/ClosedOpen Closed
    2017-11-27 psa Attached File- Added patch_householdercomplex.c, #42506

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code