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] |
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)=z, and I want to find the best fitting plane, ax + by + c = z, though these points. My points are : f(1, 1) = 2 f(1, 2) = 2.11 f(1, 3) = 3.21 f(2, 1) = 3.97 f(2, 2) = 3.76 f(2, 3) = 4.01 f(3, 1) = 3.40 f(3, 2) = 5.47 f(3, 3) = 4.35 When I put these into the equation for a plane this becomes : 1a + 1b + 1c = 2 1a + 2b + 1c = 2.11 1a + 3b + 1c = 3.21 2a + 1b + 1c = 3.97 2a + 2b + 1c = 3.76 2a + 3b + 1c = 4.01 3a + 1b + 1c = 3.40 3a + 2b + 1c = 5.47 3a + 3b + 1c = 4.35 Then I create the matrix A and the vectors x and b for the linear system Ax = 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 = 2 2 1 b = 3.76 x = 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?
#ifndef FIND_INTERPOLATION_PARAMETERS_H #define FIND_INTERPOLATION_PARAMETERS_H #include <stdio.h> #include <stdlib.h> #include <gsl/gsl_linalg.h> int find_interpolation_parameters( const size_t dim, double * point_data, double * value_data, double * parameters ); int find_interpolation_parameters_SVD( const size_t dim, const size_t num, double * point_data, double * value_data, double * parameters ); #endif
#include "find_interpolation_parameters.h" int find_interpolation_parameters_SVD( const size_t dim, const size_t num, double * point_data, double * value_data, double * parameters ) { gsl_matrix * AU = gsl_matrix_alloc( num, dim + 1 ); gsl_matrix * V = gsl_matrix_alloc( dim + 1, dim + 1 ); gsl_vector * S = gsl_vector_alloc( dim + 1 ); gsl_vector * work = gsl_vector_alloc( dim + 1 ); gsl_vector * b = gsl_vector_alloc( num ); gsl_vector * x = gsl_vector_alloc( dim + 1 ); size_t i; size_t j; for ( i = 0; i < num; ++i ) { for ( j = 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 = 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 = 0; i < dim + 1; ++i ) { parameters[i] = 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; }
#include "find_interpolation_parameters.h" int main( int argc, char * argv[] ) { const size_t dim = 2; const size_t num = 9; double point_data[] = { 1, 1, 1, 2, 1, 3, 2, 1, 2, 2, 2, 3, 3, 1, 3, 2, 3, 3 }; double value_data[] = { 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( dim, num, point_data, value_data, parameters ); printf( "%f %f %f\n", parameters[0], parameters[1], parameters[2] ); return EXIT_SUCCESS; }
Attachment:
signature.asc
Description: OpenPGP digital signature
Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|
Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |