¿Cómo comenzar a usar LAPACK en c ++?

Soy nuevo en la ciencia computacional y ya he aprendido métodos básicos para la integración, la interpolación, métodos como RK4, Numerov, etc. en c ++, pero recientemente mi profesor me pidió que aprendiera cómo usar LAPACK para resolver problemas relacionados con las matrices. Como por ejemplo encontrar valores propios de una matriz compleja. Nunca he usado bibliotecas de terceros y casi siempre escribo mis propias funciones. He estado buscando por varios días, pero no puedo encontrar ninguna guía amigable para los aficionados a lapack. Todos están escritos en palabras que no entiendo y no sé por qué el uso de funciones ya escritas debería ser tan complicado. Están llenas de palabras como zgeev, dtrsv, etc. y estoy frustrado. Solo quiero codificar algo como este pseudocódigo:

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

No sé si estoy siendo tonto o aficionado. Pero de nuevo, esto no debería ser tan difícil, ¿verdad? Ni siquiera sé si debo usar LAPACK o LAPACK ++. (Escribo códigos en c ++ y no tengo conocimiento de Python o FORTRAN) y cómo instalarlos.

No estaré de acuerdo con algunas de las otras respuestas y diré que creo que descubrir cómo usar LAPACK es importante en el campo de la computación científica.

Sin embargo, hay una gran curva de aprendizaje para usar LAPACK. Esto se debe a que está escrito en un nivel muy bajo. La desventaja de eso es que parece muy críptico y no es agradable para los sentidos. La ventaja de esto es que la interfaz es inequívoca y, básicamente, nunca cambia. Además, las implementaciones de LAPACK, como la Intel Math Kernel Library son realmente rápidas.

Para mis propios fines, tengo mis propias clases de C ++ de nivel superior que se ajustan a las subrutinas LAPACK. Muchas bibliotecas científicas también usan LAPACK debajo. A veces es más fácil usarlos, pero en mi opinión hay mucho valor en comprender la herramienta que se encuentra debajo. Con ese fin, proporcioné un pequeño ejemplo de trabajo escrito en C ++ usando LAPACK para comenzar. Esto funciona en Ubuntu, con el liblapack3paquete instalado y otros paquetes necesarios para la construcción. Probablemente se pueda usar en la mayoría de las distribuciones de Linux, pero la instalación de LAPACK y el enlace en su contra pueden variar.

Aquí está el archivo 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;
}

Esto se puede construir usando la línea de comando

g++ -o test_lapack test_lapack.cpp -llapack

Esto producirá un ejecutable llamado test_lapack. He configurado esto para leer en un archivo de entrada de texto. Aquí hay un archivo llamado que matrix.txtcontiene una matriz 3x3.

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

Para ejecutar el programa simplemente escriba

./test_lapack matrix.txt

en la línea de comando, y la salida debería ser

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

Comentarios:

• Parece confundido por el esquema de nombres para LAPACK. Una breve descripción está aquí .
• La interfaz para la subrutina DGEEV está aquí . Debería poder comparar la descripción de los argumentos allí con lo que he hecho aquí.
• Tenga en cuenta la extern "C"sección en la parte superior, y que he agregado un guión bajo a dgeev_. Esto se debe a que la biblioteca fue escrita y construida en Fortran, por lo que esto es necesario para que los símbolos coincidan al vincular. Esto depende del compilador y del sistema, por lo que si usa esto en Windows, todo tendrá que cambiar.
• Algunas personas pueden sugerir el uso de la interfaz C para LAPACK . Puede que tengan razón, pero siempre lo he hecho de esta manera.

Normalmente me resisto a decirle a la gente lo que creo que deberían hacer en lugar de responder a su pregunta, pero en este caso voy a hacer una excepción.

Lapack está escrito en FORTRAN y la API es muy similar a FORTRAN. Hay una API C para Lapack que hace que la interfaz sea un poco menos dolorosa, pero nunca será una experiencia agradable usar Lapack desde C ++.

Alternativamente, hay una biblioteca de clases de matriz de C ++ llamada Eigen que tiene muchas de las capacidades de Lapack, proporciona un rendimiento computacional comparable a las mejores implementaciones de Lapack y es muy conveniente de usar desde C ++. En particular, así es como se podría escribir su código de ejemplo usando 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;
}

Este ejemplo de problema de valor propio es un caso de prueba para la función Lapack dgeev. Puede ver el código FORTRAN y los resultados de este problema dgeev ejemplo y hacer sus propias comparaciones.

Aquí hay otra respuesta en la misma línea que la anterior.

Debería buscar en la biblioteca de álgebra lineal Armadillo C ++ .

Pros:

1. La sintaxis de la función es de alto nivel (similar a la de MATLAB). Entonces no hay DGESVmumbo-jumbo, solo X = solve( A, B )(aunque hay una razón detrás de esos nombres de funciones LAPACK de aspecto extraño ...).
2. Implementa varias descomposiciones matriciales (LU, QR, valores propios, SVD, Cholesky, etc.)
3. Es rápido cuando se usa correctamente.
4. Está bien documentado .
5. Tiene soporte para matrices dispersas (deseará examinarlas más adelante).
6. Puede vincularlo con sus bibliotecas BLAS / LAPACK súper optimizadas para un rendimiento óptimo.

Así es como se vería el código de @ BillGreene con 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;
}