3

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.

Improve this question
6
  • See the documentation at netlib.org/lapack/lapacke.html#_function_list and then scroll down to DGELS. It seems you are calling it with lots of pointers where values would suffice. Also, turn warnings of your compiler on, so it would complain if you call it with the wrong type of arguments. Commented Oct 15, 2017 at 17:31
  • 1
    As for style: it is customary to have preprocessor macros in uppercase so you won't confuse them with macros. So for #define col 3 rather use #define COL 3 Commented Oct 15, 2017 at 17:33
  • Use a debugger, and examine the values of the parameters to dgels_ after the segfault happens. Commented Oct 15, 2017 at 18:40
  • @PaulOgilvie: Yes, warnings are turned on, but gcc does not tell me any. In fact, 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 (see edited original post). Commented Oct 15, 2017 at 21:46
  • 1
    @PaulOgilvie: I tried LAPACKE_dgels(LAPACK_COL_MAJOR,transa,m,n,nrhs,A,lda,b,ldb) instead of dgels_ as it is shown in the documentation you linked. It compiles without warnings and there is no segmentation fault. However, it still doesn't work properly, which is a different question though. Commented Oct 15, 2017 at 21:57

1 Answer 1

Reset to default
1

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;
}
Improve this answer
Sign up to request clarification or add additional context in comments.

1 Comment

Thank you very much for your detailed explanation and the correct code, it solved both my original problem with dgels_ and the one I had with LAPACK_dgels.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.