This is the mail archive of the gsl-discuss@sourceware.org mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

Re: Documentation


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



#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");

}

Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]