Как начать использовать LAPACK в C ++?

Я новичок в вычислительной науке, и я уже изучил основные методы для интеграции, интерполяции, методы, такие как RK4, Numerov и т. Д. На C ++, но недавно мой профессор попросил меня научиться использовать 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, такие как Intel Math Kernel Library , очень быстрые.

Для моих собственных целей у меня есть собственные классы C ++ более высокого уровня, которые охватывают подпрограммы LAPACK. Многие научные библиотеки также используют LAPACK внизу. Иногда их проще использовать, но, на мой взгляд, очень важно понять инструмент под ним. Для этого я привел небольшой рабочий пример, написанный на C ++ и использующий LAPACK, для начала. Это работает в 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содержащий матрицу 3х3.

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 )

Комментарии:

• Вы, кажется, сбиты с толку схемой именования LAPACK. Краткое описание здесь .
• Интерфейс для подпрограммы DGEEV находится здесь . Вы должны быть в состоянии сравнить описание аргументов там с тем, что я сделал здесь.
• Обратите внимание на extern "C"раздел в верхней части, и я добавил к нему подчеркивание dgeev_. Это потому, что библиотека была написана и построена на Фортране, поэтому это необходимо, чтобы символы соответствовали при соединении. Это зависит от компилятора и системы, поэтому, если вы используете это в Windows, все это должно измениться.
• Некоторые люди могут предложить использовать интерфейс C для LAPACK . Возможно, они правы, но я всегда так делал.

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

Lapack написан на фортране, а API очень похож на фортран. Существует интерфейс C API для Lapack, который делает интерфейс немного менее болезненным, но использование Lapack из C ++ никогда не будет приятным.

В качестве альтернативы существует библиотека классов матрицы C ++ под названием Eigen, которая обладает многими возможностями 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;
}

В этом примере проблема собственных значений является тестом для функции Лапака dgeev. Вы можете просмотреть код FORTRAN и результаты для этого примера задачи dgeev и сделать свои собственные сравнения.

Вот еще один ответ в том же духе, что и выше.

Вы должны заглянуть в библиотеку линейной алгебры Armadillo C ++ .

Плюсы:

1. Синтаксис функции высокоуровневый (аналогичен синтаксису MATLAB). Так что никакого DGESVмумбо-юмбо, просто X = solve( A, B )(хотя за этими странно выглядящими именами функций LAPACK есть причина ...).
2. Реализует различные матричные разложения (LU, QR, собственные значения, SVD, Cholesky и т. Д.)
3. Это быстро при правильном использовании.
4. Это хорошо задокументировано .
5. Имеет поддержку разреженных матриц (вы захотите изучить их позже).
6. Вы можете связать его с вашими супер-оптимизированными библиотеками BLAS / LAPACK для оптимальной производительности.

Вот как будет выглядеть код @ BillGreene с Armadillo:

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