Comment commencer à utiliser LAPACK en c ++?

Je suis nouveau en science informatique et j'ai déjà appris des méthodes de base pour l'intégration, l'interpolation, des méthodes comme RK4, Numerov etc. sur c ++ mais récemment mon professeur m'a demandé d'apprendre à utiliser LAPACK pour résoudre des problèmes liés aux matrices. Comme par exemple trouver des valeurs propres d'une matrice complexe. Je n'ai jamais utilisé de bibliothèques tierces et j'écris presque toujours mes propres fonctions. Je fais des recherches depuis plusieurs jours, mais je ne trouve aucun guide amateur pour lapack. Tous sont écrits avec des mots que je ne comprends pas et je ne sais pas pourquoi l'utilisation de fonctions déjà écrites devrait être aussi compliquée. Ils sont pleins de mots comme zgeev, dtrsv, etc. et je suis frustré. Je veux juste coder quelque chose comme ce pseudo-code:

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

Je ne sais pas si je suis stupide ou amateur. Mais encore une fois, cela ne devrait pas être si difficile non? Je ne sais même pas si je dois utiliser LAPACK ou LAPACK ++. (J'écris des codes en c ++ et je ne connais pas Python ou FORTRAN) et comment les installer.

Je vais être en désaccord avec certaines des autres réponses et dire que je pense que trouver comment utiliser LAPACK est important dans le domaine de l'informatique scientifique.

Cependant, il existe une grande courbe d'apprentissage pour utiliser LAPACK. En effet, il est écrit à un niveau très bas. L'inconvénient de cela est qu'il semble très cryptique et peu agréable aux sens. L'avantage est que l'interface est sans ambiguïté et ne change fondamentalement jamais. De plus, les implémentations de LAPACK, telles que la bibliothèque Intel Math Kernel, sont très rapides.

Pour mes propres besoins, j'ai mes propres classes C ++ de niveau supérieur qui s'enroulent autour des sous-programmes LAPACK. De nombreuses bibliothèques scientifiques utilisent également LAPACK en dessous. Parfois, il est plus facile de simplement les utiliser, mais à mon avis, il est très utile de comprendre l'outil en dessous. À cette fin, j'ai fourni un petit exemple de travail écrit en C ++ à l'aide de LAPACK pour vous aider à démarrer. Cela fonctionne dans Ubuntu, avec le liblapack3package installé et d'autres packages nécessaires pour la construction. Il peut probablement être utilisé dans la plupart des distributions Linux, mais l'installation de LAPACK et sa liaison peuvent varier.

Voici le dossier 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;
}

Cela peut être construit en utilisant la ligne de commande

g++ -o test_lapack test_lapack.cpp -llapack

Cela produira un exécutable nommé test_lapack. J'ai configuré cela pour lire dans un fichier d'entrée de texte. Voici un fichier nommé matrix.txtcontenant une matrice 3x3.

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

Pour exécuter le programme, tapez simplement

./test_lapack matrix.txt

sur la ligne de commande, et la sortie doit être

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

Commentaires:

• Vous semblez découragé par le schéma de dénomination de LAPACK. Une brève description est ici .
• L'interface du sous-programme DGEEV est ici . Vous devriez pouvoir comparer la description des arguments à ce que j'ai fait ici.
• Notez la extern "C"section en haut et à laquelle j'ai ajouté un trait de soulignement dgeev_. C'est parce que la bibliothèque a été écrite et construite en Fortran, donc c'est nécessaire pour faire correspondre les symboles lors de la liaison. Cela dépend du compilateur et du système, donc si vous l'utilisez sous Windows, tout devra changer.
• Certaines personnes pourraient suggérer d'utiliser l' interface C pour LAPACK . Ils ont peut-être raison, mais je l'ai toujours fait de cette façon.

Je résiste généralement à dire aux gens ce que je pense qu'ils devraient faire plutôt que de répondre à leur question, mais dans ce cas, je vais faire une exception.

Lapack est écrit en FORTRAN et l'API est très semblable à FORTRAN. Il existe une API C pour Lapack qui rend l'interface un peu moins pénible, mais ce ne sera jamais une expérience agréable d'utiliser Lapack à partir de C ++.

Alternativement, il existe une bibliothèque de classes matricielles C ++ appelée Eigen qui possède de nombreuses capacités de Lapack, fournit des performances de calcul comparables aux meilleures implémentations Lapack et est très pratique à utiliser à partir de C ++. En particulier, voici comment votre exemple de code pourrait être écrit en utilisant 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;
}

Cet exemple de problème de valeur propre est un cas de test pour la fonction Lapack dgeev. Vous pouvez afficher le code FORTRAN et les résultats de cet exemple de problème dgeev et faire vos propres comparaisons.

Voici une autre réponse dans la même veine que ci-dessus.

Vous devriez regarder dans la bibliothèque d'algèbre linéaire Armadillo C ++ .

Avantages:

1. La syntaxe de la fonction est de haut niveau (similaire à celle de MATLAB). Donc pas de DGESVmumbo-jumbo, juste X = solve( A, B )(bien qu'il y ait une raison derrière ces noms de fonction LAPACK étrangement ...).
2. Implémente différentes décompositions matricielles (LU, QR, valeurs propres, SVD, Cholesky, etc.)
3. Il est rapide lorsqu'il est utilisé correctement.
4. C'est bien documenté .
5. Prend en charge les matrices clairsemées (vous voudrez les examiner plus tard).
6. Vous pouvez le lier à vos bibliothèques BLAS / LAPACK super optimisées pour des performances optimales.

Voici à quoi ressemblerait le code de @ BillGreene avec 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;
}