From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 2150 invoked by alias); 16 Feb 2010 19:28:13 -0000 Received: (qmail 2139 invoked by uid 22791); 16 Feb 2010 19:28:12 -0000 X-SWARE-Spam-Status: No, hits=-1.1 required=5.0 tests=AWL,BAYES_00 X-Spam-Check-By: sourceware.org Received: from mx4.oit.ohio.edu (HELO mx4.oit.ohio.edu) (132.235.250.54) by sourceware.org (qpsmtpd/0.43rc1) with ESMTP; Tue, 16 Feb 2010 19:28:07 +0000 X-IronPort-Anti-Spam-Filtered: true X-IronPort-Anti-Spam-Result: ApsEAAeBekuE6wiX/2dsb2JhbACcDrBZDoQ6iFWEWwSDFA Received: from oak3a.cats.ohiou.edu (HELO oak.cats.ohiou.edu) ([132.235.8.151]) by smtpout4.oit.ohio.edu with ESMTP; 16 Feb 2010 14:28:04 -0500 Received: from bubbles.hooperlab (crab-lab.zool.ohiou.edu [132.235.227.198]) (authenticated bits=0) by oak.cats.ohiou.edu (8.13.1/8.13.1) with ESMTP id o1GJRK6B1935378 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Tue, 16 Feb 2010 14:27:20 -0500 (EST) Message-ID: <4B7AF193.2090704@ohiou.edu> Date: Tue, 16 Feb 2010 19:28:00 -0000 From: "Kevin H. Hobbs" User-Agent: Mozilla/5.0 (X11; U; Linux x86_64; en-US; rv:1.9.1.7) Gecko/20100120 Fedora/3.0.1-1.fc11 Thunderbird/3.0.1 MIME-Version: 1.0 To: gsl-discuss@sourceware.org Subject: LinAlg SV Doc? Content-Type: multipart/signed; micalg=pgp-sha1; protocol="application/pgp-signature"; boundary="------------enigBF36B31023CC9D105634F084" 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: 2010-q1/txt/msg00020.txt.bz2 This is an OpenPGP/MIME signed message (RFC 2440 and 3156) --------------enigBF36B31023CC9D105634F084 Content-Type: multipart/mixed; boundary="------------060106070300030808050802" This is a multi-part message in MIME format. --------------060106070300030808050802 Content-Type: text/plain; charset=ISO-8859-1 Content-Transfer-Encoding: quoted-printable Content-length: 1683 I'm fairly new to using the linear algebra routines so this may be obvious to everybody else but : In gsl_linalg_SV_solve how long should b and x be? Could an example be added for this section? I believe the length of b must be M and the length of x must be N but I get really confused here, mostly around : is M the number of samples or the number of parameters (below)? I will gladly contribute the function I am working on now as an example. My problem of the moment is to find the best fitting plane for sets of 7 dimensional points, though 3-D will be enough to get me started. I have a set of 9 points from f(x,y)=3Dz, and I want to find the best fitting plane, ax + by + c =3D z, though these points. My points are : f(1, 1) =3D 2 f(1, 2) =3D 2.11 f(1, 3) =3D 3.21 f(2, 1) =3D 3.97 f(2, 2) =3D 3.76 f(2, 3) =3D 4.01 f(3, 1) =3D 3.40 f(3, 2) =3D 5.47 f(3, 3) =3D 4.35 When I put these into the equation for a plane this becomes : 1a + 1b + 1c =3D 2 1a + 2b + 1c =3D 2.11 1a + 3b + 1c =3D 3.21 2a + 1b + 1c =3D 3.97 2a + 2b + 1c =3D 3.76 2a + 3b + 1c =3D 4.01 3a + 1b + 1c =3D 3.40 3a + 2b + 1c =3D 5.47 3a + 3b + 1c =3D 4.35 Then I create the matrix A and the vectors x and b for the linear system Ax =3D b (I know x and b are getting confused): 1 1 1 2.00 1 2 1 2.11 1 3 1 3.21 2 1 1 3.97 A =3D 2 2 1 b =3D 3.76 x =3D a b c 2 3 1 4.01 3 1 1 3.40 3 2 1 5.47 3 3 1 4.35 The output I get from the attached program is : 0.983333 0.366667 0.886667 which puts a nice plane through my points. Do I have my Ms and Ns in the right places? --------------060106070300030808050802 Content-Type: text/x-chdr; name="find_interpolation_parameters.h" Content-Transfer-Encoding: quoted-printable Content-Disposition: attachment; filename="find_interpolation_parameters.h" Content-length: 448 #ifndef FIND_INTERPOLATION_PARAMETERS_H #define FIND_INTERPOLATION_PARAMETERS_H #include #include #include int find_interpolation_parameters( const size_t dim,=20 double * point_data,=20 double * value_data,=20 double * parameters ); int find_interpolation_parameters_SVD( const size_t dim,=20 const size_t num, double * point_data,=20 double * value_data,=20 double * parameters ); #endif --------------060106070300030808050802 Content-Type: text/x-csrc; name="find_interpolation_parameters_SVD.c" Content-Transfer-Encoding: quoted-printable Content-Disposition: attachment; filename="find_interpolation_parameters_SVD.c" Content-length: 1194 #include "find_interpolation_parameters.h" int find_interpolation_parameters_SVD( const size_t dim,=20 const size_t num, double * point_data,=20 double * value_data, double * parameters ) { gsl_matrix * AU =3D gsl_matrix_alloc( num, dim + 1 ); gsl_matrix * V =3D gsl_matrix_alloc( dim + 1, dim + 1 ); gsl_vector * S =3D gsl_vector_alloc( dim + 1 ); gsl_vector * work =3D gsl_vector_alloc( dim + 1 ); gsl_vector * b =3D gsl_vector_alloc( num ); gsl_vector * x =3D gsl_vector_alloc( dim + 1 ); size_t i; size_t j; for ( i =3D 0; i < num; ++i ) { for ( j =3D 0; j < dim; ++j ) { gsl_matrix_set( AU, i, j, point_data[ i * dim + j ] ); } gsl_matrix_set( AU, i, dim, 1.0 ); } for ( i =3D 0; i < num; ++i ) { gsl_vector_set( b, i, value_data[i] ); } gsl_linalg_SV_decomp( AU, V, S, work ); gsl_linalg_SV_solve( AU, V, S, b, x ); for ( i =3D 0; i < dim + 1; ++i ) { parameters[i] =3D gsl_vector_get( x, i ); } gsl_matrix_free( AU ); gsl_matrix_free( V ); gsl_vector_free( S ); gsl_vector_free( work ); gsl_vector_free( b ); gsl_vector_free( x ); return EXIT_SUCCESS; } --------------060106070300030808050802 Content-Type: text/x-csrc; name="interp_test.c" Content-Transfer-Encoding: quoted-printable Content-Disposition: attachment; filename="interp_test.c" Content-length: 600 #include "find_interpolation_parameters.h" int main( int argc, char * argv[] ) { const size_t dim =3D 2; const size_t num =3D 9; double point_data[] =3D {=20 1, 1, 1, 2, 1, 3, 2, 1, 2, 2, 2, 3, 3, 1, 3, 2, 3, 3 }; =20=20 double value_data[] =3D { 2.00, 2.11, 3.21, 3.97, 3.76, 4.01, 3.40, 5.47, 4.35 }; double parameters[3]; find_interpolation_parameters_SVD(=20 dim, num, point_data, value_data, parameters ); printf( "%f %f %f\n", parameters[0], parameters[1], parameters[2] ); return EXIT_SUCCESS; } --------------060106070300030808050802-- --------------enigBF36B31023CC9D105634F084 Content-Type: application/pgp-signature; name="signature.asc" Content-Description: OpenPGP digital signature Content-Disposition: attachment; filename="signature.asc" Content-length: 252 -----BEGIN PGP SIGNATURE----- Version: GnuPG v1.4.9 (GNU/Linux) Comment: Using GnuPG with Fedora - http://enigmail.mozdev.org/ iD8DBQFLevGYqBtEuW+gRPERAmxJAJ9Hsld1c9cDgrO6I3R6YGPlbe3iSQCff3bA zcWaKD1U8b/eRUhC6x0Hyhs= =ikNj -----END PGP SIGNATURE----- --------------enigBF36B31023CC9D105634F084--