# Passing C++ vectors to functions as arguments by reference and value

• 03-12-2014
FortranLevelC++
Passing C++ vectors to functions as arguments by reference and value
Hello. I have a program that is working very well when I pass C++ vectors as arguments to my functions by reference, but I get some compilation errors when try to make a modification. I am also posting the entire program and its output below. so that you can see what is going on. I have commented out the line that causes an error.(I apologize for some of the indentation that got corrupted when I copied the code to the browser.)

This program basically calculates the coefficients of a least square polynomial and then evaluates this polynomial at artificial data points and verifies that this actually reproduces the original data within reasonable floating point error.

The function that computes the coefficients of the least square polynomial is
Code:

` vector<double> LSPVecValued_GSL( const int,  const vector<float> &, const vector<float> &);`
and as you can see it returns a vector by value, and this vector contains the coefficients of the least square polynomial.

There is also a function that evaluates this polynomial by accepting a vector argument by reference :
Code:

`float evaluate_polynomial(double, vector<double>& ) ;`
I have also created another version of the evaluation function which accepts the same vector argument by value:
Code:

`float evaluate_polynomial_ByValue(double t, vector<double> vec_a) ;`
In the program I call the first evaluation function (whose vector argument is passed by reference) by first using an intermediate vector variable containing the coefficients, and then I pass this vector as an argument to the evaluation function, as follows:
Code:

``` vec_a = LSPVecValued_GSL(  deg,  vec_x , vec_y); for(int j=0; j< n ; j=j+20 )     {         cout<<"x["<<j<<"] = " << vec_x[j] << " ,y["<<j<<"] = " << vec_y[j] <<" , p(x["<<j<<"]) ( EVALUATED FROM REFERENCE) = "           <<  evaluate_polynomial( vec_x[j], vec_a)  << endl; // This version works without error         cout<<"x["<<j<<"] = " << vec_x[j] << " ,y["<<j<<"] = " << vec_y[j] <<" , p(x["<<j<<"]) (EVALUATED FROM VALUE) = "           <<  evaluate_polynomial_ByValue( vec_x[j], LSPVecValued_GSL(  deg,  vec_x , vec_y) ) << endl; // This also works w/o error       /*  cout<<"x["<<j<<"] = " << vec_x[j] << " ,y["<<j<<"] = " << vec_y[j] <<" , p(x["<<j<<"]) (EVALUATED FROM REFERENCE) = "           <<  evaluate_polynomial( vec_x[j],  LSPVecValued_GSL(  deg,  vec_x , vec_y))  << endl; */ // But this does not compile!!!     }```
As you can see above, I am also able to call the second evaluation function (the one whose vector argument is passed by value) directly by plugging in the function LSPVecValued_GSL"(...)" and this works without error, and this is a one step process, only one line of code is involved.

However, I get a compilation error (line number 12 that I have commented out above) if I try to plug in the function "LSPVecValued_GSL(...)" into the first evaluation function that expects a vector argument by reference. I tried to put a "&" in front ofLSPVecValued_GSL but this did not fix the bug.

What syntax is appropriate to use the first evaluation function (which accepts a vector argument by reference) if I want to plug in the vector-valued function LeastSquarePolynomial_GSL directly in the the first version of the evaluation function which expects a vector argument by reference?

MANY THANKS FOR YOUR PRECIOUS TIME!!!

HERE IS THE ENTIRE PROGRAM FOLLOWED BY ITS OUTPUT ( please see the line 96 (commented out) which causes the problem):
Code:

