LAPACK 线性方程组的解中缺少一个元素

Posted

技术标签:

【中文标题】LAPACK 线性方程组的解中缺少一个元素【英文标题】:Lacking one element in the solution to linear equations by LAPACK 【发布时间】:2014-12-17 17:12:56 【问题描述】:

我对从线性方程组得到的解感到非常困惑。我的目标是通过 lapack 中的函数求解线性方程:A*x = e。这是我的代码:

#include <iostream>
#include "/usr/include/armadillo"
#include "/usr/local/include/lapacke.h"

using namespace std;

int main()

    int n = 3;
    arma::fvec alpha( n );//define a vetor alpha with a size 3
    arma::fvec beta( n );//define a vector beta with a size 2
    alpha << 1 << 2 << 3 << arma::endr;//assign 1,2,3 to alpha;
    beta << 0.3 << 0.6 << arma::endr;//assign 0.3, 0.6 to beta;
    float a = 0.1;
    arma::fvec e = arma::zeros<arma::fvec>( n );//define a vector e with all element equal to 0;
    e( n - 1 ) = beta( n - 2 ) * beta( n - 2 ); //the last element of e equals to the square of the last element of vector beta;
    arma::fvec tri_alpha = alpha - a;
    LAPACKE_sgtsv(LAPACK_COL_MAJOR, n, 1, &( beta[ 0 ] ), &( tri_alpha[ 0 ] ), &( beta[ 0 ] ), &( e[ 0 ] ), n );
    cout << e.t() << endl;
    return 0;

向量α在对角线上,向量β在下对角线和上对角线上构造三对角矩阵,假设为T。下面是函数sgtsv的解释。

LAPACKE_sgtsv( int matrix_order, int n, int nrhs, float *dl, float *d, float *du, float *b, int ldb)

B is REAL array, dimension (LDB,NRHS),On exit, if INFO = 0, the N by NRHS solution matrix X

在我的情况下,B = e,我最终输出 e,它是(0, -0.0444, 0.1333),显然正确答案应该是(0.0148, -0.0444, 0.1333),那么第一个元素是错误的或者可能缺少,谁能帮我一个忙?谢谢。顺便说一句,我使用的图书馆是犰狳。

【问题讨论】:

【参考方案1】:

根据sgtsv() 的documentation,在求解线性方程组时,三对角矩阵A 的上(DU)和下(DL)对角线将被覆盖。

在退出时,DL 被 A 的 LU 分解中上三角矩阵 U 的第二个上对角线的 (n-2) 个元素覆盖,在 DL(1), ..., DL(n- 2)。”

还有:

在退出时,DU 被 U 的第一个上对角线的 (n-1) 个元素覆盖。

问题出现是因为beta 应该同时是上对角线DU 和下对角线DL

解决方法是复制beta

#include <iostream>

#include "/usr/include/armadillo"
#include "/usr/local/include/lapacke.h"

//#include "/usr/local/include/armadillo"
//#include "lapacke.h"

using namespace std;

int main()

    int n = 3;
    arma::fvec alpha( n );//define a vetor alpha with a size 3
    arma::fvec beta( n );//define a vector beta with a size 2
    alpha << 1 << 2 << 3 << arma::endr;//assign 1,2,3 to alpha;
    beta << 0.3 << 0.6 << arma::endr;//assign 0.3, 0.6 to beta;
    float a = 0.1;
    arma::fvec e = arma::zeros<arma::fvec>( n );//define a vector e with all element equal to 0;
    e( n - 1 ) = beta( n - 2 ) * beta( n - 2 ); //the last element of e equals to the square of the last element of vector beta;
    arma::fvec tri_alpha = alpha - a;

    arma::fvec betabis( n );//define a vector betabis with a size 2...copy of beta
    betabis=beta;
    LAPACKE_sgtsv(LAPACK_COL_MAJOR, n, 1, &( beta[ 0 ] ), &( tri_alpha[ 0 ] ), &( betabis[ 0 ] ), &( e[ 0 ] ), n );
    cout << e.t() << endl;
    return 0;

【讨论】:

以上是关于LAPACK 线性方程组的解中缺少一个元素的主要内容,如果未能解决你的问题,请参考以下文章

《Linear Algebra and Its Application》-线性方程组的解

非齐次线性方程组Ax=b的解都有哪些?

CentOS7系统上的LAPACK源码安装

用c语言编程求线性方程组的解

线性系统笔记4线性定常系统非齐次方程的解

线性系统笔记4线性定常系统非齐次方程的解