Como começar a usar o LAPACK em c ++?

10

Eu sou novo na ciência da computação e já aprendi métodos básicos para integração, interpolação, métodos como RK4, Numerov etc. em c ++, mas recentemente meu professor me pediu para aprender a usar o LAPACK para resolver problemas relacionados a matrizes. Como, por exemplo, encontrar valores próprios de uma matriz complexa. Eu nunca usei bibliotecas de terceiros e quase sempre escrevo minhas próprias funções. Estou pesquisando há vários dias, mas não consigo encontrar nenhum guia para lapack para amadores. Todos eles estão escritos em palavras que não entendo e não sei por que o uso de funções já escritas deve ser tão complicado. Eles estão cheios de palavras como zgeev, dtrsv etc. e estou frustrado. Eu só quero codificar algo como este pseudo-código:

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

Não sei se estou sendo boba ou amadora. Mas, novamente, isso não deve ser tão difícil, deveria? Eu nem sei se devo usar o LAPACK ou LAPACK ++. (Eu escrevo códigos em c ++ e não tenho conhecimento de Python ou FORTRAN) e como instalá-los.

Alireza
fonte
Talvez este exemplo seria útil: matrixprogramming.com/files/code/LAPACK
nukeguy
Se você está apenas começando, talvez seja mais fácil usar uma biblioteca mais simples como ArrayFire github.com/arrayfire/arrayfire . Você pode chamá-lo diretamente do C ++ e as APIs são mais simples e acho que ele pode fazer todas as operações que o LAPACK faz.
218 Vikram
Neste outro post, um usuário propõe seu próprio wrapper FLENS, que possui uma sintaxe muito boa que pode facilitar sua introdução ao LAPACK.
Zythos
Chamar diretamente as funções do LAPACK é muito tedioso e propenso a erros. Existem vários wrappers C ++ fáceis de usar para o LAPACK, que fornecem um uso muito mais fácil, como o Armadillo . Para o caso de uso específico da decomposição complexa de eigen, consulte a função amigável eig_gen () , que embaixo envolve essa monstruosidade LAPACK, zheev (JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, INFO), e reformata os autovalores e autovetores obtidos em representações padrão.
Hbrerkere 14/05/19

Respostas:

18

Discordo de algumas das outras respostas e digo que acredito que descobrir como usar o LAPACK é importante no campo da computação científica.

No entanto, há uma grande curva de aprendizado no uso do LAPACK. Isso ocorre porque está escrito em um nível muito baixo. A desvantagem disso é que parece muito enigmático e não agradável aos sentidos. A vantagem disso é que a interface é inequívoca e basicamente nunca muda. Além disso, as implementações do LAPACK, como a Intel Math Kernel Library, são realmente rápidas.

Para meus próprios propósitos, tenho minhas próprias classes C ++ de nível superior, que envolvem as sub-rotinas LAPACK. Muitas bibliotecas científicas também usam o LAPACK por baixo. Às vezes, é mais fácil usá-los, mas, na minha opinião, há muito valor na compreensão da ferramenta abaixo. Para esse fim, forneci um pequeno exemplo de trabalho escrito em C ++ usando LAPACK para você começar. Isso funciona no Ubuntu, com oliblapack3 pacote instalado e outros pacotes necessários para a construção. Provavelmente pode ser usado na maioria das distribuições Linux, mas a instalação do LAPACK e a vinculação podem variar.

Aqui está o arquivo 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;
}

Isso pode ser construído usando a linha de comando

g++ -o test_lapack test_lapack.cpp -llapack

Isso produzirá um executável chamado test_lapack. Eu configurei isso para ler em um arquivo de entrada de texto. Aqui está um arquivo chamado matrix.txtcontendo uma matriz 3x3.

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

Para executar o programa, basta digitar

./test_lapack matrix.txt

na linha de comando, e a saída deve ser

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

Comentários:

  • Você parece impressionado com o esquema de nomeação para LAPACK. Uma breve descrição está aqui .
  • A interface para a sub-rotina DGEEV está aqui . Você deve poder comparar a descrição dos argumentos lá com o que eu fiz aqui.
  • Note o extern "C" seção na parte superior e que adicionei um sublinhado adgeev_ . Isso ocorre porque a biblioteca foi escrita e construída no Fortran, portanto, é necessário fazer com que os símbolos correspondam ao vincular. Isso depende do compilador e do sistema; portanto, se você usá-lo no Windows, tudo terá que mudar.
  • Algumas pessoas podem sugerir o uso da interface C para o LAPACK . Eles podem estar certos, mas eu sempre fiz dessa maneira.