```/// COMPILATION SYNTAX:  g++ -I/usr/local/include -L/usr/local/lib  example4.cpp -lgsl -lgslcblas -lm /// Here the necessary libraries are libgsl, libgslcblas, and libm  // // One issue is that the double precision is very important for matrix LS computation, whereas the original data is float. // So be careful about casting the data from float to double. #include <iostream> #include <stdio.h> #include <gsl/gsl_multifit.h> #include <gsl/gsl_math.h> //#include[itex] #include<vector> using namespace std; float evaluate_polynomial(double, vector<double>& ) ; float evaluate_polynomial_ByValue(double t, vector<double> vec_a) ; void LeastSquarePolynomial_GSL( const int,  const vector<float> &, const vector<float> &, vector<double> &); vector<double> LSPVecValued_GSL( const int,  const vector<float> &, const vector<float> &); int main() {   int i, n, deg, dim_a ;   double xi, yi;   float temp;       vector<double> vec_a;   //vector<double> vec_y;   //vector<double> vec_x;   vector<float> vec_y;   vector<float> vec_x;     deg = 10 ;   n = 250 ;   for(int j=0; j<n; j++)      {     //yi = 1 + j + j*j; // y[j]=polynomial p(x)=1+x+x^2 evaluated at x[i]=i         yi=0;     //xi = (double)j;         xi= ((double)j)/n; // rescale to prevent catastrophic cancellation     //xi = ((float)j)/(5.0 +j);     //    temp = (double )xi; //    xi = 5.0*2*3.141592*j/n;     temp = (float)xi ;     vec_x.push_back(temp);     for(int k=0; k<=deg; k++)     {         yi= yi+ pow(xi,k);     //  yi = cos(xi) ;     }     temp=(float)yi;     vec_y.push_back(temp);   }   // LeastSquarePolynomial_GSL( deg,  vec_x , vec_y, vec_a );     vec_a = LSPVecValued_GSL(  deg,  vec_x , vec_y);     for(int k=0; k <= deg; k++)   {     // cout << "a["<<k<<"] = " << vec_a[k] << endl;         cout << "a["<<k<<"] = " <<  LSPVecValued_GSL(  deg,  vec_x , vec_y)[k] << endl; // k'th component of output vector of function       }   cout << "--------------------" << endl; //  vec_a = LSPVecValued_GSL(  deg,  vec_x , vec_y);     {   for(int j=0; j< n ; j=j+20 )   {         cout<<"x["<<j<<"] = " << vec_x[j] << " ,y["<<j<<"] = " << vec_y[j] <<" , p(x["<<j<<"]) ( EVALUATED FROM REFERENCE) = "           <<  evaluate_polynomial( vec_x[j], vec_a)  << endl; // This version works without error         cout<<"x["<<j<<"] = " << vec_x[j] << " ,y["<<j<<"] = " << vec_y[j] <<" , p(x["<<j<<"]) (EVALUATED FROM VALUE) = "           <<  evaluate_polynomial_ByValue( vec_x[j], LSPVecValued_GSL(  deg,  vec_x , vec_y) ) << endl; // This also works w/o error       /*  cout<<"x["<<j<<"] = " << vec_x[j] << " ,y["<<j<<"] = " << vec_y[j] <<" , p(x["<<j<<"]) (EVALUATED FROM REFERENCE) = "           <<  evaluate_polynomial( vec_x[j],  LSPVecValued_GSL(  deg,  vec_x , vec_y))  << endl; */ // But this does not compile!!!   }   }       vec_a.clear();   vec_x.clear();   vec_y.clear(); } void LeastSquarePolynomial_GSL( const int deg,  const vector<float>& vec_x , const vector<float>& vec_y, vector<double>& vec_a ) {  // this version takes float arguments for x and y values, but returns a double vec_a     //When using this float version instead of the original double precision version, numerical accuracy is lost, and therefore     // the maximum recommended polynomial degree = 10, Maximum number of data points = 500.     // (Always re-scaling 0.0 <= xi = xi/N <= 1.0)         int i, n, dim_a ;     double xi, yi, chisq;     gsl_matrix *X, *cov;     gsl_vector *y, *c;     float tempFloat;       n =vec_x.size();  // Of course vec_x.size() equals vec_y.size()     dim_a = deg+1;     double a[dim_a];         X = gsl_matrix_alloc(n,dim_a);     y = gsl_vector_alloc(n);       c = gsl_vector_alloc(dim_a);     cov = gsl_matrix_alloc(dim_a,dim_a);     for(int j=0; j<n; j++)        {         yi= (double)vec_y[j]; // cast into double the original floats..     xi= (double)vec_x[j];         gsl_vector_set(y, j, yi );             for(int k=0; k<= deg; k++)     {        gsl_matrix_set(X, j, k, pow( xi , k) );     }         }     {     gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(n,dim_a);     gsl_multifit_linear(X, y, c, cov, &chisq, work);                 gsl_multifit_linear_free(work);     }     for(int k=0; k< dim_a; k++)   {     a[k] = gsl_vector_get(c, (k) );     vec_a.push_back( a[k] );   }     gsl_matrix_free(X);   gsl_vector_free(y);   gsl_vector_free(c);   gsl_matrix_free(cov); } float evaluate_polynomial(double t, vector<double>& vec_a) // takes double arguments but returns a float value {   double  p=0.0; // accumulates in double precision   float temp;   int deg = vec_a.size() -1;   for(int k=0; k<=deg; k++)   {     p= p + vec_a[k]*pow(t,k) ;     //  p= p + k*k*pow(t,k) + k  ;     //cout << "Inside polynomial evaluation function vec_a["<<k<<"] = " << vec_a[k] << endl;   }   temp = (float)p;   return temp; } float evaluate_polynomial_ByValue(double t, vector<double> vec_a) // takes double arguments but returns a float value {   double  p=0.0; // accumulates in double precision   float temp;   int deg = vec_a.size() -1;   for(int k=0; k<=deg; k++)   {     p= p + vec_a[k]*pow(t,k) ;     //  p= p + k*k*pow(t,k) + k  ;     //cout << "Inside polynomial evaluation function vec_a["<<k<<"] = " << vec_a[k] << endl;   }   temp = (float)p;   return temp; } vector<double> LSPVecValued_GSL( const int deg,  const vector<float>& vec_x , const vector<float>& vec_y) {     vector<double> vec_a ; // note that this time vec_a is not reference but value.     // this version takes float arguments for x and y values, but returns a double vec_a     //When using this float version instead of the original double precision version, numerical accuracy is lost, and therefore     // the maximum recommended polynomial degree = 10, Maximum number of data points = 500.     // (Always re-scaling 0.0 <= xi = xi/N <= 1.0)         int i, n, dim_a ;     double xi, yi, chisq;     gsl_matrix *X, *cov;     gsl_vector *y, *c;     float tempFloat;       n =vec_x.size();  // Of course vec_x.size() equals vec_y.size()     dim_a = deg+1;     double a[dim_a];         X = gsl_matrix_alloc(n,dim_a);     y = gsl_vector_alloc(n);       c = gsl_vector_alloc(dim_a);     cov = gsl_matrix_alloc(dim_a,dim_a);     for(int j=0; j<n; j++)        {         yi= (double)vec_y[j]; // cast into double the original floats..     xi= (double)vec_x[j];         gsl_vector_set(y, j, yi );             for(int k=0; k<= deg; k++)     {        gsl_matrix_set(X, j, k, pow( xi , k) );     }         }     {     gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(n,dim_a);     gsl_multifit_linear(X, y, c, cov, &chisq, work);                 gsl_multifit_linear_free(work);     }     for(int k=0; k< dim_a; k++)   {     a[k] = gsl_vector_get(c, (k) );     vec_a.push_back( a[k] );   }     gsl_matrix_free(X);   gsl_vector_free(y);   gsl_vector_free(c);   gsl_matrix_free(cov);   return vec_a; }```
Here is the output of the program, which verifies that when the polynomial is evaluated, it does reproduce the original test data within reasonable roundoff error:
> g++ -I/usr/local/include -L/usr/local/lib LSPExampleVecValued1_GSL.cpp -lgsl -lgslcblas -lm
> ./a.out
a[0] = 1
a[1] = 0.999999
a[2] = 1.00004
a[3] = 0.999496
a[4] = 1.00292
a[5] = 0.99015
a[6] = 1.02055
a[7] = 0.972974
a[8] = 1.0219
a[9] = 0.98999
a[10] = 1.00198
--------------------
x[0] = 0 ,y[0] = 1 , p(x[0]) ( EVALUATED FROM REFERENCE) = 1
x[0] = 0 ,y[0] = 1 , p(x[0]) (EVALUATED FROM VALUE) = 1
x[20] = 0.08 ,y[20] = 1.08696 , p(x[20]) ( EVALUATED FROM REFERENCE) = 1.08696
x[20] = 0.08 ,y[20] = 1.08696 , p(x[20]) (EVALUATED FROM VALUE) = 1.08696
x[40] = 0.16 ,y[40] = 1.19048 , p(x[40]) ( EVALUATED FROM REFERENCE) = 1.19048
x[40] = 0.16 ,y[40] = 1.19048 , p(x[40]) (EVALUATED FROM VALUE) = 1.19048
x[60] = 0.24 ,y[60] = 1.31579 , p(x[60]) ( EVALUATED FROM REFERENCE) = 1.31579
x[60] = 0.24 ,y[60] = 1.31579 , p(x[60]) (EVALUATED FROM VALUE) = 1.31579
x[80] = 0.32 ,y[80] = 1.47058 , p(x[80]) ( EVALUATED FROM REFERENCE) = 1.47058
x[80] = 0.32 ,y[80] = 1.47058 , p(x[80]) (EVALUATED FROM VALUE) = 1.47058
x[100] = 0.4 ,y[100] = 1.6666 , p(x[100]) ( EVALUATED FROM REFERENCE) = 1.6666
x[100] = 0.4 ,y[100] = 1.6666 , p(x[100]) (EVALUATED FROM VALUE) = 1.6666
x[120] = 0.48 ,y[120] = 1.92248 , p(x[120]) ( EVALUATED FROM REFERENCE) = 1.92248
x[120] = 0.48 ,y[120] = 1.92248 , p(x[120]) (EVALUATED FROM VALUE) = 1.92248
x[140] = 0.56 ,y[140] = 2.26887 , p(x[140]) ( EVALUATED FROM REFERENCE) = 2.26887
x[140] = 0.56 ,y[140] = 2.26887 , p(x[140]) (EVALUATED FROM VALUE) = 2.26887
x[160] = 0.64 ,y[160] = 2.75728 , p(x[160]) ( EVALUATED FROM REFERENCE) = 2.75728
x[160] = 0.64 ,y[160] = 2.75728 , p(x[160]) (EVALUATED FROM VALUE) = 2.75728
x[180] = 0.72 ,y[180] = 3.47516 , p(x[180]) ( EVALUATED FROM REFERENCE) = 3.47516
x[180] = 0.72 ,y[180] = 3.47516 , p(x[180]) (EVALUATED FROM VALUE) = 3.47516
x[200] = 0.8 ,y[200] = 4.5705 , p(x[200]) ( EVALUATED FROM REFERENCE) = 4.5705
x[200] = 0.8 ,y[200] = 4.5705 , p(x[200]) (EVALUATED FROM VALUE) = 4.5705
x[220] = 0.88 ,y[220] = 6.29099 , p(x[220]) ( EVALUATED FROM REFERENCE) = 6.29099
x[220] = 0.88 ,y[220] = 6.29099 , p(x[220]) (EVALUATED FROM VALUE) = 6.29099
x[240] = 0.96 ,y[240] = 9.04402 , p(x[240]) ( EVALUATED FROM REFERENCE) = 9.04402
x[240] = 0.96 ,y[240] = 9.04402 , p(x[240]) (EVALUATED FROM VALUE) = 9.04402
>
• 03-12-2014
laserlight
Quote:

