Documentation

Documentation

[Reinhold Bader]

  Hello,

    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.

Regards


-- 
  Dr. Reinhold Bader

  Leibniz-Rechenzentrum, Abt. Hochleistungssysteme | Tel. +49 89 35831 8825
  Boltzmannstr. 1, 85748 Garching                  | Fax  +49 89 35831 9700

Documentation (2)

[Reinhold Bader]

  Hello,

    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.

Regards


-- 
  Dr. Reinhold Bader

  Leibniz-Rechenzentrum, Abt. Hochleistungssysteme | Tel. +49 89 35831 8825
  Boltzmannstr. 1, 85748 Garching                  | Fax  +49 89 35831 9700

Re: Documentation

[Brian Gough]

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

}

Documentation

[Alberto Manuel Brandao Simoes]

[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #1: Type: text/plain, Size: 617 bytes --]

	Hello!

	As many people knows, numexp is using GSL as a source for scientific
methods. NumExp, is being a front-end to GSL, and it is normal that people
will seek documentation on NumExp for GSL methods. That's true, too, that
GSL functions wrappers change some of the arguments, and so, GSL documentation
is not very accurate for some functions.

	All this, to ask for permition to retype some of GSL documentation on
numexp GSL module documentation (of course, with the source maintained.

	Alberto
-- 
 | Alberto Manuel Brandão Simões |
 | albie@alfarrabio.di.uminho.pt |
 | http://numexp.sourceforge.net |

Re: Documentation

[Brian Gough]

Alberto Manuel Brandao Simoes writes:
 > 	As many people knows, numexp is using GSL as a source for
 > scientific methods. NumExp, is being a front-end to GSL, and it is
 > normal that people will seek documentation on NumExp for GSL
 > methods. That's true, too, that GSL functions wrappers change some
 > of the arguments, and so, GSL documentation is not very accurate
 > for some functions.
 > 	All this, to ask for permition to retype some of GSL
 > documentation on numexp GSL module documentation (of course, with
 > the source maintained.

Yes, it's under the GNU Free Documentation License (FDL) so you can do
that provided you also use the GNU FDL for your documentation and
follow the conditions in the GNU FDL.

The conditions are normal things like including the original authors
in your authors list and adding a 'History' section to your manual
saying where the different parts came from.

You'll find the full list at http://www.gnu.org/copyleft/fdl.html in
section 4.

Brian