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

Question

我对从线性方程组得到的解决方案感到非常困惑。我的目标是求解线性方程：通过 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;
}

向量 alpha 在对角线上，

向量 β 在次对角线和超对角线上，以构造三对角矩阵，假设它是 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)，那么第一个元素是错误的或可能缺少的，谁能帮我一个忙？谢谢。顺便说一下，我正在使用的图书馆是犰狳。

Answer 1

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

在退出时，DL 被上三角矩阵 U 的第二个超对角线的 (n-2( 元素覆盖，来自 A 的 LU 分解，在 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;
}