Originally Posted by FortranLevelC++
What syntax is appropriate to use the first evaluation function (which accepts a vector argument by reference) if I want to plug in the vector-valued function LeastSquarePolynomial_GSL directly in the the first version of the evaluation function which expects a vector argument by reference?

Well, since you passed the vector by (non-const) lvalue reference, you are presumably modifying the vector or its contents within evaluate_polynomial. Therefore, you should create the vector first, and then pass it as an argument:
Code:

```vector<double> temp = LSPVecValued_GSL(deg, vec_x, vec_y); cout << "x[" << j << "] = " << vec_x[j] << " ,y[" << j << "] = " << vec_y[j]     << " , p(x[" << j << "]) (EVALUATED FROM REFERENCE) = "     << evaluate_polynomial(vec_x[j], temp) << endl; // do something with temp // ...```
If it turns out that you are not modifying the vector or its contents from within evaluate_polynomial, then you should change the parameter to be a const reference (if you do not make a copy of the vector). You will then be able to call evaluate_polynomial using the same syntax as if you used pass by value for that parameter.
• 03-13-2014
FortranLevelC++
Quote:

Originally Posted by laserlight
Well, since you passed the vector by (non-const) lvalue reference, you are presumably modifying the vector or its contents within evaluate_polynomial. Therefore, you should create the vector first, and then pass it as an argument:
Code:

