[Bug-gsl] bug in gsl_linalg_QRPT_QRsolve

[Mario Pernici]

gsl version: 1.3
Red Hat 7.3
gcc version 2.96
  
Dear Sirs,
 in linalg/qrpt.c
  
 int
 gsl_linalg_QRPT_QRsolve (const gsl_matrix * Q,...)
 {
   // ...
   else 
     {
       /* compute b' = Q^T b */
       // gsl_blas_dgemv (CblasNoTrans, 1.0, Q, b, 0.0, x);
       // becomes
       gsl_blas_dgemv (CblasTrans, 1.0, Q, b, 0.0, x);
   // ...
  }
  ---------------------------------------------------
  Example: Solve
 {{ 2  0  1 }
  { 6  2  0 }  . x = b,     where b^T = {1, 0, 1}
  {-3 -1 -1 }}

  Solution: x^T = { 1, -3. -1}
  Below I added
   gsl_matrix_transpose(q);
  to fix the bug. With this addition the output is
  1.000000        -3.000000       -1.000000
  while without it one gets
  0.829570        -2.415429       -1.466646

Best Regards
 Mario Pernici
  ----------------------------------------------------
#include <stdio.h>
#include <gsl/gsl_linalg.h>
int main() {
  double a_data[] = {2,0,1, 6,2,0, -3,-1,-1};
  double b_data[] = {1, 0, 1};
  gsl_matrix_view m = gsl_matrix_view_array(a_data, 3, 3);
  gsl_vector_view b = gsl_vector_view_array(b_data, 3);
  gsl_matrix *q, *r;
  gsl_vector *tau, *norm, *x;
  int signum, i,j;
  gsl_permutation * p = gsl_permutation_alloc (3);

  q = gsl_matrix_alloc(3,3);
  r = gsl_matrix_alloc(3,3);
  tau = gsl_vector_alloc(3);
  norm = gsl_vector_alloc(3);
  x = gsl_vector_alloc(3);

  gsl_linalg_QRPT_decomp2(&m.matrix, q, r, tau, p, &signum, norm);              
  
  gsl_matrix_transpose(q); // added to fix the bug
  gsl_linalg_QRPT_QRsolve(q, r, p, &b.vector, x);

  for(i=0; i < 3; i++) {
    printf("%lf\t", gsl_vector_get(x, i));
  }
  puts("\n");
}