Clapack_sgesv ATLAS с порядком хранения основной строки

Question

У меня проблема с методом ATLAS clapack_sgesv (соответствующий метод FORTRAN: sgesv.f), когда я пытаюсь использовать его с матрицами, хранящимися в порядке хранения основной строки.

Я использую Eigen3 для большинства задач линейной алгебры в моем приложении, но недавно начал заменять некоторые внутренние подпрограммы Eigen вызовами методов ATLAS cblas и clapack. Мое приложение должно поддерживать переключение порядка хранения матрицы на мажорную строку путем определения Eigen's EIGEN_DEFAULT_TO_ROW_MAJOR флаг. Это, конечно, работает из коробки с методами Эйгена, но требует различных путей кода для clapack_ звонки. Я столкнулся с проблемой при замене Эйгена .partialPivLu().solve() позвонить с ATLAS' clapack_sgesv метод. Вот минимальный пример кода, который иллюстрирует проблему:

#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;
}

Я собираю это с g++ -std=c++11 -Wall -Wextra -g -llapack -lcblas -latlas, Выше clapack_sgesv вызов дает те же результаты, что и решатель Эйгена, пока EIGEN_DEFAULT_TO_ROW_MAJOR не определено.

----- 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

Если я определю это, результаты ATLAS неверны.

----- 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

Моим первым подозрением было, конечно, то, что я что-то напутал с clapack_sgesv() вызов. Но кроме установки порядка хранения и переключения начальных размеров с количества строк на число столбцов, кажется, нет необходимости.

Еще одна действительно запутанная вещь, которую я заметил, заключается в следующем: когда я пытаюсь решить только для одной правой стороны, clapack_sgesv() вызов не удался с Parameter 8 to routine clapack_sgesv was incorrect, ldb must be >= MAX(N,1): ldb=1 N=4, Это не имеет никакого смысла математически.

Я подозреваю, что моя ошибка должна быть очевидной, но я ее не вижу.

Что не так с моим clapack_sgesv() вызов, который вызывает его в порядке хранения основной строки?

Answer 1

Я нашел свою ошибку. Как объясняется в FAQ по ATLAS, правая часть не обрабатывается как матрица, а представляет собой совокупность векторов-столбцов, находящихся рядом в памяти. Это не имеет значения, если порядок хранения является главным по столбцу, но это относится к порядку хранения по основной строке, поскольку элементы вектора столбца больше не являются смежными в памяти. Это работает, если всегда хранить RHS и "матрицу" решения в формате столбца.