I'm having a problem with ATLAS' method clapack_sgesv (corresponding FORTRAN method: sgesv.f) when I try to use it with matrices stored in row-major storage order.
I use Eigen3 for most linear algebra tasks in my application but have recently started to replace some internal Eigen routines with calls to ATLAS' cblas and clapack methods. My application must support switching the matrix storage order to row-major by defining Eigen's EIGEN_DEFAULT_TO_ROW_MAJOR flag. This of course works out of the box with Eigen's methods, but requires different code paths for the clapack_ calls. I've run into a problem when replacing Eigen's .partialPivLu().solve() call with ATLAS' clapack_sgesv method. Here is a minimal code example that illustrates the problem:
#include <iostream>
#define EIGEN_DEFAULT_TO_ROW_MAJOR
#include <eigen3/Eigen/Eigen>
extern "C" {
#include <clapack.h>
}
using namespace std;
int main()
{
Eigen::MatrixXf A( 4, 4 );
A << 0.680375 , 0.823295 , -0.444451 , -0.270431 ,
-0.211234 , -0.604897 , 0.10794 , 0.0268018 ,
0.566198 , -0.329554 , -0.0452059 , 0.904459 ,
0.59688 , 0.536459 , 0.257742 , 0.83239 ;
Eigen::MatrixXf B( 4, 4 );
B << 0.271423 , -0.967399 , -0.686642 , 0.997849 ,
0.434594 , -0.514226 , -0.198111 , -0.563486 ,
-0.716795 , -0.725537 , -0.740419 , 0.0258648 ,
0.213938 , 0.608353 , -0.782382 , 0.678224 ;
cout << "----- Eigen" << endl;
cout << "A = " << endl << A << endl;
cout << "B = " << endl << B << endl;
Eigen::MatrixXf X = A.partialPivLu().solve( B );
cout << "X = " << endl << X << endl;
cout << "AX = " << endl << A * X << endl;
cout << "----- ATLAS" << endl;
Eigen::VectorXi ipiv( 4 );
clapack_sgesv(
#ifdef EIGEN_DEFAULT_TO_ROW_MAJOR
CblasRowMajor,
#else
CblasColMajor,
#endif
A.rows(),
B.cols(),
A.data(),
#ifdef EIGEN_DEFAULT_TO_ROW_MAJOR
A.cols(),
#else
A.rows(),
#endif
ipiv.data(),
B.data(),
#ifdef EIGEN_DEFAULT_TO_ROW_MAJOR
B.cols()
#else
B.rows()
#endif
);
cout << "piv = " << ipiv.transpose() << endl;
cout << "LU = " << endl << A << endl;
cout << "X =" << endl << B << endl;
return 0;
}
I compile this with g++ -std=c++11 -Wall -Wextra -g -llapack -lcblas -latlas. The above clapack_sgesv call gives the same results as Eigen's solver as long as EIGEN_DEFAULT_TO_ROW_MAJOR is not defined.
----- Eigen
A =
0.680375 0.823295 -0.444451 -0.270431
-0.211234 -0.604897 0.10794 0.0268018
0.566198 -0.329554 -0.0452059 0.904459
0.59688 0.536459 0.257742 0.83239
B =
0.271423 -0.967399 -0.686642 0.997849
0.434594 -0.514226 -0.198111 -0.563486
-0.716795 -0.725537 -0.740419 0.0258648
0.213938 0.608353 -0.782382 0.678224
X =
4.29176 -3.45693 -3.46864 0.547927
-1.3688 2.04333 1.13806 0.735351
5.6716 -0.593909 -2.65158 -0.0154493
-3.69446 2.07672 1.6349 -0.0472447
AX =
0.271423 -0.967399 -0.686642 0.997849
0.434594 -0.514226 -0.198111 -0.563486
-0.716796 -0.725537 -0.740419 0.0258648
0.213938 0.608353 -0.782382 0.678224
----- ATLAS
piv = 0 2 3 3
LU =
0.680375 0.823295 -0.444451 -0.270431
0.832185 -1.01469 0.32466 1.12951
0.877281 0.183112 0.588201 0.862807
-0.310467 0.344235 -0.241085 -0.237964
X =
4.29176 -3.45694 -3.46864 0.547927
-1.3688 2.04333 1.13806 0.735351
5.6716 -0.593909 -2.65158 -0.0154493
-3.69446 2.07672 1.6349 -0.0472447
If I define it, the ATLAS' results are wrong.
----- Eigen
[... same as above ...]
----- ATLAS
piv = 1 1 3 3
LU =
0.823295 0.826405 -0.328474 -0.539844
-0.604897 0.288656 -0.595488 -0.757338
-0.329554 0.838543 1.29555 0.31797
0.536459 0.153548 1.10004 0.313854
X =
-2.21567 2.33841 -0.554441 1.45218
-2.60368 1.14776 -3.83383 1.63747
-5.05167 2.4991 -3.36881 3.08596
6.03571 -1.84576 8.32067 -4.90008
My first suspicion was of course that I had messed something up with the clapack_sgesv() call. But other than setting the storage order and switching the leading dimensions from the number of rows to the number of cols there seems nothing to be necessary.
Another really confusing thing I've noticed is the following: When I try to solve only for a single right-hand-side, the clapack_sgesv() call fails with Parameter 8 to routine clapack_sgesv was incorrect, ldb must be >= MAX(N,1): ldb=1 N=4. This does not make any sense mathematically.
I suspect that my error must be kind of obvious, but I don't see it.
What is wrong with my clapack_sgesv() call that causes it to fail in row-major storage order?
I've found my mistake. As explained in the ATLAS FAQ the right hand side is not treated like a matrix, but a collection of column-vectors that are adjacent in memory. This does not make a difference if the storage order is column-major, but it does for row-major storage order, because the elements of a column vector are no longer adjacent in memory. It works if one always stores RHS and solution "matrix" in column-major format.