Mezclando BOOST, uBLAS y CLAPACK para resolver sistemas de ecuaciones
Si bien anteriormente me entusiasmaron Newmat y MTL4, ninguna de ellas era capaz de resolver lo que necesitaba. La primera no resuelve sistemas y la segunda tardaba muchsímimo y hasta petaba en grandes sistemas (22k ecuaciones). Pero seguí investigando para utilizar las famosas librerías de resolución LAPACK o suery su versión para C: CLAPACK. Sin embargo, esta biblioteca no permite el almacenamiento eficiente de matrices dispersas/simétricas y para eso existe otra biblioteca llamada BLAS y, concretamente, la que viene ya incluida en BOOST: uBLAS.
Este post relata cómo, utilizando el almacenamiento de matrices de uBLAS y un binding sobre BOOST de CLAPACK y el propio CLAPACK, se resuelven sistemas de ecuaciones extremadamente grandes en poco tiempo.
Lo primero es tener instalado y preparado BOOST, con esto ya tenemos uBLAS a punto. Seguidamente podemos probar que esté correcto con un programita simple como este:
#include "stdafx.h"
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/io.hpp>
int main(){
using namespace boost::numeric::ublas;
mapped_matrix<double> m (5, 5, 3 * 3);
for (unsigned i = 0; i < m.size1 (); ++ i)
for (unsigned j = 0; j < m.size2 (); ++ j)
m (i, j) = 5 * i + j*5.2;
m(3,3)=10;
std::cout << m << std::endl;
matrix<double> inv(5,5);
permutation_matrix<std::size_t> pm(m.size1());
lu_factorize(m,pm);
inv.assign(identity_matrix<double>(m.size1()));
lu_substitute(m,pm,inv);
std::cout << inv << std::endl;
system("PAUSE");
return 0;
}
Lo que hace es rellenar una matriz y luego invertirla. Es la forma más rápida de invertir una matriz utilizando uBLAS que he encontrado.
Matrices dispersas
Existen tres maneras de guardar matrices dispersas con uBLAS, yo he utilizado la compressed_matrix, me pareció la más eficiente. Aquí un ejemplo:
boost::numeric::ublas::compressed_matrix<float> miMatrizDispersa(ancho,alto);
El acceso para lectura y escritura es con miMatrizDispersa(i,j), todo muy fácil, rápido y sencillo.
Añadir el binding de LAPACK
Por sí misma, la biblioteca uBLAS no resuelve sistemas, podría utilizarse para invertir una matriz y multiplicarla por un vector, resolviendo así el sistema, pero esto es extremadamente ineficiente y es muy buena idea utilizar CLAPACK. Pero para ello tenemos que añadir el binding uBLAS-LAPACK. Esto es muy sencillo:
- Colocarse en el directorio donde tengamos BOOST, típicamente: C:/Archivos de programa/boost/boost_1_44/boost/numeric/
- Crear la carpeta bindings/
- Utilizar un cliente de SubVersion para bajar todo el árbol de:
svn co http://svn.boost.org/svn/boost/sandbox/numeric_bindings-v1/boost/numeric/bindings
- Colocarse en el directorio: C:/Archivos de programa/boost/boost_1_44/libs/numeric/
- Crear la carpeta bindings/
- Utilizar un cliente de SubVersion para bajar todo el árbol de:
svn co http://svn.boost.org/svn/boost/sandbox/numeric_bindings-v1/libs/numeric/bindings
Tendremos ya los bindings instalados y fácilmente comprobable. Pero antes…
Instalar y vincular CLAPACK
Debemos bajar CLAPACK para Visual Studio y colocarlo donde queramos, por ejemplo en C:/CLAPACK-3.1.1-VisualStudio. Importante bajarse las cosas ya compiladas a ser posible. Seguidamente añadimos en Visual Studio los directorios include y lib/Win32 como directorios adicionales de inclusión y biblioteca respectivamente.
Luego en el proyecto tenemos que añadir las bibliotecas adicionales: clapackd.lib y BLASd.lib en el vinculador.
Resolver un sistema de ecuaciones
Nos falta una cosa importante para que toda esta mezcla funcione, hay que añadir lo siguiente al principio del fichero donde se incluyan las cosas importantes:
#define BOOST_NUMERIC_BINDINGS_USE_CLAPACK //¡muy importante para que todo funcione! #include <boost/numeric/bindings/traits/ublas_matrix.hpp> #include <boost/numeric/bindings/lapack/gesv.hpp> #include <boost/numeric/ublas/matrix_sparse.hpp> #include <boost/numeric/ublas/vector.hpp> #include <boost/numeric/ublas/io.hpp>
Con este define y los includes que necesitemos ya podemos rellenar la matriz y resolver así:
ublas::matrix<float, ublas::column_major> A(ancho,alto); //importante el column_major para compatibilidad con LAPACK ublas::matrix<float, ublas::column_major> b(alto,1); //ponemos una matriz alto x 1 para hacerlo más sencillo, pero seguramente con un vector //rellenar A y b [...] bindings::lapack::gesv(A,b); //la solución se encontrará en b
Así hemos conseguido finalmente resolver el sistema Ax=b y sorprenderá la velocidad con que lo hace CLAPACK, el cuello de botella de aquí es la inserción de datos en una matriz dispersa.
