Cómo aprender cuda en un mes…

Si bien hablé en otro post sobre la biblioteca Newmat, he seguido buscando alternativas pues Newmat no maneja matrices dispersas. He encontrado tres bibliotecas interesantes: SparseLib++, CSparse y MTL4.

De ellas decidí probar MTL4 para guardar una matriz de 1.400 millones de elementos y resolver un sistema también muy grande. He tenido bastantes problemillas al usarlo con Visual Studio 2005, pero al final parece que funciona.

Primeramente, MTL2 lo he descartado por obsoleto y he tirado directamente a MTL4, bajando todo el código del servidor subversion que proporcionan. Las instrucciones están muy bien detalladas y básicamente vienen a decir:

• Descarga e instala Boost C++
• Descarga MTL4 del subversion (con TortoiseSVN para Windows)
• Añadir los directorios de inclusión adicionales en Visual Studio:
  • El de boost, típicamente: C:\Archivos de programa\boost\boost_1_44
  • El de MTL4, típicamente: C:\mtl4
• Incluir las librerías necesarias cuando vayamos a usar MTL

#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/itl/itl.hpp>

Y ahora dejo unas muestras de código para hacer matrices, invertirlas…

Matrices densas

dense2D<float>   A(ancho, alto);
A[0][4]=55; //ponemos a 55 el valor de la posición 0,4
float v=A[0][4]; //sacamos el valor 0,4 para guardarlo en v
std::cout << "A is \n" << A; //imprimir una matriz es muy fácil

Matrices dispersas

Éstas son algo distintas en su manejo, no se pueden rellenar de la misma manera, necesitan un inserter.

compressed2D<double>   B(ancho, alto);
if(true){
  matrix::inserter<compressed2D<double>> ins(B); //crear un inserter sobre B
  ins[0][4] << 55;  //ponemos a 55 el valor de la posición 0,4
}
float v=B[0][4]; //la extracción es igual
std::cout << "B is \n" << B;

Nótese el if(true), eso lo he hecho porque la creación del inserter debe ser dentro de un bloque. Esto es así porque hasta que el inserter no se destruya, la matriz no podrá ser accesible y dará una excepción al intentar hacer B[0][4] o cualquier otro acceso. Se puede meter dentro de un bucle o función para que el inserter se destruya una vez hecho su trabajo.

Resolver un sistema

MTL sólo permite resolver sistemas con matrices densas, una pena. Se deben seguir estos pasos:

dense_vector<float> solucion(tam);
dense2D<float> A(tam,tam);
dense_vector<float> b(tam);

solucion=lu_solve(A,b);

Obviamente, teniendo en cuenta que todos los tamaños permitan la resolución del sistema.