scalapack - Entendiendo las llamadas de LAPACK en C++ con un ejemplo simple
libreria blas (3)
Soy un principiante con la interfaz de LAPACK y C ++ / Fortran. Necesito resolver problemas de ecuaciones lineales y valores propios utilizando LAPACK / BLAS en Mac OS-X Lion. OS-X Lion proporciona bibliotecas optimizadas BLAS y LAPACK (en / usr / lib) y estoy vinculando estas bibliotecas en lugar de descargarlas desde netlib.
Mi programa (reproducido a continuación) está compilando y ejecutándose bien, pero me está dando respuestas incorrectas. He investigado en la web y Stackoverflow y el problema puede tener que ver con la forma en que C ++ y Fortran almacenan las matrices en diferentes formatos (fila mayor frente a columna principal). Sin embargo, como verá en mi ejemplo, la matriz simple para mi ejemplo debería verse idéntica en C ++ y fortran. De todos modos aquí va.
Digamos que queremos resolver el siguiente sistema lineal:
x + y = 2
x - y = 0
La solución es (x, y) = (1,1). Ahora traté de resolver esto usando Lapack como sigue
// 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);
}
Este código fue compilado de la siguiente manera:
g++ -Wall -llapack -lblas lapacktest.cpp
Sin embargo, al ejecutar esto, obtengo la solución como (-2,2) que obviamente es incorrecta. He intentado todas las combinaciones de reorganización de fila / columna de mi matriz a . También observe que la representación matricial de a debe ser idéntica en los formatos de fila y columna. Creo que hacer que este sencillo ejemplo funcione me ayudará a comenzar con LAPACK y cualquier ayuda será apreciada.
Para aquellos que no quieren molestarse con Accelerate Framework, les proporciono el código de Stephen Canon (gracias a él, por supuesto) con nada más que un enlace de biblioteca puro
// 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);
}
Y sobre el manual, hay una versión completa en PDF disponible en el sitio web de Intel. Aquí hay una muestra de su documentación HTML.
Si desea utilizar LAPACK de C ++, es posible que desee echar un vistazo a FLENS . Define las interfaces de bajo y alto nivel para LAPACK, pero también reimplementa algunas funciones de LAPACK.
Con la interfaz FLENS-LAPACK de bajo nivel, puede usar sus propios tipos de matriz / vector (si tienen un diseño de memoria compatible con LAPACK). Tu llamada de dgetrf se vería así:
info = lapack::getrf(NoTrans, dim, nrhs, a.begin(), LDA, ipiv);
y para los dgetrs
lapack::getrs(NoTrans, dim, nrhs, a.begin(), LDA, ipiv, b.begin(), LDB);
Por lo tanto, las funciones FLENS-LAPACK de bajo nivel están sobrecargadas con respecto a los tipos de elementos. En consecuencia, la función LAPACK sgetrs , dgetrs , cgetrs , zgetrs encuentran en la interfaz de bajo nivel de FLENS-LAPACK lapack::getrs . También pasa los parámetros por valor / referencia y no como puntero (por ejemplo, LDA lugar de &LDA ).
Si utiliza los tipos de matriz FLENS, puede codificarlo como
info = lapack::trf(NoTrans, A, ipiv);
if (info==0) {
lapack::trs(NoTrans, A, ipiv, b);
}
O simplemente utiliza la función dgesv controlador LAPACK
lapack::sv(NoTrans, A, ipiv, b);
Aquí una lista de las funciones del controlador FLENS-LAPACK .
Descargo de responsabilidad: Sí, FLENS es mi bebé! Eso significa que codifiqué alrededor del 95% y cada línea de código valió la pena.
dgetrf factorizar la matriz (llamando a dgetrf ) antes de poder resolver el sistema utilizando dgetrs . Alternativamente, puede usar la rutina dgesv , que realiza ambos pasos por usted.
Por cierto, no necesita declarar las interfaces usted mismo, están en los encabezados de Accelerate:
// 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);
}