Comprendre les appels LAPACK en C ++ avec un exemple simple

Je suis un débutant avec LAPACK et C++/Fortran interfaçage. J'ai besoin pour résoudre des équations linéaires et les valeurs propres problèmes à l'aide de LAPACK/BLAS sur Mac OS X Lion. OS X Lion fournit optimisé BLAS et LAPACK (bibliothèques dans /usr/lib) et je suis reliant ces bibliothèques au lieu de les télécharger à partir de netlib.

Mon programme (reproduit ci-dessous) est la compilation et l'exécution de l'amende, mais il me donne de mauvaises réponses. J'ai fait des recherches sur le web et de Stackoverflow et le problème peut avoir à traiter que du C++ et Fortran stocker des tableaux dans différents formats (ligne principale vs Colonne principale). Cependant, comme vous le verrez dans mon exemple, le tableau simple pour mon exemple, devraient être identiques, en C++ et fortran. De toute façon, ici, va.

Disons que nous voulons résoudre le système linéaire suivant:

x + y = 2

x - y = 0

La solution est (x,y) = (1,1). Maintenant, j'ai essayé de résoudre ce à l'aide de Lapack comme suit

//LAPACK test code

#include<iostream>
#include<vector>

using namespace std;
extern "C" void dgetrs(char *TRANS, int *N, int *NRHS, double *A, 
                      int *LDA, int *IPIV, double *B, int *LDB, int *INFO );

int main()
{
    char trans = 'N';
    int dim = 2;    
    int nrhs = 1;
    int LDA = dim;
    int LDB = dim;
    int info;

    vector<double> a, b;

    a.push_back(1);
    a.push_back(1);
    a.push_back(1);
    a.push_back(-1);

    b.push_back(2);
    b.push_back(0);

    int ipiv[3];


    dgetrs(&trans, &dim, &nrhs, & *a.begin(), &LDA, ipiv, & *b.begin(), &LDB, &info);


    std::cout << "solution is:";    
    std::cout << "[" << b[0] << ", " << b[1] << ", " << "]" << std::endl;
    std::cout << "Info = " << info << std::endl; 

    return(0);
}

Ce code a été réalisé de la façon suivante:

g++ -Wall -llapack -lblas lapacktest.cpp

Sur l'exécution de cette, cependant, je reçois la solution (-2,2), qui est évidemment faux. J'ai essayé toute combinaison de ligne/colonne ré-arrangement de ma matrice a. Observez également la représentation matricielle de a doit être identique dans la rangée et de la colonne formats. Je pense que l'obtention de ce simple exemple de travail permettra de me commencé avec LAPACK et toute aide sera appréciée.

source d'informationauteur RDK

c++fortran lapack

Vous devez facteur la matrice (en appelant dgetrf) avant de pouvoir résoudre le système à l'aide de dgetrs. Alternativement, vous pouvez utiliser le dgesv routine, qui assure à la fois les démarches pour vous.

Par le façon, vous n'avez pas besoin de déclarer les interfaces de vous-même, ils sont dans l'Accélérer en-têtes:

//LAPACK test code

#include <iostream>
#include <vector>
#include <Accelerate/Accelerate.h>

using namespace std;

int main()
{
    char trans = 'N';
    int dim = 2;    
    int nrhs = 1;
    int LDA = dim;
    int LDB = dim;
    int info;

    vector<double> a, b;

    a.push_back(1);
    a.push_back(1);
    a.push_back(1);
    a.push_back(-1);

    b.push_back(2);
    b.push_back(0);

    int ipiv[3];

    dgetrf_(&dim, &dim, &*a.begin(), &LDA, ipiv, &info);
    dgetrs_(&trans, &dim, &nrhs, & *a.begin(), &LDA, ipiv, & *b.begin(), &LDB, &info);


    std::cout << "solution is:";    
    std::cout << "[" << b[0] << ", " << b[1] << ", " << "]" << std::endl;
    std::cout << "Info = " << info << std::endl; 

    return(0);
}

Pour ceux qui ne veulent pas s'embêter avec les Accélérer Cadre, je donne le code de Stephen Canon (merci à lui, bien sûr) avec rien, mais pure de la bibliothèque de liaison

//LAPACK test code
//compile with: g++ main.cpp -llapack -lblas -o testprog
#include <iostream>
#include <vector>
using namespace std;
extern "C" void dgetrf_(int* dim1, int* dim2, double* a, int* lda, int* ipiv, int* info);
extern "C" void dgetrs_(char *TRANS, int *N, int *NRHS, double *A, int *LDA, int *IPIV, double *B, int *LDB, int *INFO );
int main()
{
char trans = 'N';
int dim = 2;
int nrhs = 1;
int LDA = dim;
int LDB = dim;
int info;
vector<double> a, b;
a.push_back(1);
a.push_back(1);
a.push_back(1);
a.push_back(-1);
b.push_back(2);
b.push_back(0);
int ipiv[3];
dgetrf_(&dim, &dim, &*a.begin(), &LDA, ipiv, &info);
dgetrs_(&trans, &dim, &nrhs, & *a.begin(), &LDA, ipiv, & *b.begin(), &LDB, &info);
std::cout << "solution is:";
std::cout << "[" << b[0] << ", " << b[1] << ", " << "]" << std::endl;
std::cout << "Info = " << info << std::endl;
return(0);
}

Et sur le manuel, il ya une pleine version PDF disponible au site web d'Intel. Voici un échantillon de leur documentation HTML.

http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/GUID-A02DB70F-9704-42A4-9071-D409D783D911.htm

2

Si vous souhaitez utiliser LAPACK à partir de C++, vous pouvez avoir un look un FLENS. Il définit de bas et de haut niveau des interfaces de LAPACK mais aussi de ré-implémente certaines LAPACK fonctions.

Avec le faible niveau de FLENS-LAPACK interface, vous pouvez utiliser votre propre matrice/vecteur types (si ils ont un LAPACK conforme disposition de la mémoire). Votre appel de dgetrf ressembleraient à ça:
```
info = lapack::getrf(NoTrans, dim, nrhs, a.begin(), LDA, ipiv);
```
et pour dgetrs
```
lapack::getrs(NoTrans, dim, nrhs, a.begin(), LDA, ipiv, b.begin(), LDB);
```
De sorte que le faible niveau de FLENS-LAPACK fonctions sont surchargées à l'égard du type d'élément. Par conséquent, LAPACK fonction sgetrsdgetrscgetrszgetrs sont dans l'interface de bas niveau de FLENS-LAPACK lapack::getrs. Vous aussi passer des paramètres par valeur/référence et non en tant que pointeur (par exemple LDA au lieu de &LDA).

Si vous utilisez le FLENS la matrice des types que vous pouvez coder en tant que
```
info = lapack::trf(NoTrans, A, ipiv);
if (info==0) {
lapack::trs(NoTrans, A, ipiv, b);
}
```
Ou vous venez d'utiliser la LAPACK fonction de pilote dgesv
```
lapack::sv(NoTrans, A, ipiv, b);
```
Ici un liste de FLENS-LAPACK fonctions du pilote.

Avertissement: Oui, FLENS est mon bébé! Cela signifie que j'ai codé environ 95% de celui-ci et chaque ligne de code en valait la peine.

Vous devez vous connecter pour publier un commentaire.