I am a beginner with LAPACK and C++/Fortran interfacing. I need to solve linear equations and eigenvalues problems using LAPACK/BLAS on Mac OS-X Lion. OS-X Lion provides optimized BLAS and LAPACK libraries (in /usr/lib) and I am linking these libraries instead of downloading them from netlib.
My program (reproduced below) is compiling and running fine, but it is giving me wrong answers. I have researched in the web and Stackoverflow and the issue may have to deal with how C++ and Fortran store arrays in differing formats (row major vs Column major). However, as you will see in my example, the simple array for my example should look identical in C++ and fortran. Anyway here goes.
Lets say we want to solve the following linear system:
x + y = 2
x - y = 0
The solution is (x,y) = (1,1). Now I tried to solve this using Lapack as follows
// LAPACK test code
#include<iostream>
#include<vector>
using namespace std;
extern "C" void dgetrs(char *TRANS, int *N, int *NRHS, double *A,
int *LDA, int *IPIV, double *B, int *LDB, int *INFO );
int main()
{
char trans = 'N';
int dim = 2;
int nrhs = 1;
int LDA = dim;
int LDB = dim;
int info;
vector<double> a, b;
a.push_back(1);
a.push_back(1);
a.push_back(1);
a.push_back(-1);
b.push_back(2);
b.push_back(0);
int ipiv[3];
dgetrs(&trans, &dim, &nrhs, & *a.begin(), &LDA, ipiv, & *b.begin(), &LDB, &info);
std::cout << "solution is:";
std::cout << "[" << b[0] << ", " << b[1] << ", " << "]" << std::endl;
std::cout << "Info = " << info << std::endl;
return(0);
}
This code was compiled as follows:
g++ -Wall -llapack -lblas lapacktest.cpp
On running this, however, I get the solution as (-2,2) which is obviously wrong. I have tried all combination of row/column re-arrangement of my matrix a
. Also observe the matrix representation of a
should be identical in row and column formats. I think getting this simple example to work will get me started with LAPACK and any help will be appreciated.
LAPACK ("Linear Algebra Package") is a standard software library for numerical linear algebra. It provides routines for solving systems of linear equations and linear least squares, eigenvalue problems, and singular value decomposition.
BLAS (Basic Linear Algebra Subprograms) is a library of vector, vector-vector, matrix-vector and matrix-matrix operations. LAPACK, a library of dense and banded matrix linear algebra routines such as solving linear systems, the eigenvalue- and singular value decomposition.
LAPACK is written in Fortran 90 and provides routines for solving systems of simultaneous linear equations, least-squares solutions of linear systems of equations, eigenvalue problems, and singular value problems.
You need to factor the matrix (by calling dgetrf
) before you can solve the system using dgetrs
. Alternatively, you can use the dgesv
routine, which does both steps for you.
By the way, you don't need to declare the interfaces yourself, they are in the Accelerate headers:
// LAPACK test code
#include <iostream>
#include <vector>
#include <Accelerate/Accelerate.h>
using namespace std;
int main()
{
char trans = 'N';
int dim = 2;
int nrhs = 1;
int LDA = dim;
int LDB = dim;
int info;
vector<double> a, b;
a.push_back(1);
a.push_back(1);
a.push_back(1);
a.push_back(-1);
b.push_back(2);
b.push_back(0);
int ipiv[3];
dgetrf_(&dim, &dim, &*a.begin(), &LDA, ipiv, &info);
dgetrs_(&trans, &dim, &nrhs, & *a.begin(), &LDA, ipiv, & *b.begin(), &LDB, &info);
std::cout << "solution is:";
std::cout << "[" << b[0] << ", " << b[1] << ", " << "]" << std::endl;
std::cout << "Info = " << info << std::endl;
return(0);
}
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With