Hola, estos días me he puesto a aprender en serio las librerías GSL junto con C++ en mi Ubuntu. Por ahí van rulando unos apuntes en inglés sobre dichas librerías, son muy extensos. Pero la verdad es que vale la pena utilizar estas librerías.
En dichos apuntes hay un ejemplo de cómo multiplicar dos matrices, pero lo hace para un caso particular; cosa que no me agrada, es muy poco matemático. Así que me he puesto a hacerlo para el caso general. La verdad es que me ha costado más de lo que creía por obcecarme a hacerlo como viene ahí. Después de pensar junto con la ayuda de mi amigo JamesMR (( y de preguntar en Ubuntu-Es )) , con el cual aprendimos a hacerlo de dos maneras diferentes poniendo punteros a las matrices y a los size_t; pues al final lo he conseguido como yo quería 
Os pongo el código del programa y explicaciones sobre ello.
#include<iostream>
#include<gsl/gsl_matrix.h>
#include<gsl/gsl_blas.h>
using namespace std;
int main(int argc, char **argv)
{
// Declaramos las variables que serán el número de filas y columnas de las matrices
size_t m,n,m1,n1;
//Declaramos los contadores
size_t i,j;
cout<<”\nIntroduce el número de filas de la matriz A “;
cin>>m;
cout<<”\nIntroduce el número de columnas de la matriz A “;
cin>>n;
cout<<”\nIntroduce el número de filas de la matriz B “;
cin>>m1;
cout<<”\nIntroduce el número de columnas de la matriz B “;
cin>>n1;
//Declaramos las matrices de forma dinámica con gsl, así podremos liberar memoria cuando acabemos
gsl_matrix * A = gsl_matrix_alloc(m,n);
gsl_matrix * B = gsl_matrix_alloc(m1,n1);
gsl_matrix * C = gsl_matrix_alloc(m,n1);
cout<<”\nIntroduce los elementos de la matriz A:”<<endl;
for(i=0; i<m*n; i++)
{
cout<<”A(“<<i+1<<”)”<<”=”;
cin>>A->data[i];
}
cout<<”\nIntroduce los elementos de la matriz B:”<<endl;
for(i=0; i<m1*n1; i++)
{
cout<<”B(“<<i+1<<”)”<<”=”;
cin>>B->data[i];
}
//Con esta orden multiplicamos A*B y guardamos el resultado en C
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, A, B, 0.0, C);
cout<<”\nEl resultado de multiplicar A por B es:”<<endl;
for(i=0; i<m; i++)
{
for(j=0;j<n1; j++)
{
cout<<gsl_matrix_get(C,i,j)<<”\t”;
}
cout<<endl;
}
cout<<endl;
// Liberamos Memoria
gsl_matrix_free(A);
gsl_matrix_free(B);
gsl_matrix_free(C);
return(0);
}
Bueno ese es el código. Hay que tener en cuenta un par de cosas. La primera es que GSL toma las matrices mediante la implementación de la struct gsl_matrix, la cual ya viene implementada en GSL con la forma:
typedef struct
{
size_t size1;
size_t size2;
size_tda tda;
double * data;
gsl_block * block;
int owner;
} gsl_matrix;
De esa struct lo más importante es que size1 es el número de filas, size2 es el número de columnas y data es un array que nos permite acceder a los elementos de la matriz de forma dinámica. De esa forma en el código que he introducido antes en donde pone:
for(i=0; i<m*n; i++)
{
cout<<”A(“<<i+1<<”)”<<”=”;
cin>>A->data[i];
}
En la línea cin>>A->data[i]; lo que hacemos es pedir por terminal los elementos de la matriz A con cin e introducirlos en data[m*n]. De igual forma lo hacemos con el segundo bucle for para la matriz B.
Bien, ahora veamos con más detenimiento la línea en la que ordenamos a nuestro programa que multiplique A por B y guarde el resultado en C, que es ésta:
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, A, B, 0.0, C);
El código general de la orden gsl_blas_dgemm es
int gsl_blas_dgemm(CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB, double alpha, const gsl_matrix * A, const gsl_matrix * B, double beta, gsl_matrix * C)
Con esta orden podemos hacer en general la operación C = alfa*op(A)*op(B)+beta*C, donde op(A) u op(B) pueden ser A, A transpuesta, A conjugada para TransA=CblasNoTrans, CblasTrans, CblasConjTrans. De igual forma para la op(B). Y alpha y beta son escalares de tipo real doble.
Por tanto esta orden nos permite multiplicar A por B si ponemos CblasNoTrans, alpha=1.0, beta=0.0, que es lo que hemos puesto en nuestro código. Pero también nos permite multiplicar la transpuesta de A por la transpuesta de B, o la transpuesta de A por B, etc.. Eso sí, sólo nos lo hará si el número de columnas de A coincide con el número de filas de B.
Si queréis probarlo podéis guardarlo con el nombre matprod.cpp y compilarlo con la orden:
g++ matprod.cpp -o matprod -lgsl -lgslcblas -lm
./matprod
La verdad es que con éste código tenemos un programa en el que podemos con una sóla línea multiplicar matrices (sin utilizar los for ), pero además, nos permite obtener mucha precisión en matrices con números reales; con lo cual con un poco de maña podemos hacer un algoritmo que nos resuelva sistemas de ecuaciones lineales numéricos (( con decimales )) obteniendo unos resultados bastante aceptables ( con una aceptable precisión ). Aunque en breve os mostraré cómo resolver dichos sistemas con los típicos algoritmos pivotales pero utilizando GSL.
Enlace Relacionado-> Instalar GSL en Ubuntu
Enlace Relacionado->Web De GSL
Saludos
