Bagaimana cara mulai menggunakan LAPACK di c ++?

10

Saya baru mengenal ilmu komputasi dan saya sudah mempelajari metode dasar untuk integrasi, interpolasi, metode seperti RK4, Numerov dll pada c ++ tetapi baru-baru ini profesor saya meminta saya untuk belajar bagaimana menggunakan LAPACK untuk menyelesaikan masalah yang berkaitan dengan matriks. Seperti misalnya menemukan nilai eigen dari matriks kompleks. Saya tidak pernah menggunakan perpustakaan pihak ketiga dan saya hampir selalu menulis fungsi saya sendiri. Saya telah mencari-cari selama beberapa hari tetapi tidak dapat menemukan panduan ramah-amatir untuk Lapack. Semuanya ditulis dengan kata-kata yang saya tidak mengerti dan saya tidak tahu mengapa menggunakan fungsi yang sudah ditulis harus rumit ini. Mereka penuh dengan kata-kata seperti zgeev, dtrsv, dll. Dan saya frustrasi. Saya hanya ingin kode seperti pseudo-code ini:

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

Saya tidak tahu apakah saya bodoh atau amatir. Tetapi sekali lagi, ini seharusnya tidak terlalu sulit bukan? Saya bahkan tidak tahu apakah saya harus menggunakan LAPACK atau LAPACK ++. (Saya menulis kode dalam c ++ dan tidak memiliki pengetahuan tentang Python atau FORTRAN) dan cara menginstalnya.

Alireza
sumber
Mungkin contoh ini akan berguna: matrixprogramming.com/files/code/LAPACK
nukeguy
Jika Anda baru memulai, mungkin akan lebih mudah menggunakan perpustakaan yang lebih sederhana seperti ArrayFire github.com/arrayfire/arrayfire . Anda dapat langsung menyebutnya dari C ++ dan API lebih sederhana dan saya pikir itu bisa melakukan semua operasi yang dilakukan LAPACK.
Vikram
Dalam postingan lain ini pengguna mengusulkan FLENS pembungkusnya sendiri, yang memiliki sintaks yang sangat bagus yang dapat memudahkan pengenalan Anda ke LAPACK.
Zythos
Memanggil fungsi LAPACK secara langsung sangat membosankan dan rawan kesalahan. Ada beberapa pembungkus C ++ yang mudah digunakan untuk LAPACK yang memberikan penggunaan lebih mudah, seperti Armadillo . Untuk kasus penggunaan spesifik dari dekomposisi eigen kompleks, lihat fungsi eig_gen () yang mudah digunakan, yang di bawahnya membungkus kebodohan, zheev (JOBZ, UPLO, N, A, LDA, W, KERJA, LWORK, RWORK, INFO), dan memformat ulang nilai eigen yang diperoleh dan vektor eigen menjadi representasi standar.
hbrerkere

Jawaban:

18

Saya akan tidak setuju dengan beberapa jawaban lain dan mengatakan bahwa saya percaya bahwa mencari tahu cara menggunakan LAPACK adalah penting dalam bidang komputasi ilmiah.

Namun, ada kurva belajar yang besar untuk menggunakan LAPACK. Ini karena ini ditulis pada tingkat yang sangat rendah. Kerugian dari itu adalah bahwa itu tampak sangat samar, dan tidak menyenangkan indra. Keuntungannya adalah bahwa antarmuka tidak ambigu dan pada dasarnya tidak pernah berubah. Selain itu, implementasi LAPACK, seperti Perpustakaan Intel Math Kernel sangat cepat.

Untuk keperluan saya sendiri, saya memiliki kelas C ++ tingkat tinggi saya sendiri yang membungkus subrutin LAPACK. Banyak perpustakaan ilmiah juga menggunakan LAPACK di bawahnya. Terkadang lebih mudah untuk hanya menggunakannya, tetapi menurut saya ada banyak nilai dalam memahami alat di bawahnya. Untuk itu, saya telah memberikan contoh kerja kecil yang ditulis dalam C ++ menggunakan LAPACK untuk membantu Anda memulai. Ini berfungsi di Ubuntu, dengan liblapack3paket diinstal, dan paket lain yang diperlukan untuk membangun. Ini mungkin dapat digunakan di sebagian besar distribusi Linux, tetapi instalasi LAPACK dan menghubungkannya dapat bervariasi.

Ini filenya 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;
}

Ini dapat dibangun menggunakan baris perintah

g++ -o test_lapack test_lapack.cpp -llapack

Ini akan menghasilkan nama yang dapat dieksekusi test_lapack. Saya telah mengatur ini untuk dibaca dalam file input teks. Berikut adalah file bernama yang matrix.txtmengandung matriks 3x3.

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

Untuk menjalankan program cukup ketik

./test_lapack matrix.txt

di baris perintah, dan output seharusnya

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

Komentar:

  • Anda sepertinya terlempar oleh skema penamaan untuk LAPACK. Deskripsi singkat ada di sini .
  • Antarmuka untuk subroutine DGEEV ada di sini . Anda harus dapat membandingkan deskripsi argumen di sana dengan apa yang saya lakukan di sini.
  • Perhatikan extern "C"bagian di atas, dan saya telah menambahkan garis bawah ke dgeev_. Itu karena perpustakaan ditulis dan dibangun di Fortran, jadi ini perlu untuk membuat simbol cocok ketika menghubungkan. Ini adalah kompiler dan bergantung pada sistem, jadi jika Anda menggunakan ini pada Windows, semuanya harus berubah.
  • Beberapa orang mungkin menyarankan menggunakan antarmuka C ke LAPACK . Mereka mungkin benar, tetapi saya selalu melakukannya dengan cara ini.