```vector<double> temp = LSPVecValued_GSL(deg, vec_x, vec_y); cout << "x[" << j << "] = " << vec_x[j] << " ,y[" << j << "] = " << vec_y[j]     << " , p(x[" << j << "]) (EVALUATED FROM REFERENCE) = "     << evaluate_polynomial(vec_x[j], temp) << endl; // do something with temp // ...```
If it turns out that you are not modifying the vector or its contents from within evaluate_polynomial, then you should change the parameter to be a const reference (if you do not make a copy of the vector). You will then be able to call evaluate_polynomial using the same syntax as if you used pass by value for that parameter.

LaserLight,

Many thanks! Actually, I was not intending to modify the arguments anyway. I have made the arguments of the function constants as you recommended, and this time it is working. However, this time I don't understand why the compiler is not giving an error message when the arguments are constants. Why would the error message go away when I make the argument constant?

After all, I still cannot understand the difference between the two step evaluation like this:

Code:

```vec_a = LSPVecValued_GSL(  deg,  vec_x, vec_y); cout<<  evaluate_polynomial(vec_x[j], vec_a)  << endl;  // which compiles and works without error```
and the one line version like this:

Code:

`cout <<  evaluate_polynomial(vec_x[j],  LSPVecValued_GSL(  deg,  vec_x , vec_y))  << endl;  // compiler gives an error message`
I imagined that the only difference should be that within the for loop calling evaluate_function many times will be computationally inefficient, but this was just a way of testing the syntax of my code at this stage. Otherwise, I understand that if the function tries to modify constant arguments then there will be an error.

Anyway, here is the program and its output with both direct and indirect function calls, and all versions are working:

Code:

```/// COMPILATION SYNTAX:  g++ -I/usr/local/include -L/usr/local/lib  example4.cpp -lgsl -lgslcblas -lm /// Here the necessary libraries are libgsl, libgslcblas, and libm  // // One issue is that the double precision is very important for matrix LS computation, whereas the original data is float. // So be careful about casting the data from float to double. #include <iostream> #include <stdio.h> #include <gsl/gsl_multifit.h> #include <gsl/gsl_math.h> //#include[itex] #include<vector> using namespace std; // float evaluate_polynomial(double, vector<double>& ) ; float evaluate_polynomial( const double, const vector<double>& ) ; float evaluate_polynomial_ByValue(double , vector<double> ) ; void LeastSquarePolynomial_GSL( const int,  const vector<float> &, const vector<float> &, vector<double> &); vector<double> LSPVecValued_GSL( const int,  const vector<float> &, const vector<float> &); int main() {   int i, n, deg, dim_a ;       double xi, yi;       float temp;       vector<double> vec_a;   //vector<double> vec_y;   //vector<double> vec_x;   vector<float> vec_y;   vector<float> vec_x;     deg = 10 ;   n = 250 ;     for(int j=0; j<n; j++)      {     //yi = 1 + j + j*j; // y[j]=polynomial p(x)=1+x+x^2 evaluated at x[i]=i         yi=0;     //xi = (double)j;         xi= ((double)j)/n; // rescale to prevent catastrophic cancellation     //xi = ((float)j)/(5.0 +j);     //    temp = (double )xi; //    xi = 5.0*2*3.141592*j/n;     temp = (float)xi ;     vec_x.push_back(temp);     for(int k=0; k<=deg; k++)     {         yi= yi+ pow(xi,k);     //  yi = cos(xi) ;     }     temp=(float)yi;     vec_y.push_back(temp);   }   // LeastSquarePolynomial_GSL( deg,  vec_x , vec_y, vec_a );     vec_a = LSPVecValued_GSL(  deg,  vec_x , vec_y);     for(int k=0; k <= deg; k++)   {     // cout << "a["<<k<<"] = " << vec_a[k] << endl;         cout << "a["<<k<<"] = " <<  LSPVecValued_GSL(  deg,  vec_x , vec_y)[k] << endl; // k'th component of output vector of function       }   cout << "--------------------" << endl;       {     for(int j=0; j< n ; j=j+20 )     {         cout<<"x["<<j<<"] = " << vec_x[j] << " ,y["<<j<<"] = " << vec_y[j] <<" , p(x["<<j<<"]) ( EVALUATED FROM REFERENCE INDIRECTLY) = "           <<  evaluate_polynomial( vec_x[j], vec_a)  << endl; // This version works without error         cout<<"x["<<j<<"] = " << vec_x[j] << " ,y["<<j<<"] = " << vec_y[j] <<" , p(x["<<j<<"]) (EVALUATED FROM VALUE) = "           <<  evaluate_polynomial_ByValue( vec_x[j], LSPVecValued_GSL(  deg,  vec_x , vec_y) ) << endl; // This also works w/o error     cout<<"x["<<j<<"] = " << vec_x[j] << " ,y["<<j<<"] = " << vec_y[j] <<" , p(x["<<j<<"]) (EVALUATED FROM REFERENCE DIRECTLY) = "           <<  evaluate_polynomial( vec_x[j],  LSPVecValued_GSL(  deg,  vec_x , vec_y))  << endl;  // But this does not compile!!!     }   }       vec_a.clear();   vec_x.clear();   vec_y.clear(); } void LeastSquarePolynomial_GSL( const int deg,  const vector<float>& vec_x , const vector<float>& vec_y, vector<double>& vec_a ) {  // this version takes float arguments for x and y values, but returns a double vec_a     //When using this float version instead of the original double precision version, numerical accuracy is lost, and therefore     // the maximum recommended polynomial degree = 10, Maximum number of data points = 500.     // (Always re-scaling 0.0 <= xi = xi/N <= 1.0)         int i, n, dim_a ;     double xi, yi, chisq;     gsl_matrix *X, *cov;     gsl_vector *y, *c;     float tempFloat;       n =vec_x.size();  // Of course vec_x.size() equals vec_y.size()     dim_a = deg+1;     double a[dim_a];         X = gsl_matrix_alloc(n,dim_a);     y = gsl_vector_alloc(n);       c = gsl_vector_alloc(dim_a);     cov = gsl_matrix_alloc(dim_a,dim_a);     for(int j=0; j<n; j++)        {         yi= (double)vec_y[j]; // cast into double the original floats..     xi= (double)vec_x[j];         gsl_vector_set(y, j, yi );             for(int k=0; k<= deg; k++)     {        gsl_matrix_set(X, j, k, pow( xi , k) );     }         }     {     gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(n,dim_a);     gsl_multifit_linear(X, y, c, cov, &chisq, work);                 gsl_multifit_linear_free(work);     }     for(int k=0; k< dim_a; k++)   {     a[k] = gsl_vector_get(c, (k) );     vec_a.push_back( a[k] );   }     gsl_matrix_free(X);   gsl_vector_free(y);   gsl_vector_free(c);   gsl_matrix_free(cov); } //float evaluate_polynomial(double t, vector<double>& vec_a) // takes double arguments but returns a float value float evaluate_polynomial( const double t, const vector<double>& vec_a) // takes double arguments but returns a float value {   double  p=0.0; // accumulates in double precision   float temp;   int deg = vec_a.size() -1;   for(int k=0; k<=deg; k++)   {     p= p + vec_a[k]*pow(t,k) ;     //  p= p + k*k*pow(t,k) + k  ;     //cout << "Inside polynomial evaluation function vec_a["<<k<<"] = " << vec_a[k] << endl;   }   temp = (float)p;   return temp; } float evaluate_polynomial_ByValue(double t, vector<double> vec_a) // takes double arguments but returns a float value {   double  p=0.0; // accumulates in double precision   float temp;   int deg = vec_a.size() -1;   for(int k=0; k<=deg; k++)   {     p= p + vec_a[k]*pow(t,k) ;     //  p= p + k*k*pow(t,k) + k  ;     //cout << "Inside polynomial evaluation function vec_a["<<k<<"] = " << vec_a[k] << endl;   }   temp = (float)p;   return temp; } vector<double> LSPVecValued_GSL( const int deg,  const vector<float>& vec_x , const vector<float>& vec_y) {     vector<double> vec_a ; // note that this time vec_a is not reference but value.     // this version takes float arguments for x and y values, but returns a double vec_a     //When using this float version instead of the original double precision version, numerical accuracy is lost, and therefore     // the maximum recommended polynomial degree = 10, Maximum number of data points = 500.     // (Always re-scaling 0.0 <= xi = xi/N <= 1.0)         int i, n, dim_a ;     double xi, yi, chisq;     gsl_matrix *X, *cov;     gsl_vector *y, *c;     float tempFloat;       n =vec_x.size();  // Of course vec_x.size() equals vec_y.size()     dim_a = deg+1;     double a[dim_a];         X = gsl_matrix_alloc(n,dim_a);     y = gsl_vector_alloc(n);       c = gsl_vector_alloc(dim_a);     cov = gsl_matrix_alloc(dim_a,dim_a);     for(int j=0; j<n; j++)        {         yi= (double)vec_y[j]; // cast into double the original floats..     xi= (double)vec_x[j];         gsl_vector_set(y, j, yi );             for(int k=0; k<= deg; k++)     {        gsl_matrix_set(X, j, k, pow( xi , k) );     }         }     {     gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(n,dim_a);     gsl_multifit_linear(X, y, c, cov, &chisq, work);                 gsl_multifit_linear_free(work);     }     for(int k=0; k< dim_a; k++)   {     a[k] = gsl_vector_get(c, (k) );     vec_a.push_back( a[k] );   }     gsl_matrix_free(X);   gsl_vector_free(y);   gsl_vector_free(c);   gsl_matrix_free(cov);   return vec_a; }```

