我是计算科学的新手，我已经学习了 C++ 上的积分、插值、RK4、Numerov 等方法的基本方法，但最近我的教授让我学习如何使用 LAPACK 解决与矩阵相关的问题。例如寻找复杂矩阵的特征值。我从未使用过第三方库，而且我几乎总是编写自己的函数。我已经搜索了几天，但找不到任何适合业余爱好者的 lapack 指南。所有这些都是用我不懂的文字写的，我不知道为什么使用已经写好的函数会这么复杂。他们充满了zgeev，dtrsv等词，我很沮丧。我只想编写类似这样的伪代码：
#include <lapack:matrix>
int main(){
LapackComplexMatrix A(n,n);
for...
for...
cin>>A(i,j);
cout<<LapackEigenValues(A);
return 0;
}
我不知道我是傻还是业余。但是，这不应该那么难吗？我什至不知道我应该使用 LAPACK 还是 LAPACK++。（我用 C++ 编写代码，不了解 Python 或 FORTRAN）以及如何安装它们。
我将不同意其他一些答案，并说我相信弄清楚如何使用 LAPACK在科学计算领域很重要。
但是，使用 LAPACK 有很大的学习曲线。这是因为它是在非常低的级别上编写的。这样做的缺点是它看起来很神秘，而且感觉不愉快。它的优点是界面明确，基本不会改变。此外，LAPACK 的实施，例如英特尔数学内核库，速度非常快。
出于我自己的目的，我有自己的更高级别的 C++ 类，这些类包含 LAPACK 子例程。许多科学图书馆也在底层使用 LAPACK。有时只使用它们会更容易，但在我看来，理解下面的工具有很多价值。为此，我提供了一个使用 LAPACK 用 C++ 编写的小型工作示例来帮助您入门。这适用于 Ubuntu，安装了liblapack3软件包，以及其他必要的构建软件包。它可能可以在大多数 Linux 发行版中使用，但是 LAPACK 的安装和链接可能会有所不同。
这是文件test_lapack.cpp
#include <iostream>
#include <fstream>
using namespace std;
// dgeev_ is a symbol in the LAPACK library files
extern "C" {
extern int dgeev_(char*,char*,int*,double*,int*,double*, double*, double*, int*, double*, int*, double*, int*, int*);
}
int main(int argc, char** argv){
// check for an argument
if (argc<2){
cout << "Usage: " << argv[0] << " " << " filename" << endl;
return -1;
}
int n,m;
double *data;
// read in a text file that contains a real matrix stored in column major format
// but read it into row major format
ifstream fin(argv[1]);
if (!fin.is_open()){
cout << "Failed to open " << argv[1] << endl;
return -1;
}
fin >> n >> m; // n is the number of rows, m the number of columns
data = new double[n*m];
for (int i=0;i<n;i++){
for (int j=0;j<m;j++){
fin >> data[j*n+i];
}
}
if (fin.fail() || fin.eof()){
cout << "Error while reading " << argv[1] << endl;
return -1;
}
fin.close();
// check that matrix is square
if (n != m){
cout << "Matrix is not square" <<endl;
return -1;
}
// allocate data
char Nchar='N';
double *eigReal=new double[n];
double *eigImag=new double[n];
double *vl,*vr;
int one=1;
int lwork=6*n;
double *work=new double[lwork];
int info;
// calculate eigenvalues using the DGEEV subroutine
dgeev_(&Nchar,&Nchar,&n,data,&n,eigReal,eigImag,
vl,&one,vr,&one,
work,&lwork,&info);
// check for errors
if (info!=0){
cout << "Error: dgeev returned error code " << info << endl;
return -1;
}
// output eigenvalues to stdout
cout << "--- Eigenvalues ---" << endl;
for (int i=0;i<n;i++){
cout << "( " << eigReal[i] << " , " << eigImag[i] << " )\n";
}
cout << endl;
// deallocate
delete [] data;
delete [] eigReal;
delete [] eigImag;
delete [] work;
return 0;
}
这可以使用命令行构建
g++ -o test_lapack test_lapack.cpp -llapack
这将生成一个名为test_lapack. 我已将其设置为读取文本输入文件。这是一个名为matrix.txt包含 3x3 矩阵的文件。
3 3
-1.0 -8.0 0.0
-1.0 1.0 -5.0
3.0 0.0 2.0
要运行程序，只需键入
./test_lapack matrix.txt
在命令行，输出应该是
--- Eigenvalues ---
( 6.15484 , 0 )
( -2.07742 , 3.50095 )
( -2.07742 , -3.50095 )
评论：
extern "C"顶部的部分，我在dgeev_. 这是因为该库是用 Fortran 编写和构建的，因此在链接时需要使符号匹配。这取决于编译器和系统，所以如果你在 Windows 上使用它，它都必须改变。这是与上述相同的另一个答案。
您应该查看Armadillo C++ 线性代数库。
优点：
DGESV胡说八道，只是X = solve( A, B )（尽管那些看起来很奇怪的 LAPACK 函数名称背后有一个原因......）。这是@BillGreene 的代码在犰狳中的样子：
#include <iostream>
#include <armadillo>
using namespace std;
using namespace arma;
int main()
{
const int k = 4;
mat A = zeros<mat>(k,k) // mat == Mat<double>
// with the << operator...
A <<
0.35 << 0.45 << -0.14 << -0.17 << endr
0.09 << 0.07 << -0.54 << 0.35 << endr
-0.44 << -0.33 << -0.03 << 0.17 << endr
0.25 << -0.32 << -0.13 << 0.11 << endr;
// but using an initializer list is faster
A = { {0.35, 0.45, -0.14, -0.17},
{0.09, 0.07, -0.54, 0.35},
{-0.44, -0.33, -0.03, 0.17},
{0.25, -0.32, -0.13, 0.11} };
cx_vec eigval; // eigenvalues may well be complex
cx_mat eigvec;
// eigenvalue decomposition for general dense matrices
eig_gen(eigval, eigvec, A);
std::cout << eigval << std::endl;
return 0;
}
我通常拒绝告诉人们我认为他们应该做什么，而不是回答他们的问题，但在这种情况下，我将破例。
Lapack 是用 FORTRAN 编写的，API 非常类似于 FORTRAN。Lapack 有一个 C API，可以让界面稍微不那么痛苦，但从 C++ 中使用 Lapack 永远不会是一种愉快的体验。
或者，有一个名为Eigen的 C++ 矩阵类库，它具有 Lapack 的许多功能，提供与更好的 Lapack 实现相当的计算性能，并且在 C++ 中使用非常方便。特别是，这里是使用 Eigen 编写示例代码的方式
#include <iostream>
using std::cout;
using std::endl;
#include <Eigen/Eigenvalues>
int main()
{
const int n = 4;
Eigen::MatrixXd a(n, n);
a <<
0.35, 0.45, -0.14, -0.17,
0.09, 0.07, -0.54, 0.35,
-0.44, -0.33, -0.03, 0.17,
0.25, -0.32, -0.13, 0.11;
Eigen::EigenSolver<Eigen::MatrixXd> es;
es.compute(a);
Eigen::VectorXcd ev = es.eigenvalues();
cout << ev << endl;
}
此示例特征值问题是 Lapack 函数的测试用例
dgeev。您可以查看此问题
dgeev 示例的 FORTRAN 代码和结果并进行自己的比较。