用C语言中的LAPACKE dgels_分割故障

Question

我想用LAPACK在C中的dgels_来解决最小二乘问题-Ax-b->min，但是我得到了一个分割错误错误(我知道有类似的问题，但是代码是完全不同的，答案不适用于我的问题)。我已经找到了这个问题，当dgels_被执行时，它是正确的。

代码：

#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;
}

首先打印“测试1”，然后出现分段错误。

编辑:我以前尝试过用值编译它，即使用dgels_(transa，m，n，nrhs，A，lda，b，ldb，work，0，info)。然而，我得到了大量的错误消息：

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,

然后我试着查找一个示例程序，并找到了这一个 (使用sgesv，但我认为它可能与dgels类似)。在重写dgels_(transa，m，n，nrhs，A，lda，b，ldb，work，0，info)到dgels_(&transa，&m，&n，&nrhs，A，&lda，b，&ldb，work，0，&info)之后，我认为这是正确的方法。

Answer 1

0不是lwork of 丁凝胶这一论点的法律价值。实际上，它必须是一个指向c中整数的指针，因为所有参数在fortran中都是通过引用调用的，而所有参数在c中是通过值调用的。而且，lwork的值指定数组work的长度，该数组必须大于1：LWORK >= max( 1, MN + max( MN, NRHS ) )。唯一的例外是lwork=-1：在这种特殊情况下，函数返回到work[0]中的数组work的最佳大小。

例如，可以尝试以下几行：

    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]);

有两个对dgels_的调用:第一个调用返回work的正确大小。然后，分配work，并再次调用dgels_执行最小二乘最小化。

结果可能与常规C开发人员所期望的不同：多维数组的Fortran和C维数顺序是不同的。、.The、Lapacke包装器为函数LAPACKE_dgels()和LAPACKE_dgels_work()提供了参数、LAPACK_COL_MAJOR/LAPACK_ROW_MAJOR和transa。通过执行测试，始终确保对这些参数使用正确的值！尝试不同的值如果失败..。

下面是一个基于您的示例代码，展示了通过Lapacke接口调用dgels的不同方法。它可以由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;
}