LedHead
fonte
3
Muito do que você está procurando pode ser encontrado com um Googlage rápido. Talvez você não tenha certeza do que procurar. O Netlib é o goleiro do LAPACK. A documentação pode ser encontrada aqui . Esta página possui uma tabela útil das principais funcionalidades do LAPACK. Alguns dos mais importantes são (1) resolver sistemas de equações, (2) problemas de autovalores, (3) decomposições de valores singulares e (4) fatorações QR. Você entendeu o manual para DGEEV?
LedHead
11
Eles são apenas interfaces diferentes para a mesma coisa. LAPACK é o original. Está escrito em Fortran, portanto, para usá-lo, você precisa jogar alguns jogos para fazer a compilação cruzada do C / C ++ funcionar, como mostrei. Eu nunca usei o LAPACKE, mas parece que é um invólucro em C bastante fino sobre o LAPACK que evita esse negócio de compilação cruzada, mas ainda é de baixo nível. O LAPACK ++ parece ser um invólucro C ++ de nível ainda mais alto, mas acho que nem é mais suportado (alguém me corrija se eu estiver errado).
LedHead
11
Não conheço nenhuma coleção de códigos específica. Mas se você pesquisar no Google algum dos nomes de sub-rotina LAPACK, invariavelmente encontrará uma pergunta antiga em um dos sites StackExchange.
LedHead
11
@AlirezaHashemi A propósito, o motivo pelo qual você deve fornecer a matriz WORK é porque, em regra, o LAPACK não aloca memória nas sub-rotinas. Se estamos usando o LAPACK, provavelmente estamos usando grandes quantidades de memória e a alocação de memória é cara, por isso faz sentido permitir que as rotinas de chamada sejam responsáveis ​​pela alocação de memória. Como o DGEEV requer memória para armazenar quantidades intermediárias, precisamos fornecer esse espaço de trabalho para ele.
LedHead
11
Entendi. E eu escrevi com sucesso meu primeiro código para calcular autovalores de uma matriz complexa usando zgeev. E já está fazendo mais! Obrigado!
Alireza
7

Normalmente, resisto a dizer às pessoas o que acho que elas deveriam fazer, em vez de responder à pergunta, mas, neste caso, vou abrir uma exceção.

Lapack é escrito em FORTRAN e a API é muito parecida com FORTRAN. Existe uma API C para o Lapack que torna a interface um pouco menos dolorosa, mas nunca será uma experiência agradável usar o Lapack do C ++.

Como alternativa, existe uma biblioteca de classes de matriz C ++ chamada Eigen que possui muitos dos recursos do Lapack, fornece desempenho computacional comparável às melhores implementações do Lapack e é muito conveniente de usar no C ++. Em particular, eis como seu código de exemplo pode ser escrito 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 exemplo de problema de autovalor é um caso de teste para a função Lapack dgeev. Você pode visualizar o código FORTRAN e os resultados desse exemplo de problema dgeev e fazer suas próprias comparações.

Bill Greene
fonte
Obrigado pela sua resposta e explicação! Vou tentar esta biblioteca e escolher a que melhor se adequa às minhas necessidades.
Alireza
Oh, eles sobrecarregam operator,! Nunca vi isso feito na prática :-)
Wolfgang Bangerth
11
Na verdade, essa operator,sobrecarga é mais interessante / melhor do que parece à primeira vista. É usado para inicializar matrizes. As entradas que inicializam a matriz podem ser constantes escalares, mas também podem ser matrizes ou sub-matrizes definidas anteriormente. Muito parecido com o MATLAB. Desejo a minha capacidade C ++ programação era bom o suficiente para implementar algo que sofisticado me ;-)
Bill Greene
7

Aqui está outra resposta na mesma linha que a anterior.

Você deve procurar na biblioteca de álgebra linear do Armadillo C ++ .

Prós:

  1. A sintaxe da função é de alto nível (semelhante à do MATLAB). Portanto, nada de DGESVbobagem, apenas X = solve( A, B )(embora haja uma razão por trás desses nomes de função LAPACK de aparência estranha ...).
  2. Implementa várias decomposições de matriz (LU, QR, autovalores, SVD, Cholesky, etc.)
  3. É rápido quando usado corretamente.
  4. Está bem documentado .
  5. Tem suporte para matrizes esparsas (você precisará examiná-las posteriormente).
  6. Você pode vinculá-lo às suas bibliotecas BLAS / LAPACK super otimizadas para obter o desempenho ideal.

Veja como o código do @ BillGreene ficaria com o Tatu:

#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;
}
GoHokies
fonte
Obrigado pela sua resposta e explicação! Vou tentar esta biblioteca e escolher a que melhor se adequa às minhas necessidades.
Alireza