LedHead
sumber
3
Banyak hal yang Anda cari dapat ditemukan dengan beberapa Googlage cepat. Mungkin Anda tidak yakin apa yang harus dicari. Netlib adalah penjaga LAPACK. Dokumentasi dapat ditemukan di sini . Halaman ini memiliki tabel praktis tentang fungsi utama LAPACK. Beberapa yang penting adalah (1) menyelesaikan sistem persamaan, (2) masalah nilai eigen, (3) dekomposisi nilai singular, dan (4) faktorisasi QR. Apakah Anda memahami manual untuk DGEEV?
LedHead
1
Mereka semua hanyalah antarmuka yang berbeda untuk hal yang sama. LAPACK adalah aslinya. Ini ditulis dalam Fortran, jadi untuk menggunakannya Anda harus memainkan beberapa game untuk membuat kompilasi silang dari C / C ++ berfungsi, seperti yang saya tunjukkan. Saya tidak pernah menggunakan LAPACKE, tetapi sepertinya itu adalah bungkus C yang cukup tipis di atas LAPACK yang menghindari bisnis kompilasi silang ini, tetapi masih tingkat yang cukup rendah. LAPACK ++ tampaknya menjadi pembungkus C ++ tingkat yang lebih tinggi, tapi saya rasa itu tidak didukung lagi (seseorang mengoreksi saya jika saya salah).
LedHead
1
Saya tidak tahu pengumpulan kode tertentu. Tetapi jika Anda Google salah satu nama subrutin LAPACK, Anda akan selalu menemukan pertanyaan lama di salah satu situs StackExchange.
LedHead
1
@AlirezaHashemi Ngomong-ngomong, alasan Anda harus menyediakan array WORK adalah karena LAPACK tidak mengalokasikan memori apa pun di dalam subrutinnya. Jika kita menggunakan LAPACK, kemungkinan besar kita menggunakan memori yang banyak, dan mengalokasikan memori itu mahal, jadi masuk akal untuk membiarkan rutinitas panggilan bertanggung jawab atas alokasi memori. Karena DGEEV membutuhkan memori untuk menyimpan jumlah setengah jadi, kami harus menyediakan ruang kerja untuk itu.
LedHead
1
Oke. Dan saya berhasil menulis kode pertama saya untuk menghitung nilai eigen dari matriks kompleks menggunakan zgeev. Dan sudah melakukan lebih banyak! Terima kasih!
Alireza
7

Saya biasanya menolak memberi tahu orang apa yang menurut saya harus mereka lakukan daripada menjawab pertanyaan mereka, tetapi dalam kasus ini saya akan membuat pengecualian.

Lapack ditulis dalam FORTRAN dan API sangat mirip FORTRAN. Ada API C untuk Lapack yang membuat antarmuka sedikit lebih menyakitkan tetapi tidak akan pernah menjadi pengalaman yang menyenangkan untuk menggunakan Lapack dari C ++.

Atau, ada perpustakaan kelas matriks C ++ yang disebut Eigen yang memiliki banyak kemampuan Lapack, memberikan kinerja komputasi yang sebanding dengan implementasi Lapack yang lebih baik, dan sangat mudah digunakan dari C ++. Secara khusus, berikut adalah bagaimana kode contoh Anda dapat ditulis menggunakan 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;
}

Contoh masalah nilai eigen ini adalah kasus uji untuk fungsi Lapack dgeev. Anda dapat melihat kode FORTRAN dan hasil untuk masalah ini contoh dgeev dan membuat perbandingan Anda sendiri.

Bill Greene
sumber
Terima kasih atas jawaban dan penjelasan Anda! Saya akan mencoba perpustakaan ini dan memilih yang paling sesuai dengan kebutuhan saya.
Alireza
Oh, mereka membebani operator,! Belum pernah melihat itu dilakukan dalam praktik sebenarnya :-)
Wolfgang Bangerth
1
Sebenarnya, operator,kelebihan itu lebih menarik / lebih baik daripada yang pertama kali muncul. Ini digunakan untuk menginisialisasi matriks. Entri yang menginisialisasi matriks dapat berupa konstanta skalar tetapi juga dapat merupakan matriks atau sub-matriks yang telah ditentukan sebelumnya. Sangat mirip MATLAB. Semoga kemampuan pemrograman C ++ saya cukup baik untuk mengimplementasikan sesuatu yang canggih sendiri ;-)
Bill Greene
7

Berikut jawaban lain dengan nada yang sama seperti di atas.

Anda harus melihat ke perpustakaan aljabar linier Armadillo C ++ .

Pro:

  1. Sintaks fungsi adalah tingkat tinggi (mirip dengan MATLAB). Jadi tidak ada DGESVomong kosong, hanya X = solve( A, B )(meskipun ada alasan di balik nama-nama fungsi LAPACK yang tampak aneh ...).
  2. Menerapkan berbagai dekomposisi matriks (LU, QR, nilai eigen, SVD, Cholesky, dll.)
  3. Hal ini cepat bila digunakan dengan benar.
  4. Itu didokumentasikan dengan baik .
  5. Memiliki dukungan untuk matriks yang jarang (Anda akan ingin melihatnya nanti).
  6. Anda dapat menautkannya dengan perpustakaan BLAS / LAPACK super-optimal Anda untuk kinerja yang optimal.

Begini tampilan kode @ BillGreene dengan 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;
}
GoHokies
sumber
Terima kasih atas jawaban dan penjelasan Anda! Saya akan mencoba perpustakaan ini dan memilih yang paling sesuai dengan kebutuhan saya.
Alireza