LAPACK sgelss使用

Question

我想使用LAPACK框架中的sgelss函数来求解过确定的线性方程组。

minimize || A*x-b||

为此，我选择使用LAPACK中的sgelss函数，头文件可以在这里找到：SGELSS。

我的问题是我接缝错误地使用它，因为结果不符合。

举个例子，我用这些材料：

double a[] =
    {
       25,5,1,
       16,4,1,
        9,3,1,
        4,2,1,
        1,1,1,
    };
double b[] =
     {
       31,21,13,7,3
     }

然后我会打电话：

sgelss_(&(5), &(3), &(1), a, (5), b, &(5), 3, -1, &(3), malloc(3000), &(3000), &output);

[由于sgelss期望大多数输入变量不是普通变量，而是我使用的指针&(.)这里，为了跳过实例化变量，然后引用它们的内存位置，以便这里的可读性！]

我希望之后的b是1,1,1，XXXXX，因为我的输入a和b是这样设置的。不幸的是，情况并非如此。

我也尝试过旋转a(切换行和列)，但没有成功。

Answer 1

在您发布的代码中有很多问题；我假设它是从您正在做的事情中抽象出来的，因为它显然无法按原样编译。有两个主要问题，我怀疑其中一个或两个都是造成实际代码中的问题的原因：

首先，您尝试将sgelss_与double数据结合使用。sgelss_在float (单精度)上运行.如果您有双精度数据，则需要使用dgelss_ .

其次，您的参数描述了一个5x3矩阵，但请记住，LAPACK使用了矩阵元素的列主要排序。这意味着代码所描述的矩阵是：

25  1  2 
 5  9  1
 1  3  1
16  1  1
 4  4  1

不知怎的，我怀疑这是你真正想要的矩阵。你更有可能想要矩阵：

25  5  1
16  4  1
 9  3  1
 4  2  1
 1  1  1

我继续编写了一个工作版本，假设这是您实际上试图使用的矩阵：

#include <Accelerate/Accelerate.h>
#include <stdio.h>
#include <stdlib.h>

int main(int argc, char *argv[]) {
    double a[] = /* column major storage order! */
    {
        25, 16,  9,  4,  1,
         5,  4,  3,  2,  1,
         1,  1,  1,  1,  1,
    };
    double b[] =
    {
        31,21,13,7,3
    };

    // Setup parameters
    __CLPK_integer m = 5;
    __CLPK_integer n = 3;
    __CLPK_integer nrhs = 1;
    __CLPK_integer lda = 5;
    __CLPK_integer ldb = 5;
    double *s = malloc(3 * sizeof*s);
    double rcond = -1.0f; // use machine precision
    __CLPK_integer rank;
    __CLPK_integer info;
    // Query correct worksize
    double worksize;
    __CLPK_integer lwork = -1;
    dgelss_(&m, &n, &nrhs, a, &lda, b, &ldb, s, &rcond, &rank, &worksize, &lwork, &info);
    // Allocate workspace
    lwork = worksize;
    double *work = malloc(lwork * sizeof *work);
    // Do computation
    dgelss_(&m, &n, &nrhs, a, &lda, b, &ldb, s, &rcond, &rank, work, &lwork, &info);
    // Free workspace
    free(work);
    // Print result vector
    for (int i=0; i<3; ++i)
        printf("%g\t", b[i]);
    printf("\n");
    return 0;
}