我对从线性方程组中获得的解决方案感到非常困惑。我的目标是解决线性方程:通过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;
}
vector alpha在对角线上, vector beta在次对角线和超对角线上,以构造三对角矩阵(假定为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)
,那么第一个元素是错误的或可能缺少,有人可以帮我吗?谢谢。顺便说一句,我正在使用的图书馆是 Armadillo 。 最佳答案
根据sgtsv()
的documentation,当求解线性方程时,三对角矩阵DU
的上对角线(DL
)和下对角线(A
)将被覆盖。
和:
之所以出现此问题,是因为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;
}