bugGNU 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

Sat 14 Oct 2023 07:18:52 PM UTC, original submission:  

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:

> (7.00, 0.00)(8.00, 0.00)(9.00, 0.00)(0.14, 0.00)(0.86, 0.00)(1.71, 0.00)(0.57, 0.00)(0.50, 0.00)(0.00, 0.00)


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

Larry Darryl Lee Jr. <llee454>

 

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

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

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

    No changes have been made to this item

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code