I want to solve the least square problem |Ax-b|->min with LAPACK's dgels_ in C, but I'm getting a segmentation fault error (I am aware that there was a similar problem, but the code is quite different and the answers do not apply to my problem). I already located the problem, it is right when dgels_ is executed.
Code:
#include <stdio.h>
#include <lapacke.h>
#define COL 3
#define ROW 4
int main()
{
char transa ='N';
lapack_int m, n, nrhs, lda, ldb, info;
m=ROW;
n=COL;
nrhs=1;
lda=ROW;
ldb=ROW;
double work[COL];
double A [COL*ROW] =
{ 1.1, 4.2, 1.7,
2.5, 2.1, 2.8,
3.4, 4.2, 8.5,
4.4, 5.2, 7.8 };
double b[ROW] =
{ 1.5,
2.1,
3.8,
3.4 };
printf("test 1 \n");
dgels_(&transa, &m, &n, &nrhs, A, &lda, b, &ldb, work, 0, &info);
printf("test 2 \n");
return 0;
}
The first "test 1" is printed then there is the segmentation fault error.
Edit: I tried to compile it before with values, i.e. using dgels_(transa, m, n, nrhs, A, lda, b, ldb, work, 0, info). However, then I get lots of error messages:
lapack_error.c: In function ‘main’:
lapack_error.c:32:13: warning: passing argument 1 of ‘dgels_’ makes pointer from integer without a cast [-Wint-conversion]
dgels_(transa, m, n, nrhs, A, lda, b, ldb, work, 0, info);
^~~~~~
In file included from /usr/include/lapacke.h:143:0,
from lapack_error.c:2:
/usr/include/lapacke.h:14793:6: note: expected ‘char *’ but argument is of type ‘char’
void LAPACK_dgels( char* trans, lapack_int* m, lapack_int* n, lapack_int* nrhs,
^
lapack_error.c:32:21: warning: passing argument 2 of ‘dgels_’ makes pointer from integer without a cast [-Wint-conversion]
dgels_(transa, m, n, nrhs, A, lda, b, ldb, work, 0, info);
^
In file included from /usr/include/lapacke.h:143:0,
from lapack_error.c:2:
/usr/include/lapacke.h:14793:6: note: expected ‘int *’ but argument is of type ‘int’
void LAPACK_dgels( char* trans, lapack_int* m, lapack_int* n, lapack_int* nrhs,
^
lapack_error.c:32:24: warning: passing argument 3 of ‘dgels_’ makes pointer from integer without a cast [-Wint-conversion]
dgels_(transa, m, n, nrhs, A, lda, b, ldb, work, 0, info);
^
In file included from /usr/include/lapacke.h:143:0,
from lapack_error.c:2:
/usr/include/lapacke.h:14793:6: note: expected ‘int *’ but argument is of type ‘int’
void LAPACK_dgels( char* trans, lapack_int* m, lapack_int* n, lapack_int* nrhs,
^
lapack_error.c:32:27: warning: passing argument 4 of ‘dgels_’ makes pointer from integer without a cast [-Wint-conversion]
dgels_(transa, m, n, nrhs, A, lda, b, ldb, work, 0, info);
^~~~
In file included from /usr/include/lapacke.h:143:0,
from lapack_error.c:2:
/usr/include/lapacke.h:14793:6: note: expected ‘int *’ but argument is of type ‘int’
void LAPACK_dgels( char* trans, lapack_int* m, lapack_int* n, lapack_int* nrhs,
^
lapack_error.c:32:36: warning: passing argument 6 of ‘dgels_’ makes pointer from integer without a cast [-Wint-conversion]
dgels_(transa, m, n, nrhs, A, lda, b, ldb, work, 0, info);
^~~
In file included from /usr/include/lapacke.h:143:0,
from lapack_error.c:2:
/usr/include/lapacke.h:14793:6: note: expected ‘int *’ but argument is of type ‘int’
void LAPACK_dgels( char* trans, lapack_int* m, lapack_int* n, lapack_int* nrhs,
^
lapack_error.c:32:44: warning: passing argument 8 of ‘dgels_’ makes pointer from integer without a cast [-Wint-conversion]
dgels_(transa, m, n, nrhs, A, lda, b, ldb, work, 0, info);
^~~
In file included from /usr/include/lapacke.h:143:0,
from lapack_error.c:2:
/usr/include/lapacke.h:14793:6: note: expected ‘int *’ but argument is of type ‘int’
void LAPACK_dgels( char* trans, lapack_int* m, lapack_int* n, lapack_int* nrhs,
^
lapack_error.c:32:58: warning: passing argument 11 of ‘dgels_’ makes pointer from integer without a cast [-Wint-conversion]
dgels_(transa, m, n, nrhs, A, lda, b, ldb, work, 0, info);
^~~~
In file included from /usr/include/lapacke.h:143:0,
from lapack_error.c:2:
/usr/include/lapacke.h:14793:6: note: expected ‘int *’ but argument is of type ‘int’
void LAPACK_dgels( char* trans, lapack_int* m, lapack_int* n, lapack_int* nrhs,
I then tried and look up a sample program and found this one (using sgesv but I figured it might be similar with dgels). And after rewriting dgels_(transa, m, n, nrhs, A, lda, b, ldb, work, 0, info) to dgels_(&transa, &m, &n, &nrhs, A, &lda, b, &ldb, work, 0, &info) there were no more error messages, so I thought this was the right way.
0 is not a legal value for the argument lwork of dgels. Indeed, it must be a pointer to an integer in c, since all arguments are called by reference in fortran, while all arguments are called by value in c. Moreover, the value of lwork specifies the length of the array work, which must be larger than 1: LWORK >= max( 1, MN + max( MN, NRHS ) ). The only exception is lwork=-1: in this particular case, the function returns to optimal size for the array work in work[0].
For instance, the following lines can be tried:
lapack_int lwork=n;
if (m<n){lwork=m;}
lwork=-1;
double query;
dgels_(&transa, &m, &n, &nrhs, A, &lda, b, &ldb, &query, &lwork, &info);
lwork=(int)query;
printf("the optimal size is %d \n",lwork);
work=malloc(lwork*sizeof(double));
if(work==NULL){fprintf(stderr,"maloc failed\n");exit(1);}
dgels_(&transa, &m, &n, &nrhs, A, &lda, b, &ldb, work, &lwork, &info);
printf("test 2 \n");
printf("%g %g %g\n",b[0],b[1],b[2]);
There are two calls to dgels_ : the first returns the correct size of work. Then, work is allocated and dgels_ is called again to perform the least square minimization.
The result can be different from what is expected by a regular C developper: Fortran and C order of dimensions for multidimensional arrays are different.The Lapacke wrapper provide the functions LAPACKE_dgels() and LAPACKE_dgels_work() with arguments LAPACK_COL_MAJOR/LAPACK_ROW_MAJOR and transa. Always make sure that you use the correct value for these arguments by performing a test! Try different values if it fails...
Here is a sample code based on yours showing different ways to call dgels through the Lapacke interface. It can be compiled by gcc main.c -o main -llapacke -llapack -lblas -lm -Wall.
#include <stdio.h>
#include <lapacke.h>
#define COL 3
#define ROW 4
int main()
{
char transa ='N';
lapack_int m, n, nrhs, lda, ldb, info;
m=ROW;
n=COL;
nrhs=1;
lda=ROW;
ldb=ROW;
double* work;
double A [COL*ROW] =
{ 1.1, 4.2, 1.7, 2.5,
2.1, 2.8, 3.4, 4.2,
8.5, 4.4, 5.2, 7.8 };
double b[ROW] =
{ 39.3,
27.4,
29.3,
42.1 };
printf("test 1 \n");
lapack_int lwork=n;
if (m<n){lwork=m;}
lwork=-1;
double query;
dgels_(&transa, &m, &n, &nrhs, A, &lda, b, &ldb, &query, &lwork, &info);
lwork=(int)query;
printf("the optimal size is %d \n",lwork);
work=malloc(lwork*sizeof(double));
if(work==NULL){fprintf(stderr,"maloc failed\n");exit(1);}
dgels_(&transa, &m, &n, &nrhs, A, &lda, b, &ldb, work, &lwork, &info);
printf("test 2 \n");
printf("%g %g %g\n",b[0],b[1],b[2]);
free(work);
double Aa [COL*ROW] =
{ 1.1, 4.2, 1.7, 2.5,
2.1, 2.8, 3.4, 4.2,
8.5, 4.4, 5.2, 7.8 };
double bb[ROW] =
{ 39.3,
27.4,
29.3,
42.1 };
//with lapacke
info= LAPACKE_dgels(LAPACK_COL_MAJOR,transa, m, n, nrhs, Aa, lda, bb, ldb);
printf("%g %g %g\n",bb[0],bb[1],bb[2]);
//
double Aaa [COL*ROW] =
{ 1.1, 4.2, 1.7,
2.5, 2.1, 2.8,
3.4, 4.2, 8.5,
4.4, 5.2, 7.8 };
double bbb[ROW] =
{ 16.3,
17.9,
45.8,
46 };
//with lapacke
info= LAPACKE_dgels(LAPACK_COL_MAJOR,'T', n, m, nrhs, Aaa, n, bbb, ldb);
printf("%g %g %g\n",bbb[0],bbb[1],bbb[2]);
double Aaaa [COL*ROW] =
{ 1.1, 4.2, 1.7,
2.5, 2.1, 2.8,
3.4, 4.2, 8.5,
4.4, 5.2, 7.8 };
double bbbb[ROW] =
{ 16.3,
17.9,
45.8,
46 };
// it is still possible to allocate the buffer yourself once for all if LAPACKE_dgels_work is used.
// it can be useful if dgels is used many times, using the same transa, m,n,nrhs,lda,ldb but different A or b.
info= LAPACKE_dgels_work(LAPACK_COL_MAJOR,'T', n, m, nrhs, Aaaa, n, bbbb, ldb,&query,-1);
lwork=(int)query;
printf("the optimal size is %d \n",lwork);
work=malloc(lwork*sizeof(double));
if(work==NULL){fprintf(stderr,"maloc failed\n");exit(1);}
info= LAPACKE_dgels_work(LAPACK_COL_MAJOR,'T', n, m, nrhs, Aaaa, n, bbbb, ldb,work,lwork);
free(work);
printf("%g %g %g\n",bbbb[0],bbbb[1],bbbb[2]);
return 0;
}