使用 Lapack 的 dgesv 求解线性方程组答案

【问题标题】：Solving a system of linear equations using Lapack's dgesv使用 Lapack 的 dgesv 求解线性方程组
【发布时间】：2012-10-09 10:00:40
【问题描述】：

我想使用 C++ 中的 Lapack 包求解线性方程组。我打算像this 一样使用here 中的例程来实现它，即dgesv。

这是我的代码：

unsigned short int n=profentries.size();
double **A, *b, *x;
A= new double*[n];
b=new double[n];
x=new double[2];
// Generate the matrices for my equation
for(int i = 0; i < n; ++i)
{
    A[i] = new double[2];
}

for(int i=0; i<n; i++)
{
    A[i][0]=5;
    A[i][1]=1;
    b[i]=3;
}

std::cout<< "printing matrix A:"<<std::endl;
for(int i=0; i<n; i++)
{
    std::cout<< A[i][0]<<" " <<A[i][1]<<std::endl;
}
// Call the LAPACK solver
    x = new double[n];//probably not necessary
dgesv(A, b, 2, x); //wrong result for x!
std::cout<< "printing vector x:"<<std::endl;

/*prints 
3
3
but that's wrong, the solution is (0.6, 0)!
*/
for(int i=0; i<2; i++) 
{
    std::cout<< x[i]<<std::endl;
}

我有以下问题：

dgesv 怎么会计算包含元素 {3, 3} 的向量 x？解决方案应该是 {0.6, 0}（用 matlab 检查过）。

问候

编辑：dgesv 可能适用于方阵。下面我的解决方案展示了如何使用 dgels 解决超定系统。

【问题讨论】：

Argh，你使用的那些界面太可怕了。请帮自己一个忙，使用最新 lapack 版本附带的官方 C 绑定：netlib.org/lapack/lapacke.html
对应的函数似乎是 LAPACKE_dgels，但它既不修改通过引用传递的向量 x 也不返回向量 x，我看对了吗？那么，如何得到解x呢？
我会说 LAPACKE_dgesv 是 Fortran 例程 DGESV 的包装器，是什么让您认为它应该是 LAPACKE_dgels？根据文档，在输出时，解决方案 x 存储到 B 中。
@janneb 我采纳了所有建议，但我并没有真正得到的是如何获得包含“未知数”的解向量 x 。我的 A 有 19 行/2 列，b 有 19 行/1 列，所以 x 应该有 2 行/1 列。这些值似乎不在 b 中！
对于这些维度，您的系统是超定的，并且您不能使用 dgesv（即 dgesv 要求 A 为正方形）。如果您想要最小二乘意义上的近似解（即 IIRC，matlab '\' 运算符对非平方 A 所做的工作），请改用 dgels。

标签： c++ linear-algebra lapack equation-solving

【解决方案1】：

Lapack 假设矩阵以列主要格式存储。它不适用于矩阵 A 的指针对指针格式。矩阵的所有元素都需要连续存储在内存中。

尝试这样声明 A：

double *A;
A = new double[2*n];

还有for循环：

for(int i=0; i<n; i++)
{
    A[i+0*n]=5;
    A[i+1*n]=1;
    b[i]=3;
}

【讨论】：

好的，我不知道 Lapack 期望列主要格式。修改了您的评论后，我收到以下消息：错误：“double *”类型的参数与“double **”类型的参数不兼容。 dgesv 实际上是否需要指针到指针（即二维输入）？
我假设您使用的是常规的 lapack C 包装器。看来您正在使用的代码实际上正在为您进行转换。查看 dgesv.h 中的 dgesv_ctof() 函数
@janneb 提供的链接指出，根据参数，Lapack 可以接受行和列主要格式。编辑：这是指我的界面，而不是官方的 lapack c 绑定！
有很多不同的 Lapack 绑定。关于行优先，官方声明：“请注意，使用行优先排序可能需要比列优先排序更多的内存和时间，因为例程必须将行优先顺序转换为底层 LAPACK 所需的列优先顺序例行公事。”

【解决方案2】：

问题曾经是行主要/列主要的混淆，而且显然 dgesv 不适合非方阵。

可以使用 dgels (#include "lapacke.h") 以最小二乘方式求解超定线性方程组。

值得注意的是，对我来说乍一看并不明显：解向量（通常用 x 表示）然后存储在 b 中。

我的示例系统由一个矩阵 A 组成，矩阵 A 在其第一列中包含值 1-10，以及一个向量 b，其值为 {i|1

    double a[10][2]={1,1,2,1,3,1,4,1,5,1,6,1,7,1,8,1,9,1,10,1};
double b[10][1]={1,4,9,16,25,36,49,64,81,100};
lapack_int info,m,n,lda,ldb,nrhs;
int i,j;

// description here: http://www.netlib.org/lapack/double/dgesv.f
m = 10;
n = 2;
nrhs = 1;
lda = 2;
ldb = 1;

for(int i=0; i<m; i++)
{
    for(int k=0; k<lda; k++)
    {
        std::cout << a[i][k]<<" ";
    }
    std::cout <<"" << std::endl;
}
std::cout<< "printing vector b:"<<std::endl;
for(int i=0; i<m; i++) 
{
    std::cout<< *b[i]<<std::endl;
}
std::cout<< "\nStarting calculation..."<<std::endl;

info = LAPACKE_dgels(LAPACK_ROW_MAJOR,'N',m,n,nrhs,*a,lda,*b,ldb);

std::cout<< "Done."<<std::endl;

std::cout<< " "<<std::endl;
std::cout<< "printing vector b:"<<std::endl;
for(int i=0; i<2; i++) 
{
    std::cout<< *b[i]<<std::endl;
}
std::cout<< "Used values: "<< m << ", " <<  n << ", " << nrhs << ", " << lda << ", " << ldb << std::endl;

【讨论】：