bug-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

## [Bug-gsl] bug in gsl_linalg_QRPT_QRsolve

 From: Mario Pernici Subject: [Bug-gsl] bug in gsl_linalg_QRPT_QRsolve Date: Sat, 26 Apr 2003 09:15:25 +0200 (CEST)

```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");
}

```

reply via email to

 [Prev in Thread] Current Thread [Next in Thread]
• [Bug-gsl] bug in gsl_linalg_QRPT_QRsolve, Mario Pernici <=