* Re: Documentation
2008-08-21 19:40 ` Documentation (2) Reinhold Bader
@ 2008-08-26 14:54 ` Brian Gough
0 siblings, 0 replies; 3+ messages in thread
From: Brian Gough @ 2008-08-26 14:54 UTC (permalink / raw)
To: Bader; +Cc: gsl-discuss
[-- Attachment #1: Type: text/plain, Size: 1199 bytes --]
At Thu, 21 Aug 2008 19:35:12 +0200,
Reinhold Bader wrote:
> the routine fgsl_linalg_QR_update appears to be incorrectly
> documented. Comments in the source as well as my own testing
> says that the update formula is
>
> Q'R' = Q ( R + w v^T )
>
> and not Q'R' = Q R + w v^T, as the 1.10 manual presently says.
Thanks, you're right, I've corrected the manual to match the comment
in the source.
> the routine gsl_linalg_QRPT_update appears to be also
> incorrectly documented; in this case I think the comments in the
> source are also partially incorrect. My own testing implies
> that the update formula is
>
> Q'R' = Q ( R + w v^T ) P
>
> and not Q'R' = Q R + w v^T, as the 1.10 manual presently says,
> or
> Q'R' = Q ( R + w v^T P), as mentioned in the qrpt source
> file.
Yes, the manual is wrong. I think the comment in the source is right
though -- the test program below reproduces Q(R + w v^T P) (let me
know if I'm missing something here!).
I've added some tests for gsl_linalg_QRPT_update, since they were
missing, and also extended it to handle rectangular matrices.
Thanks for the bug report.
--
Brian Gough
[-- Attachment #2: qrpt.c --]
[-- Type: text/x-csrc, Size: 2533 bytes --]
#include <stdio.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
int main ()
{
/* a =
0.918312 0.891045 0.575062 0.548165
0.046109 0.638387 0.365323 0.680073
0.816039 0.903395 0.227423 0.499313
0.184944 0.199351 0.674672 0.625001
[q,r,p] = qr(a)
w =
0.164588
0.502226
0.655178
0.055821
v =
0.97662
0.27795
0.39237
0.19906
*/
double q_data[] = { -0.621218, 0.166224, -0.417979, -0.641678,
-0.445071, 0.045123, 0.884341, -0.133478,
-0.629828, -0.394970, -0.200523, 0.638048,
-0.138983, 0.902404, -0.054994, 0.404138 };
double r_data[] = { -1.43435, -0.75684, -1.13066, -1.04456,
0.00000, 0.63107, -0.00069, 0.48859,
0.00000, 0.00000, -0.51686, 0.23780,
0.00000, 0.00000, 0.00000, 0.12865 };
double p_data[] = { 0, 0, 1, 0, /* p = 1,2,0,3 */
1, 0, 0, 0,
0, 1, 0, 0,
0, 0, 0, 1 } ;
double w_data[] = { 0.164588,
0.502226,
0.655178,
0.055821 } ;
double v_data[] = { 0.97662,
0.27795,
0.39237,
0.19906 } ;
gsl_matrix_view q = gsl_matrix_view_array (q_data, 4, 4);
gsl_matrix_view r = gsl_matrix_view_array (r_data, 4, 4);
gsl_vector_view w = gsl_vector_view_array (w_data, 4);
gsl_vector_view v = gsl_vector_view_array (v_data, 4);
gsl_permutation * p = gsl_permutation_alloc(4);
p->data[0] = 1; p->data[1] = 2; p->data[2] = 0; p->data[3] = 3;
gsl_matrix * z = gsl_matrix_alloc(4,4);
gsl_linalg_QRPT_update (&q.matrix, &r.matrix, p, &w.vector, &v.vector);
printf("q2="); gsl_matrix_fprintf(stdout, &q.matrix, "%g");
printf("r2="); gsl_matrix_fprintf(stdout, &r.matrix, "%g");
/* q*(r+w*v'*p)
ans =
0.799757 0.446193 0.597557 0.482787
0.783298 0.569891 0.555280 0.783854
0.792831 0.071342 0.427551 0.420130
0.315217 0.838238 0.592062 0.707928
*/
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans,
1.0, &q.matrix, &r.matrix,
0.0, z);
printf("z="); gsl_matrix_fprintf(stdout, z, "%g");
}
^ permalink raw reply [flat|nested] 3+ messages in thread