From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 20331 invoked by alias); 26 Aug 2008 14:54:17 -0000 Received: (qmail 20314 invoked by uid 22791); 26 Aug 2008 14:54:16 -0000 X-Spam-Check-By: sourceware.org Received: from mail.network-theory.co.uk (HELO mail.network-theory.co.uk) (66.199.228.187) by sourceware.org (qpsmtpd/0.31) with ESMTP; Tue, 26 Aug 2008 14:49:38 +0000 Date: Tue, 26 Aug 2008 14:54:00 -0000 Message-ID: From: Brian Gough To: Bader@lrz.de Cc: gsl-discuss@sourceware.org Subject: Re: Documentation In-Reply-To: <48ADA750.8030500@lrz.de> <48ADC452.30707@lrz.de> References: <48ADA750.8030500@lrz.de> User-Agent: Wanderlust/2.14.0 (Africa) Emacs/22.2 Mule/5.0 (SAKAKI) MIME-Version: 1.0 (generated by SEMI 1.14.6 - "Maruoka") Content-Type: multipart/mixed; boundary="Multipart_Tue_Aug_26_15:47:27_2008-1" X-Message-Mac: 7ba17b1dc24372055ce4c93a61d346b7 Mailing-List: contact gsl-discuss-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sourceware.org X-SW-Source: 2008-q3/txt/msg00021.txt.bz2 --Multipart_Tue_Aug_26_15:47:27_2008-1 Content-Type: text/plain; charset=US-ASCII Content-length: 1199 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 --Multipart_Tue_Aug_26_15:47:27_2008-1 Content-Type: text/x-csrc; charset=US-ASCII Content-Disposition: attachment; filename="qrpt.c" Content-Transfer-Encoding: quoted-printable Content-length: 2520 #include #include #include int main () { /* a =3D 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] =3D qr(a) w =3D 0.164588 0.502226 0.655178 0.055821 v =3D 0.97662 0.27795 0.39237 0.19906 =20=20=20=20=20 */ double q_data[] =3D { -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[] =3D { -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[] =3D { 0, 0, 1, 0, /* p =3D 1,2,0,3 */ 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1 } ; double w_data[] =3D { 0.164588, 0.502226, 0.655178, 0.055821 } ; double v_data[] =3D { 0.97662,=20=20=20 0.27795,=20=20=20 0.39237,=20=20=20 0.19906 } ; gsl_matrix_view q =3D gsl_matrix_view_array (q_data, 4, 4); gsl_matrix_view r =3D gsl_matrix_view_array (r_data, 4, 4); gsl_vector_view w =3D gsl_vector_view_array (w_data, 4); gsl_vector_view v =3D gsl_vector_view_array (v_data, 4); gsl_permutation * p =3D gsl_permutation_alloc(4); p->data[0] =3D 1; p->data[1] =3D 2; p->data[2] =3D 0; p->data[3] =3D 3; gsl_matrix * z =3D gsl_matrix_alloc(4,4); gsl_linalg_QRPT_update (&q.matrix, &r.matrix, p, &w.vector, &v.vector); printf("q2=3D"); gsl_matrix_fprintf(stdout, &q.matrix, "%g"); printf("r2=3D"); gsl_matrix_fprintf(stdout, &r.matrix, "%g"); /* q*(r+w*v'*p) ans =3D 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=3D"); gsl_matrix_fprintf(stdout, z, "%g"); } --Multipart_Tue_Aug_26_15:47:27_2008-1--