I already found a similiar topic
"Problems with results for SGELS for least squares in C++"
for determining a "straight least squares line":
I have a couple of 3D points (x_i, y_i, z_i) and I set up the following matrices like described in the topic mentioned above:
A = (1, x_i, y_i)
b = (z_i)
What I need are the coefficients (u, v, w) for the line z = v*x + w*y + u. Afaik, using DGELS these three coefficients are stored in the
output of b[0] = u, b[1] = v, b[2] = w, right? Thus, e.g. setting x1=0, y1=1 and x2=1, y2=1 I obtain two points laying on the straight least squares line.
I was checking my 3D points with the line obtained, and it looks "good", but cross-checking it with a different approach of calculations I also obtained
another straight line that is significantly different. Here is the code I use:
- Code: Select all
#include <iostream>
#include <iomanip>
using namespace std;
extern "C" void dgels_(char *trans, int *m, int *n, int *nrhs,
double *a, int *lda, double *b, int *ldb,
double *work, int *lwork, int *info);
int main(int argv, char** argc)
{
char trans = 'N';
int m = 5;
int n = 3;
int nrhs = 1;
int lda = m; // LDA >= max(1,M)
int ldb = m; // LDB >= MAX(1,M,N)
int lwork = 8; // LWORK >= max( 1, min(M, N) + max( min(M, N), NRHS ) )
int info;
double work[8]; // max(1, LWORK)
double a[15];
double b[5];
// for dlassq
int tmp = m - n;
int incx = 1;
double scale;
double sumsq;
// first column
a[0] = 1;
a[1] = 1;
a[2] = 1;
a[3] = 1;
a[4] = 1;
// cartn_x column
a[5] = 0.535;
a[6] = 2.165;
a[7] = 2.393;
a[8] = 3.222;
a[9] = 0.961;
// cartn_y column
a[10] = -0.503;
a[11] = -0.420;
a[12] = 2.527;
a[13] = 3.152;
a[14] = 4.609;
// cartn_z column
b[0] = -11.613;
b[1] = -8.169;
b[2] = -5.757;
b[3] = -2.087;
b[4] = 0.614;
cout << "input A:" << endl;
for(int i=0; i < m; i++)
{
for(int j=0; j < n*m; j = j + m)
{
cout << setprecision(5) << setw(10) << a[j + i];
}
cout << ""<< endl;
}
cout << "------------------------------------------------------------\n"<< endl;
dgels_(&trans,&m,&n,&nrhs,a,&lda,b,&ldb,work,&lwork,&info);
if (info == 0)
{
cout << "output (after dgels): first column, second column, ..."<< endl;
for(int i=0; i < m; i++)
{
cout << "i: " << i << " " << "no col: " << n << setprecision(5) << setw(8) << b[i] << endl;
}
}
else cout << "Problem!!! info = " << info << "\n";
cout << "------------------------------------------------------------\n"<< endl;
return 0;
}
Output is:
b[0] = -9.9631
b[1] = 0.45284
b[2] = 1.9864
Is this the right approach of determining the straight line coordinates?
Thanks a lot in advance,
Gerd
