Jak zacząć korzystać z LAPACK w c ++?

Jestem nowy w informatyce i nauczyłem się już podstawowych metod integracji, interpolacji, metod takich jak RK4, Numerov itp. Na c ++, ale ostatnio mój profesor poprosił mnie o nauczenie się, jak używać LAPACK do rozwiązywania problemów związanych z macierzami. Jak na przykład znajdowanie wartości własnych złożonej macierzy. Nigdy nie korzystałem z bibliotek stron trzecich i prawie zawsze piszę własne funkcje. Szukałem od kilku dni, ale nie znalazłem przyjaznego dla amatora przewodnika po lapacku. Wszystkie są napisane słowami, których nie rozumiem i nie wiem, dlaczego korzystanie z już napisanych funkcji powinno być tak skomplikowane. Są pełne słów takich jak zgeev, dtrsv itp. I jestem sfrustrowany. Chcę tylko zakodować coś takiego jak ten pseudo-kod:

#include <lapack:matrix>
int main(){
  LapackComplexMatrix A(n,n);
  for...
   for...
    cin>>A(i,j);
  cout<<LapackEigenValues(A);
  return 0;
}

Nie wiem, czy jestem głupi, czy amatorski. Ale znowu, to nie powinno być takie trudne, prawda? Nie wiem nawet, czy powinienem używać LAPACK lub LAPACK ++. (Piszę kody w c ++ i nie znam znajomości Python ani FORTRAN) i jak je zainstalować.

Nie zgadzam się z niektórymi innymi odpowiedziami i powiem, że uważam, że zastanawianie się, jak korzystać z LAPACK, jest ważne w dziedzinie informatyki naukowej.

Istnieje jednak duża krzywa uczenia się przy korzystaniu z LAPACK. Jest tak, ponieważ jest napisane na bardzo niskim poziomie. Wadą tego jest to, że wydaje się bardzo tajemnicze i nieprzyjemne dla zmysłów. Zaletą jest to, że interfejs jest jednoznaczny i zasadniczo nigdy się nie zmienia. Ponadto implementacje LAPACK, takie jak Intel Math Kernel Library, są naprawdę szybkie.

Na własne potrzeby mam własne klasy C ++ wyższego poziomu, które otaczają podprogramy LAPACK. Wiele bibliotek naukowych używa również LAPACK poniżej. Czasami łatwiej jest po prostu z nich korzystać, ale moim zdaniem zrozumienie narzędzia pod spodem ma dużą wartość. W tym celu podałem mały działający przykład napisany w C ++ przy użyciu LAPACK, aby zacząć. Działa to w Ubuntu, z liblapack3zainstalowanym pakietem i innymi niezbędnymi pakietami do budowania. Prawdopodobnie można go używać w większości dystrybucji Linuksa, ale instalacja LAPACK i łączenie z nim może się różnić.

Oto plik 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;
}

Można to zbudować za pomocą wiersza poleceń

g++ -o test_lapack test_lapack.cpp -llapack

Spowoduje to utworzenie pliku wykonywalnego o nazwie test_lapack. Skonfigurowałem to do odczytu w pliku wejściowym tekstu. Oto plik o nazwie matrix.txtzawierający matrycę 3x3.

3 3
-1.0 -8.0  0.0
-1.0  1.0 -5.0
 3.0  0.0  2.0

Aby uruchomić program, po prostu wpisz

./test_lapack matrix.txt

w wierszu polecenia, a wynik powinien być

--- Eigenvalues ---
( 6.15484 , 0 )
( -2.07742 , 3.50095 )
( -2.07742 , -3.50095 )

Komentarze:

• Wygląda na zaskoczonego schematem nazewnictwa LAPACK. Krótki opis jest tutaj .
• Interfejs podprogramu DGEEV znajduje się tutaj . Powinieneś być w stanie porównać opis argumentów tam z tym, co tutaj zrobiłem.
• Zwróć uwagę na extern "C"sekcję u góry i do której dodałem podkreślenie dgeev_. Wynika to z faktu, że biblioteka została napisana i zbudowana w Fortranie, dlatego konieczne jest dopasowanie symboli podczas łączenia. Jest to zależne od kompilatora i systemu, więc jeśli używasz tego w systemie Windows, wszystko będzie musiało się zmienić.
• Niektóre osoby mogą sugerować użycie interfejsu C do LAPACK . Mogą mieć rację, ale zawsze tak robiłem.

Zwykle powstrzymuję się od mówienia ludziom, co moim zdaniem powinni zrobić, zamiast odpowiadania na ich pytania, ale w tym przypadku zrobię wyjątek.

Lapack jest napisany w FORTRAN, a API jest bardzo podobne do FORTRAN. Lapack posiada C API, który sprawia, że ​​interfejs jest nieco mniej bolesny, ale korzystanie z Lapacka z C ++ nigdy nie będzie przyjemne.

Alternatywnie istnieje biblioteka klas macierzy C ++ o nazwie Eigen, która ma wiele możliwości Lapacka, zapewnia wydajność obliczeniową porównywalną z lepszymi implementacjami Lapacka i jest bardzo wygodna w użyciu z C ++. W szczególności oto, w jaki sposób można napisać przykładowy kod za pomocą 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;
}

Ten przykładowy problem wartości własnej jest przypadkiem testowym dla funkcji Lapacka dgeev. Możesz przeglądać kod FORTRAN i wyniki dla tego problemu dgeev i tworzyć własne porównania.

Oto kolejna odpowiedź w tym samym stylu co powyżej.

Powinieneś zajrzeć do biblioteki algebry liniowej Armadillo C ++ .

Plusy:

1. Składnia funkcji jest na wysokim poziomie (podobna do MATLAB). Więc nie ma DGESVmumbo-jumbo, tylko X = solve( A, B )(chociaż istnieje powód tych dziwnie wyglądających nazw funkcji LAPACK ...).
2. Realizuje różne dekompozycje macierzy (LU, QR, wartości własne, SVD, Cholesky itp.)
3. Jest szybki, gdy jest właściwie stosowany.
4. Jest dobrze udokumentowany .
5. Obsługuje rzadkie matryce (będziesz chciał zajrzeć do nich później).
6. Możesz połączyć go z superoptymalizowanymi bibliotekami BLAS / LAPACK, aby uzyskać optymalną wydajność.

Oto jak wyglądałby kod @ BillGreene w 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;
}