Wie starte ich LAPACK in c ++?

Ich bin neu in der Computerwissenschaft und habe bereits grundlegende Methoden für Integration, Interpolation, Methoden wie RK4, Numerov usw. in C ++ gelernt. Kürzlich hat mich mein Professor gebeten, zu lernen, wie man LAPACK zur Lösung von Problemen im Zusammenhang mit Matrizen verwendet. Wie zum Beispiel das Finden von Eigenwerten einer komplexen Matrix. Ich habe noch nie Bibliotheken von Drittanbietern verwendet und schreibe fast immer meine eigenen Funktionen. Ich habe mehrere Tage lang gesucht, kann aber keinen amateurfreundlichen Leitfaden für Lapack finden. Alle von ihnen sind in Worten geschrieben, die ich nicht verstehe, und ich weiß nicht, warum die Verwendung bereits geschriebener Funktionen so kompliziert sein sollte. Sie sind voller Wörter wie zgeev, dtrsv usw. und ich bin frustriert. Ich möchte nur so etwas wie diesen Pseudocode codieren:

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

Ich weiß nicht, ob ich albern oder amateurhaft bin. Aber auch das sollte nicht so schwer sein, oder? Ich weiß nicht einmal, ob ich LAPACK oder LAPACK ++ verwenden soll. (Ich schreibe Codes in C ++ und habe keine Kenntnisse über Python oder FORTRAN) und wie man sie installiert.

Ich werde einigen der anderen Antworten nicht zustimmen und sagen, dass ich glaube, dass es im Bereich des wissenschaftlichen Rechnens wichtig ist , herauszufinden, wie man LAPACK verwendet .

Die Verwendung von LAPACK ist jedoch mit einer großen Lernkurve verbunden. Dies liegt daran, dass es auf einer sehr niedrigen Ebene geschrieben ist. Der Nachteil davon ist, dass es sehr kryptisch und für die Sinne nicht angenehm erscheint. Der Vorteil ist, dass die Schnittstelle eindeutig ist und sich grundsätzlich nie ändert. Darüber hinaus sind Implementierungen von LAPACK wie der Intel Math Kernel Library sehr schnell.

Für meine eigenen Zwecke habe ich meine eigenen übergeordneten C ++ - Klassen, die sich um LAPACK-Subroutinen drehen. Viele wissenschaftliche Bibliotheken verwenden auch LAPACK darunter. Manchmal ist es einfacher, sie einfach zu verwenden, aber meiner Meinung nach ist es sehr wertvoll, das darunter liegende Tool zu verstehen. Zu diesem Zweck habe ich ein kleines Arbeitsbeispiel bereitgestellt, das in C ++ mit LAPACK geschrieben wurde, um Ihnen den Einstieg zu erleichtern. Dies funktioniert in Ubuntu mit dem liblapack3installierten Paket und anderen zum Erstellen erforderlichen Paketen. Es kann wahrscheinlich in den meisten Linux-Distributionen verwendet werden, aber die Installation von LAPACK und die Verknüpfung damit können variieren.

Hier ist die Datei 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;
}

Dies kann über die Befehlszeile erstellt werden

g++ -o test_lapack test_lapack.cpp -llapack

Dies erzeugt eine ausführbare Datei mit dem Namen test_lapack. Ich habe dies so eingerichtet, dass es eine Texteingabedatei einliest. Hier ist eine Datei mit dem Namen matrix.txt3x3-Matrix.

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

Um das Programm auszuführen, geben Sie einfach ein

./test_lapack matrix.txt

an der Kommandozeile, und die Ausgabe sollte sein

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

Bemerkungen:

• Das Namensschema für LAPACK scheint Sie zu stören. Eine kurze Beschreibung finden Sie hier .
• Die Schnittstelle für das DGEEV-Unterprogramm befindet sich hier . Sie sollten in der Lage sein, die Beschreibung der Argumente dort mit dem zu vergleichen, was ich hier getan habe.
• Beachten Sie den extern "C"Abschnitt oben, dem ich einen Unterstrich hinzugefügt habe dgeev_. Das liegt daran, dass die Bibliothek in Fortran geschrieben und erstellt wurde. Dies ist also erforderlich, damit die Symbole beim Verknüpfen übereinstimmen. Dies ist vom Compiler und vom System abhängig. Wenn Sie dies unter Windows verwenden, muss sich alles ändern.
• Einige Leute schlagen möglicherweise vor, die C-Schnittstelle für LAPACK zu verwenden . Sie mögen Recht haben, aber ich habe es immer so gemacht.

Normalerweise weigere ich mich, den Leuten zu sagen, was sie tun sollen, anstatt ihre Frage zu beantworten, aber in diesem Fall mache ich eine Ausnahme.

Lapack ist in FORTRAN geschrieben und die API ist sehr FORTRAN-ähnlich. Es gibt eine C-API für Lapack, die die Benutzeroberfläche etwas weniger schmerzhaft macht, aber es wird niemals eine angenehme Erfahrung sein, Lapack aus C ++ zu verwenden.

Alternativ gibt es eine C ++ - Matrixklassenbibliothek namens Eigen , die viele der Funktionen von Lapack bietet, eine mit den besseren Lapack-Implementierungen vergleichbare Rechenleistung bietet und von C ++ aus sehr bequem zu verwenden ist. Hier erfahren Sie insbesondere, wie Ihr Beispielcode mit Eigen geschrieben werden kann

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

Dieses beispielhafte Eigenwertproblem ist ein Testfall für die Lapack-Funktion dgeev. Sie können den FORTRAN-Code und die Ergebnisse für dieses Problem anzeigen und Ihre eigenen Vergleiche anstellen .

Hier ist eine weitere Antwort in der gleichen Weise wie oben.

Sie sollten in die Armadillo C ++ - Bibliothek für lineare Algebra schauen .

Vorteile:

1. Die Funktionssyntax ist auf hoher Ebene (ähnlich der von MATLAB). Also kein DGESVHokuspokus, nur X = solve( A, B )(obwohl es einen Grund für diese seltsam aussehenden LAPACK-Funktionsnamen gibt ...).
2. Implementiert verschiedene Matrixzerlegungen (LU, QR, Eigenwerte, SVD, Cholesky usw.)
3. Es ist schnell, wenn es richtig verwendet wird.
4. Es ist gut dokumentiert .
5. Unterstützt spärliche Matrizen (Sie werden diese später untersuchen wollen).
6. Sie können es mit Ihren superoptimierten BLAS / LAPACK-Bibliotheken verknüpfen, um eine optimale Leistung zu erzielen.

So würde der Code von @ BillGreene mit Armadillo aussehen:

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