LAPACK dgetrs vs dgesv
LAPACK dgetrs vs dgesv
在 LAPACK 文档中,
dgesv 求解线性方程组 AX=B。
dgetrs 使用 DGETRF 计算的 LU 分解求解线性方程组 AX=B、A T X=B或 A H X=B。
有什么区别?另外,我制作了以下C/C++程序,两者都给出了不同的结果。为什么会这样?
#include <iostream>
#include <iomanip>
#include <vector>
#include <math.h>
#include <time.h>
#include "stdafx.h"
using namespace std;
extern "C" {
// LU decomoposition of a general matrix
void dgetrf_(int* M, int *N, double* A, int* lda, int* IPIV, int* INFO);
//// generate inverse of a matrix given its LU decomposition
//void dgetri_(int* N, double* A, int* lda, int* IPIV, double* WORK, int* lwork, int* INFO);
void dgetrs_(char* C, int* N, int* NRHS, double* A, int* LDA, int* IPIV, double* B, int* LDB, int* INFO);
void dgesv_(int *n, int *nrhs, double *a, int *lda, int *ipiv, double *b, int *ldb, int *info);
}
void solvelineq(double* A, double* B, int N)
{
int *IPIV = new int[N + 1];
int LWORK = N*N;
double *WORK = new double[LWORK];
int INFO;
char C = 'N';
int NRHS = 1;
dgetrf_(&N, &N, A, &N, IPIV, &INFO);
dgetrs_(&C, &N, &NRHS, A, &N, IPIV, B, &N, &INFO);
//dgetri_(&N, A, &N, IPIV, WORK, &LWORK, &INFO);
delete IPIV;
delete WORK;
}
double comparematrices(double* A, double* B, int N)
{
double sum = 0.;
for (int i = 0; i < N; ++i)
sum += fabs(A[i] - B[i]);
return sum;
}
int main() {
int dim;
std::cout << "Enter Dimensions of Matrix : n";
std::cin >> dim;
clock_t c1, c2;
c1 = clock();
vector<double> a(dim*dim);
vector<double> b(dim);
vector<double> c(dim);
vector<int> ipiv(dim);
srand(1); // seed the random # generator with a known value
double maxr = (double)RAND_MAX;
for (int r = 0; r < dim; r++) { // set a to a random matrix, i to the identity
for (int c = 0; c < dim; c++) {
a[r + c*dim] = rand() / maxr;
}
b[r] = rand() / maxr;
c[r] = b[r];
}
c2 = clock();
cout << endl << dim << " x " << dim << " matrix generated in : " << double(c2 - c1) / CLK_TCK << " seconds " << endl;
// C is identical matrix to B because we are solving by two methods.
c1 = clock();
solvelineq(&*a.begin(), &*c.begin(), dim);
c2 = clock();
cout << endl << " dgetrf_ and dgetrs_ completed in : " << double(c2 - c1) / CLK_TCK << " seconds " << endl;
vector<double> a1(a);
vector<double> b1(b);
int info;
int one = 1;
c1 = clock();
dgesv_(&dim, &one, &*a.begin(), &dim, &*ipiv.begin(), &*b.begin(), &dim, &info);
c2 = clock();
cout << endl << " dgesv_ completed in : " << double(c2 - c1) / CLK_TCK << " seconds " << endl;
cout << "info is " << info << endl;
double eps = 0.;
c1 = clock();
for (int i = 0; i < dim; ++i)
{
double sum = 0.;
for (int j = 0; j < dim; ++j)
sum += a1[i + j*dim] * b[j];
eps += fabs(b1[i] - sum);
}
c2 = clock();
cout << endl << " dgesv_ check completed in : " << double(c2 - c1) / CLK_TCK << " seconds " << endl;
cout << "check is " << eps << endl;
cout << "nFinal Matrix 1 : " << endl;
for (int i = 0; i < dim; ++i)
cout << b[i] << endl;
cout << "nFinal Matrix 2 : " << endl;
for (int i = 0; i < dim; ++i)
cout << c[i] << endl;
double diff;
c1 = clock();
diff = comparematrices(&*b.begin(), &*c.begin(), dim);
c2 = clock();
cout << endl << " Both functions compared in : " << double(c2 - c1) / CLK_TCK << " seconds " << endl;
cout << "Sum of difference is : " << diff << endl;
return 0;
}
我的结果 :输入矩阵的尺寸:5
5 x 5 矩阵生成时间:0 秒
dgetrf_和dgetrs_在 0.001 秒内完成
dgesv_完成时间:0 秒信息 0
dgesv_ 检查完成时间:0 秒检查是 2.02009e-15
最终矩阵 1 :-2.976294.839482.00769-1.98887-5.15561
最终矩阵 2 :-1.406222.290290.480829-1.635970.71962
两个函数的比较时间: 0 秒
差额总和为 : 11.8743
调用
dgesv
相当于调用dgetrf
和dgetrs
。dgesv
的源代码非常简单。它主要包含对dgetrf
和dgetrs
函数的调用。
示例中的结果不同,因为dgetrf
更改了A
矩阵。在 lapack 的参考文献中,它指出:
A is DOUBLE PRECISION array, dimension (LDA,N)
On entry, the M-by-N matrix to be factored.
On exit, the factors L and U from the factorization
A = P*L*U; the unit diagonal elements of L are not stored.
在您的示例中,第二次使用不同的A
矩阵。既然你打算使用lapack
,我建议你也尝试使用blas
例程(daxpy
,dnrm2
)。
我创建了一个比较dgesv
、dgetrf
和dgetrs
的简单示例
#include <iostream>
#include <vector>
#include <stdlib.h>
using namespace std;
extern "C" {
void daxpy_(int* n,double* alpha,double* dx,int* incx,double* dy,int* incy);
double dnrm2_(int* n,double* x, int* incx);
void dgetrf_(int* M, int *N, double* A, int* lda, int* IPIV, int* INFO);
void dgetrs_(char* C, int* N, int* NRHS, double* A, int* LDA, int* IPIV, double* B, int* LDB, int* INFO);
void dgesv_(int *n, int *nrhs, double *a, int *lda, int *ipiv, double *b, int *ldb, int *info);
}
void print(const char* head, int m, int n, double* a){
cout << endl << head << endl << "---------" << endl;
for(int i=0; i<m; ++i){
for(int j=0; j<n; ++j){
cout << a[i+j*m]<< ' ';
}
cout << endl;
}
cout << "---------" << endl;
}
int main() {
int dim;
std::cout << "Enter Dimensions of Matrix : n";
std::cin >> dim;
vector<double> a(dim*dim);
vector<double> b(dim);
srand(1); // seed the random # generator with a known value
double maxr = (double)RAND_MAX;
for (int r = 0; r < dim; r++) { // set a to a random matrix, i to the identity
for (int i = 0; i < dim; i++) {
a[r + i*dim] = rand() / maxr;
}
b[r] = rand() / maxr;
}
int info;
int one = 1;
char N = 'N';
vector<int> ipiv(dim);
vector<double> a1(a);
vector<double> b1(b);
dgesv_(&dim, &one, a1.data(), &dim, ipiv.data(), b1.data(), &dim, &info);
print("B1",5,1,b1.data());
vector<double> a2(a);
vector<double> b2(b);
dgetrf_(&dim, &dim, a2.data(), &dim, ipiv.data(), &info);
dgetrs_(&N, &dim, &one, a2.data(), &dim, ipiv.data(), b2.data(), &dim, &info);
print("B2",5,1,b2.data());
double d_m_one = -1e0;
daxpy_(&dim,&d_m_one,b2.data(),&one,b1.data(),&one);
cout << endl << "diff = " << dnrm2_(&dim,b1.data(),&one) << endl;
return 0;
}
在我的PC中运行dim=5
情况,它给出了:
B1
---------
0.782902
3.35317
4.40269
-5.10512
-1.26615
---------
B2
---------
0.782902
3.35317
4.40269
-5.10512
-1.26615
---------
diff = 0
相关文章:
- 在VS代码中交叉编译Windows与Linux上的MinGW的SDL程序
- 如何为模板化对象创建模板向量?VS正在投掷C3203
- 数据成员SFINAE的C++17测试:gcc vs clang
- 为什么在Windows上的VS 2019和Clang 9中"size_t"在没有标题的情况下工作
- 在for循环中使用auto vs decltype(vec.size())来处理字符串的向量
- 正在VS调试器中监视映射条目
- Confusion: decltype vs std::function
- 将IBM Rhapsody模型集成到VS 2019中
- VS Code "command":"make"与终端窗口中的命令行"make"不同
- 使用VS Code和CMake Tools运行自定义命令
- 修改 VS Code 中的默认C++代码段
- 如何使用c++在VS 2019上运行SQL查询
- vs 2015 constexpr变量不恒定,但与2019相比还好吗
- 完美前进使用 std::forward vs RefRefCast
- 从VS 2015更新3更新到VS2015更新3 d后浮点计算行为不同的原因
- VS 2015 链接错误 无法构建依赖于 libcurl 的项目
- consteval wrapper vs. source_location
- VS Code C++:不准确的系统包括路径错误(wchar.h,boost/lambda/lambda.hpp)
- QStringList vs list<shared_ptr<QString>> 性能比较C++
- LAPACK dgetrs vs dgesv