- Java 双重比较
- java - 比较器与 Apache BeanComparator
- Objective-C 完成 block 导致额外的方法调用?
- database - RESTful URI 是否应该公开数据库主键?
在 LAPACK 文档中,
dgesv 求解一般线性方程组 AX=B。
dgetrs 使用 DGETRF 计算的 LU 分解求解一般线性方程组 AX=B、AT X=B 或 AH 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
在 0 秒内生成 5 x 5 矩阵
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
。source code dgesv
非常简单。它主要包含对dgetrf
和dgetrs
函数的调用。
示例中的结果不同,因为 dgetrf
更改了 A
矩阵。在 lapack's reference它指出:
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;
}
在我运行 dim=5
的 PC 中,它给出了:
B1
---------
0.782902
3.35317
4.40269
-5.10512
-1.26615
---------
B2
---------
0.782902
3.35317
4.40269
-5.10512
-1.26615
---------
diff = 0
关于c++ - LAPACK dgetrs 与 dgesv,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/36063993/
如果我有一个 1,000 x 1,000 的方阵,Lapack 可以计算这个矩阵的特征向量和特征值吗?如果可以,需要多长时间?对于 10,000 x 10,000 矩阵甚至 1,000,000 x 1
我正在使用 LAPACK 的 zheev(在英特尔 MKL 中)。我得到了 int INFO=99。我一直在互联网上搜索这对应的内容,但找不到包含所有整数错误代码及其含义列表的文档。 有没有人有指向
通过以下方式安装 lapack 后: yum install lapack lapack-devel 并重新启动 httpd 服务我仍然得到 Fatal error: Class 'Lapack' n
Lapack 很可能没有任何计算行列式的例程。但是,我们可以使用 LU、QR 或 SVD 分解来计算它。我更喜欢使用 LU 分解。现在 lapack 使用一些 dgetrf 子程序将矩阵 A 分解为带
根据官方用户指南,sgelsd用于解决最小二乘问题 min_x || b - Ax ||_2 并允许矩阵 A 为矩形且秩亏。并且根据sgelsd源码中的接口(interface)描述,b作为输入输出参
我用 Eigen 和 实现了一段代码。我希望 Eigen 使用 BLAS 和 LAPACK 。 我看过here ,这是可能的,但我不知道如何或将这些值/指令放在代码中的位置。 我必须在某处指定值 EI
简介:我用 C++ 开发了一个应用程序,它在 Windows 上使用 LAPACK(LAPACKE) 和 MPI。在 Windows 中工作正常(编译和链接通过 Code::Blocks IDE 处理
我最近从 Linux 切换到 Mac OS。我需要 BLAS 和 LAPACK 来做一些计算。通过查看 BLAS 的维基百科,我了解到这两个库已经在 Mac OS 中实现。不过,据说 Apple's
我基于以下链接为我的 Visual Studio 2008 构建了 LAPACKE 的 DLL 和库: http://icl.cs.utk.edu/lapack-for-windows/lapack/
我很好奇用于在 MATLAB 中计算 SVD 的 DGESVD 函数。据我从 Gene H. Golub 和 Charles F. Van Loan 的“矩阵计算”中可以看出,使用了两种可能的双对角化
我是数值线性代数的新手,我刚刚开始使用 LAPACK 和 BLAS。 是否有可以在打包存储和完整存储之间复制/转换对称矩阵的例程? 我找到了 dtrttp ,我可以用它来将 double 全对称矩阵转
根据我的理解,需要进行分解/因式分解(LU、QR、Cholesky 等),然后基于因式分解进行矩阵逆计算。还有其他方法可以解决这个问题吗(我试图弄清楚我是否可以坚持使用 CULAtools 试用版中免
我正在尝试编写一个函数,该函数可以为代表性不足的方程组生成单一解(例如,描述该系统的矩阵宽大于高)。为了做到这一点,我一直在 LAPACK 文档中寻找一种将矩阵行归约到它的归约梯队形式的方法,类似于
我是使用 LAPACK 例程的新手,所以我对它们不是很了解,我想在并行化循环(openmp)中使用它们。 我使用 Ubuntu 14.04LTS 并使用我的包管理器安装了 LAPACK。安装的版本是:
我在允许我求逆矩阵的 c 代码中使用 LAPACK。更准确地说,我使用 dgetrf_ 然后使用 dgetri_ 进行反演。 但是当我处理大矩阵并且我不知道矩阵是否可逆时,我浪费了很多时间来反转不可逆
我想使用Fortran和LAPACK对角化一个实对称矩阵。 LAPACK基本上提供了两个例程,一个例程在完整矩阵上运行,另一个例程在打包存储中的矩阵上运行。虽然后者肯定会使用较少的内存,但我想知道关于
我正在使用 LAPACK 库中的 DSYEV 和 DSYEVD 来查找特征值和特征向量(编译语法: gfortran -llapack )。但是,我发现特定矩阵的错误特征值(-0.44,0.35,0.
我想计算对称矩阵的特征值,并希望使用 C++ 中英特尔 MKL 库中的 LAPACKE_dsyev 函数来计算。但我对矩阵需要传递的形式有点困惑。 来自文档https://software.intel
我使用链接:./configure --with-blas-incdir="-L/home/moritz/build/CoinIpopt_test/ThirdParty/openblas/includ
我想安装最新版本的numpy(Python的数值库),版本(v1.6.1)还没有在Ubuntu Oneiric repositories .当我继续手动安装它时,我阅读了 INSTALL numpy
我是一名优秀的程序员,十分优秀!