GNU Scientific Library - Bugs: bug #64777, gsl_linalg_complex_LU_decomp...
You are not allowed to post comments on this tracker with your current authentication level.
bug #64777: gsl_linalg_complex_LU_decomp returns incorrect results
Submitter: | Larry Darryl Lee Jr. <llee454> | ||
Submitted: | Sat 14 Oct 2023 07:18:52 PM UTC | ||
Category: | Runtime error | Severity: | 3 - Normal |
Operating System: | MacOS | Status: | None |
Assigned to: | None | Open/Closed: | Open |
Release: | 2.7.1 |
No files currently attached
Depends on the following items: None found
Items that depend on this one: None found
Carbon-Copy List
No changes have been made to this item
Powered by Savane 3.13-02a9.
Corresponding source code
here appears to be a bug in gsl_linalg_complex_LU_decomp. It returns incorrect results even for matrices that only contain real values.
Demo:
```
#include <stdio.h>
#include <gsl/gsl_linalg.h>
#include <gsl_complex.h>
#include <gsl_complex_math.h>
int main (void)
{
gsl_matrix_complex* m = gsl_matrix_complex_alloc (3, 3);
gsl_matrix_complex_set (m, 0, 0, (gsl_complex){ .dat = {1, 0} });
gsl_matrix_complex_set (m, 0, 1, (gsl_complex){ .dat = {2, 0} });
gsl_matrix_complex_set (m, 0, 2, (gsl_complex){ .dat = {3, 0} });
gsl_matrix_complex_set (m, 1, 0, (gsl_complex){ .dat = {4, 0} });
gsl_matrix_complex_set (m, 1, 1, (gsl_complex){ .dat = {5, 0} });
gsl_matrix_complex_set (m, 1, 2, (gsl_complex){ .dat = {6, 0} });
gsl_matrix_complex_set (m, 2, 0, (gsl_complex){ .dat = {7, 0} });
gsl_matrix_complex_set (m, 2, 1, (gsl_complex){ .dat = {8, 0} });
gsl_matrix_complex_set (m, 2, 2, (gsl_complex){ .dat = {9, 0} });
int signum = 0;
gsl_permutation* perm = gsl_permutation_calloc (3);
gsl_linalg_complex_LU_decomp (m, perm, &signum);
for (size_t i = 0; i < 3; i ++) {
for (size_t j = 0; j < 3; j ++) {
printf ("(%0.2f, %0.2f)",
GSL_REAL (gsl_matrix_complex_get (m, i, j)),
GSL_IMAG (gsl_matrix_complex_get (m, i, j)));
}
}
return 0;
}
```
I compiled and ran this with:
```
gcc -I. -I/usr/local/Cellar/gsl/2.7.1/include/gsl -L /usr/local/Cellar/gsl/2.7.1/lib -lgsl -lgslcblas bug.c -o bug
./bug
```
The output is:
The correct LU decomposition is:
```
(%i3) get_lu_factors (lu_factor (matrix ([1, 2, 3], [4, 5, 6], [7, 8, 9])));
[ 1 0 0 ] [ 1 0 0 ] [ 1 2 3 ]
[ ] [ ] [ ]
(%o3) [[ 0 1 0 ], [ 4 1 0 ], [ 0 - 3 - 6 ]]
[ ] [ ] [ ]
[ 0 0 1 ] [ 7 2 1 ] [ 0 0 0 ]
```
Please let me know if I made some mistake or, if this is a true bug, what the timeline for fixing it is.
Thanks,
- Larry