在 LAPACK 的线性方程解中缺少一个元素
Lacking one element in the solution to linear equations by LAPACK
我对从线性方程组得到的解决方案感到非常困惑。我的目标是求解线性方程:通过 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)
,那么第一个元素是错误的或可能缺少的,谁能帮我一个忙?谢谢。顺便说一下,我正在使用的图书馆是犰狳。
根据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;
}
相关文章:
- lower_bound()返回最后一个元素
- 使用std::transform将一个范围的元素添加到另一个范围中
- 我想访问std::unique_ptr中的一个特定元素
- 基于范围的 for 循环:迭代使用一个元素扩展的向量
- 使用运算符 [] 引用 std::vector 上最后一个元素时出现问题<>
- 删除映射和分割错误中的一个过去结束元素
- Lower_bound不适用于具有 3 个元素的向量的最后一个元素
- 查找数组中第一个最小值和最后一个最大值元素之间的算术平均值
- 检查 2D 网格的某个元素是否与另一个元素共享对角线、水平线或垂直线
- 如何创建一个所有行大小不同的 2D 数组,并且用户将指定每行将有多少个元素?
- 仅显示链表的最后一个元素
- 为什么在 std::map 上移动无法将元素从一个映射移动到另一个映射
- 查找最小的下一个更大的元素
- push_back通过自行创建的对象获取最后一个元素的向量
- 使用 map.end() 访问 map 的最后一个元素
- 如何知道地图中的最后一个元素是否被删除?
- 选择一个元素而不是一个对象的数组的原因
- std::矢量保存 只推送最后一个元素
- 为什么这个选择排序算法仍然切换一个元素,当它已经是其他元素中最小的元素时?
- 动态分配列表 - 创建一个函数,用于删除所有包含偶数值的元素