The output of the program:
>
> g++ -I/usr/local/include -L/usr/local/lib LSPExampleVecValued1_GSL.cpp -lgsl -lgslcblas -lm
> ./a.out
a[0] = 1
a[1] = 0.999999
a[2] = 1.00004
a[3] = 0.999496
a[4] = 1.00292
a[5] = 0.99015
a[6] = 1.02055
a[7] = 0.972974
a[8] = 1.0219
a[9] = 0.98999
a[10] = 1.00198
--------------------
x[0] = 0 ,y[0] = 1 , p(x[0]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 1
x[0] = 0 ,y[0] = 1 , p(x[0]) (EVALUATED FROM VALUE) = 1
x[0] = 0 ,y[0] = 1 , p(x[0]) (EVALUATED FROM REFERENCE DIRECTLY) = 1
x[20] = 0.08 ,y[20] = 1.08696 , p(x[20]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 1.08696
x[20] = 0.08 ,y[20] = 1.08696 , p(x[20]) (EVALUATED FROM VALUE) = 1.08696
x[20] = 0.08 ,y[20] = 1.08696 , p(x[20]) (EVALUATED FROM REFERENCE DIRECTLY) = 1.08696
x[40] = 0.16 ,y[40] = 1.19048 , p(x[40]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 1.19048
x[40] = 0.16 ,y[40] = 1.19048 , p(x[40]) (EVALUATED FROM VALUE) = 1.19048
x[40] = 0.16 ,y[40] = 1.19048 , p(x[40]) (EVALUATED FROM REFERENCE DIRECTLY) = 1.19048
x[60] = 0.24 ,y[60] = 1.31579 , p(x[60]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 1.31579
x[60] = 0.24 ,y[60] = 1.31579 , p(x[60]) (EVALUATED FROM VALUE) = 1.31579
x[60] = 0.24 ,y[60] = 1.31579 , p(x[60]) (EVALUATED FROM REFERENCE DIRECTLY) = 1.31579
x[80] = 0.32 ,y[80] = 1.47058 , p(x[80]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 1.47058
x[80] = 0.32 ,y[80] = 1.47058 , p(x[80]) (EVALUATED FROM VALUE) = 1.47058
x[80] = 0.32 ,y[80] = 1.47058 , p(x[80]) (EVALUATED FROM REFERENCE DIRECTLY) = 1.47058
x[100] = 0.4 ,y[100] = 1.6666 , p(x[100]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 1.6666
x[100] = 0.4 ,y[100] = 1.6666 , p(x[100]) (EVALUATED FROM VALUE) = 1.6666
x[100] = 0.4 ,y[100] = 1.6666 , p(x[100]) (EVALUATED FROM REFERENCE DIRECTLY) = 1.6666
x[120] = 0.48 ,y[120] = 1.92248 , p(x[120]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 1.92248
x[120] = 0.48 ,y[120] = 1.92248 , p(x[120]) (EVALUATED FROM VALUE) = 1.92248
x[120] = 0.48 ,y[120] = 1.92248 , p(x[120]) (EVALUATED FROM REFERENCE DIRECTLY) = 1.92248
x[140] = 0.56 ,y[140] = 2.26887 , p(x[140]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 2.26887
x[140] = 0.56 ,y[140] = 2.26887 , p(x[140]) (EVALUATED FROM VALUE) = 2.26887
x[140] = 0.56 ,y[140] = 2.26887 , p(x[140]) (EVALUATED FROM REFERENCE DIRECTLY) = 2.26887
x[160] = 0.64 ,y[160] = 2.75728 , p(x[160]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 2.75728
x[160] = 0.64 ,y[160] = 2.75728 , p(x[160]) (EVALUATED FROM VALUE) = 2.75728
x[160] = 0.64 ,y[160] = 2.75728 , p(x[160]) (EVALUATED FROM REFERENCE DIRECTLY) = 2.75728
x[180] = 0.72 ,y[180] = 3.47516 , p(x[180]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 3.47516
x[180] = 0.72 ,y[180] = 3.47516 , p(x[180]) (EVALUATED FROM VALUE) = 3.47516
x[180] = 0.72 ,y[180] = 3.47516 , p(x[180]) (EVALUATED FROM REFERENCE DIRECTLY) = 3.47516
x[200] = 0.8 ,y[200] = 4.5705 , p(x[200]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 4.5705
x[200] = 0.8 ,y[200] = 4.5705 , p(x[200]) (EVALUATED FROM VALUE) = 4.5705
x[200] = 0.8 ,y[200] = 4.5705 , p(x[200]) (EVALUATED FROM REFERENCE DIRECTLY) = 4.5705
x[220] = 0.88 ,y[220] = 6.29099 , p(x[220]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 6.29099
x[220] = 0.88 ,y[220] = 6.29099 , p(x[220]) (EVALUATED FROM VALUE) = 6.29099
x[220] = 0.88 ,y[220] = 6.29099 , p(x[220]) (EVALUATED FROM REFERENCE DIRECTLY) = 6.29099
x[240] = 0.96 ,y[240] = 9.04402 , p(x[240]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 9.04402
x[240] = 0.96 ,y[240] = 9.04402 , p(x[240]) (EVALUATED FROM VALUE) = 9.04402
x[240] = 0.96 ,y[240] = 9.04402 , p(x[240]) (EVALUATED FROM REFERENCE DIRECTLY) = 9.04402
>
• 03-13-2014
laserlight
Quote:

Originally Posted by FortranLevelC++
However, this time I don't understand why the compiler is not giving an error message when the arguments are constants. Why would the error message go away when I make the argument constant?

A temporary, like that returned by LSPVecValued_GSL, cannot be bound to a non-const lvalue reference. However, it can be bound to a const lvalue reference, in which case its lifetime is extended to be the lifetime of the reference.

By the way, I suggest writing evaluate_polynomial as something like:
Code:

```float evaluate_polynomial(const double t, const vector<double>& vec_a) {     double p = 0.0; // accumulates in double precision     const vector<double>::size_type vec_a_size = vec_a.size();     for (vector<double>::size_type k = 0; k < vec_a_size; ++k)     {         p += vec_a[k] * pow(t, k);     }     return static_cast<float>(p); }```
• 03-13-2014
FortranLevelC++
Many thanks!
• 03-13-2014
grumpy
Quote:

Originally Posted by laserlight
By the way, I suggest writing evaluate_polynomial as something like:
Code:

```float evaluate_polynomial(const double t, const vector<double>& vec_a) {     double p = 0.0; // accumulates in double precision     const vector<double>::size_type vec_a_size = vec_a.size();     for (vector<double>::size_type k = 0; k < vec_a_size; ++k)     {         p += vec_a[k] * pow(t, k);     }     return static_cast<float>(p); }```

Personally, I'd specify the return type as double, write a version that uses an templated iterator, avoid the use of pow(), myself. Something like
Code:

```template<class Iter> double evaluate_polynomial(double t, Iter begin, Iter end) {       double result = 0.0, tpower = 1.0;       while (begin != end)       {             result += *begin * tpower;             tpower *= t;             ++begin;       }       return result; } double evaluate_polynomial(double t, const std::vector<double> & vec_a) {     return evaluate_polynomial(t, vec_a.begin(), vec_a.end()); }```
In C++11, I'd consider using a range-based for loop, as an alternative to template machinery (I'm still of two minds which is better - both have strengths and weaknesses).
• 03-13-2014
phantomotap
Quote:

Personally, I'd specify the return type as double, write a version that uses an templated iterator, avoid the use of pow(), myself.
O_o

If you are going that route, you should use templates and traits to determine the result and intermediate representation.

With C++11, the mechanic should be an accumulator.

Soma
• 03-13-2014
grumpy
Quote:

Originally Posted by phantomotap
If you are going that route, you should use templates and traits to determine the result and intermediate representation.

Yes, indeed. I wasn't going to write such a thing on the fly in a forum though.

Quote:

Originally Posted by phantomotap
With C++11, the mechanic should be an accumulator.

I would use the word "can" rather than "should". I'm not disagreeing with you, but I consider it best to make a decision on approach rather than blindly follow such a rule.
• 03-13-2014
phantomotap
O_o

Indeed. I was referring to this "flavor" of algorithm not attempting to make a rule.

Soma
• 03-13-2014
rcgldr
If interested, here is a link to a document (.rtf) of an alternate method of doing a polynomial least squares fit that uses orthogonal polynomials, eliminating the need to invert a matrix (this improves the evaluation of higher order polynomials). It includes example C code.

http://rcgldr.net/misc/opls.rtf
• 03-13-2014
rcgldr
Quote:

Originally Posted by rcgldr
If interested, here is a link to a document (.rtf) of an alternate method of doing a polynomial least squares fit that uses orthogonal polynomials, It includes example C code.

The example code is really psuedo code, using an index of [-1][] on one of the variables. Here is example code, it's old. Note that f[] is the y values, where y = f(x).

Code:

```static int m, n;                /* input values */ static double x[nmax]; static double f[nmax]; static double w[nmax]; static double b[mmax];          /* generated values */ static double A[mmax]; static double B[mmax]; static double L[mmax]; static double W[mmax]; static double p2[nmax]; static double p1[nmax]; static double p0[nmax]; static double c[mmax];          /* coefficients for y(x) */ static double z[nmax];          /* calculated f[] */ static double xi, zi, vr; static double D0, D1;          /* constants */ static double *pp2;            /* pointers to logical p2, p1, p0 */ static double *pp1; static double *pp0; static double *ppx; /*------------------------------------------------------*/ /*      polyf          poly fit                        */ /*      in:            x[], f[], w[], n, m            */ /*      out:            b[], A[], B[], L[], W[]        */ /*------------------------------------------------------*/ static void polyf() { register int i, j;     D0 = (double)0.;            /* init */     D1 = (double)1.;     pp2 = p2;     pp1 = p1;     pp0 = p0;     j = 0;     A[j] = D0;                  /* calc A, p[j], p[j-1], L, W */     L[j] = D0;                  /* note A[0] not used */     W[j] = D0;     for(i = 0; i < n; i++){         pp0[i] = D1;         pp1[i] = D0;         L[j] += w[i];         W[j] += w[i]*f[i];}     B[0] = D0;     b[j] = W[j]/L[j];     for(j = 1; j <= m; j++){         ppx = pp2;              /* save old p[j], p[j-1] */         pp2 = pp1;         pp1 = pp0;         pp0 = ppx;         A[j] = D0;              /* calc A */         for(i = 0; i < n; i++){             A[j] += w[i]*x[i]*pp1[i]*pp1[i]/L[j-1];}         L[j] = D0;              /* calc p[j], L, W */         W[j] = D0;         for(i = 0; i < n; i++){             pp0[i] = (x[i]-A[j])*pp1[i]-B[j-1]*pp2[i];             L[j] += w[i]*pp0[i]*pp0[i];             W[j] += w[i]*f[i]*pp0[i];}         B[j] = L[j]/L[j-1];    /* calc B[], b[] */         b[j] = W[j]/L[j];} } /*------------------------------------------------------*/ /*      gcoef  generate coefficients                  */ /*      in:    b[], A[], B[]                          */ /*      out:    c[]                                    */ /*      uses:  p0[], p1[], p2[]                        */ /*------------------------------------------------------*/ static void gcoef() { register int i, j;     for(i = 0; i <= m; i++){            /* init */         c[i] = p2[i] = p1[i] = p0[i] = 0.;}     p0[0] = D1;     c[0] += b[0]*p0[0];     for(j = 1; j <= m; j++){            /* generate coefs */         p2[0] = p1[0];         p1[0] = p0[0];         p0[0] = -A[j]*p1[0]-B[j-1]*p2[0];         c[0] += b[j]*p0[0];         for(i = 1; i <= j; i++){             p2[i] = p1[i];             p1[i] = p0[i];             p0[i] = p1[i-1]-A[j]*p1[i]-B[j-1]*p2[i];             c[i] += b[j]*p0[i];}} }```
• 03-13-2014
King Mir
Quote:

Originally Posted by FortranLevelC++
However, this time I don't understand why the compiler is not giving an error message when the arguments are constants. Why would the error message go away when I make the argument constant?

Laserlight answered this partially, but there is a larger reason. Logically a parameter passed by non-const reference is an output parameter, who's later state reflects the result of the function. A temporary cannot be accessed after it's passed to the function, so it cannot be a good output parameter. A temporary can be a const ref or value parameter, serving as input. In C++11 it can also be a rvalue reference, in which case you're explicitly consuming the value, and that object is not intended to be accessed after the call.

As a quirk, a temporary can be a non-const this parameter to functions and operators.
• 03-14-2014
FortranLevelC++
Quote:

Originally Posted by rcgldr
If interested, here is a link to a document (.rtf) of an alternate method of doing a polynomial least squares fit that uses orthogonal polynomials, eliminating the need to invert a matrix (this improves the evaluation of higher order polynomials). It includes example C code.

http://rcgldr.net/misc/opls.rtf

Many thanks. I will investigate this under the microscope.
• 03-14-2014
grumpy
Quote:

Originally Posted by King Mir
A temporary cannot be accessed after it's passed to the function, so it cannot be a good output parameter.

Philosophically (albeit not in C++) a temporary can be a good output parameter if the caller wishes to discard that output.
• 03-15-2014
King Mir
Quote:

Originally Posted by grumpy
Philosophically (albeit not in C++) a temporary can be a good output parameter if the caller wishes to discard that output.

Arguably, yes. The counterpoint would be if when encountering a temporary used that way, is it more likely to be a garbage value or programmer error. If the latter, then more verbose ways of discarding output are better.