[C] Multiplicación de matrices sparse por vectores

[Regexp]

Buenas chicos, les traigo un laboratorio de C, que hice en un curso de cómputo científico, espero les sea de utilidad a alguno.

Matriz sparse: es una matriz cuyos elementos no nulos figuran una minoría entre todos. Por ejemplo, una matriz 1000x1000 tiene 10⁶ elementos, supongamos que de los cuales sólo el 10% son no-nulos (10⁵). A la hora de resolver diversos problemas se hace útil almacenarla de una manera especial teniendo en cuenta esto.

(Pueden encontrar un poco más de información del tema en la wikipedia en inglés, también hay papers en google académico)

El programa que se presenta pide un archivo el cual tiene la siguiente estructura

indices de fila de elementos no nulos
indices de columna de elementos no nulos
valores no nulos de la matriz asociados al conjunto de indices
valores de un vector b
.... (así sucesivamente)

Se pide un número de matrices a leer. La función leeMatriz lee una matriz específica en el archivo (1 si es la primera, 2 si es la segunda, etc). Se supone que cada elemento de la matriz tiene al menos un elemento no nulo en cada fila y columna. La función matrizPorVector devuelve el puntero hacia el comienzo del arreglo que contiene los valores del resultado de la multiplicación.

Sin más, acá el código

main.c

Código (c) [Seleccionar] Expandir

// Programa principal
//  Se encarga de establecer interfaz con el usuario.

#include <stdio.h>
#include <stdlib.h>

//Prototipo
int leeMatriz(
char *nomarch, int nummat,
int **ifil, int **icol, double **xval, double **b,
int *nfil, int *ncol);

double *matrizPorVector(int *ifil, int *icol, double *xval, double *b, int nza, int nfil);

int main()
{
char nombre[100];
char salida[] = "resultados.txt";
FILE *archivo;
int num;
double *xval;
double *b;
int *ifil, *icol;
int nfil, ncol;
int retorno, i, k;
double *resultado;

//Se pide al usuario que introduzca los datos para la lectura
printf("Introduce nombre de archivo: ");
scanf("%s", nombre);

printf("Introduce numero de matrices en el archivo: ");
scanf("%d", &num);

//Abrimos el archivo de salida
if ( (archivo = fopen(salida, "w+")) == NULL)
{
printf("No se ha podido abrir el archivo de salida %s \n", salida);
exit(EXIT_FAILURE);
}

//Leemos la cantidad de matrices que el usuario ingresó
for(k=1; k<=num; k++)
{
//Leemos la k-ésima matriz
retorno = leeMatriz(nombre, k, &ifil, &icol, &xval, &b, &nfil, &ncol);

//Damos un mensaje alusivo
fprintf(archivo, "La matriz numero %d\r\n", k);

//Imprimimos en el archivo de salida, dependiendo de lo que retornó leeMatriz
switch(retorno)
{

case -2:
fprintf(archivo, "Los tamaños de los vectores no son compatibles\r\n");
break;

case -1:
fprintf(archivo, "El tamaño del vector no coincide con el numero de columnas\r\n");
break;

default:
fprintf(archivo, "Los tamaños coinciden. El vector resultado es:\r\n");

//Procesamos los datos para obtener el resultado
resultado = matrizPorVector(ifil, icol, xval, b, retorno, nfil);
for(i=0; i<nfil; i++)
fprintf(archivo, "%lf   ", *(resultado+i) );

//Liberamos la memoria del vector resultado solo en este caso
free(resultado);
fprintf(archivo, "\r\n");
break;
}
//Liberamos la memoria reservada para volver a llamar
//de manera segura la función malloc dentro de las funciones

//Liberamos la memoria de los vectores que generan el par matriz-vector
//luego de esto, quedan apuntando a NULL, tal cosa es lo que espera
//la funcion leeMatriz, al final del programa esto se hace de último
//y por lo tanto nos aseguramos de que liberamos toda la memoria dinamica
//que se utilizó a lo largo del programa
free(ifil);
free(icol);
free(xval);
free(b);

fprintf(archivo, "\r\n");
}

//Cerramos el archivo de salida
fclose(archivo);


return 0;
}

matrizPorVector.c

Código (c) [Seleccionar] Expandir

// Funcion que realiza la multiplicacion de una matriz sparse por un vector
// Argumentos.
// Entrada:
// int *ifil, vector de indices de filas
//  int *icol, vector de indices de columnas
//  int nza, cantidad de valores no nulos
//  int nfil, orden de filas
// double *xval, *b, los valores no nulos y el vector

// Salida (por return)
// double *, una dirección de memoria, el cual será el comienzo
//  del vector resultado de la multiplicación

#include <stdio.h>
#include <stdlib.h>

double *matrizPorVector(int *ifil, int *icol, double *xval, double *b, int nza, int nfil)
{
unsigned k; // unsigned i;
double *resultado;

if ( ( resultado = (double *) malloc ( nfil*( sizeof(double) ) ) ) == NULL )
{
printf("No hubo memoria disponible para la reserva de un resultado\n");
exit(EXIT_FAILURE);
}

//Inicializamos los valores del vector resultado
for(k=0; k<nfil; k++)
*(resultado + k) = 0.0f;

// Cada k tal que 0<=k<nza cumple que
// El elemento en la posicion *(ifil+k), *(icol+k) es *(xval+k)
// Sólo hay que recordar que el conjunto de indices en los
// vectores ifil, icol comienzan desde 1, no desde cero
// por lo tanto debemos restarle 1 cuando utilizamos este
// conjunto de indices

// Como resultado obtenemos un for más compacto, y además
// un algoritmo mucho más rapido para matrices muy grandes
for(k=0; k<nza; k++)
*(resultado+(*(ifil+k))-1) += (*(b - 1 + (*(icol+k))))*(*(xval+k));


// Retornamos la direccion de memoria del comienzo del arreglo
return resultado;
}


leeMatriz.c

Código (c) [Seleccionar] Expandir

// Funcion que lee matriz desde archivo de texto.
// Argumentos.
// Entrada:
// char *nomarch : nombre de archivo
// int nummat : numero de matriz a leer

// Salida (por referencia)
// int *ifil : vector de indices de filas
// int *icol : vector indices de columnas
// double *xval : vector de valores asociado a la matriz en los indices de ifil e icol
// double *b : vector asociado al par nummat de matrices-vectores
// int *nfil : numero de filas
// int *ncol : numero de columnas

// Salida (por return)
// -2, si los vectores xval, ifil e icol no tienen el mismo tamaño
// -1, el tamaño de b no corresponde con el nuero de columnas de la matriz
// a, si todos los valores se corresponden, donde a es el numero de elementos no nulos de la matriz


#include <stdio.h>
#include <stdlib.h>

int leeMatriz(
char *nomarch, int nummat,
int **ifil, int **icol, double **xval, double **b,
int *nfil, int *ncol)

{
int a = 0; // Valor de retorno
int nbytes; // Identificador para el retorno de fscanf */
int nval, nb; // Numero de valores no nulos de la matriz, y numero de valores del vector */
char c;
// Caracter que se lee después de los valores */
int foo, aux;
// Un contador y una variable auxiliars */
int bar[4]; // Variables logicas auxiliares

FILE *archivo;

// Validamos la correcta apertura del archivo
if ( (archivo = fopen(nomarch, "r")) == NULL)
{
printf("No se ha podido abrir el archivo %s \n", nomarch);
exit(EXIT_FAILURE);
}




// Saltamos las lineas necesarias hasta la matriz que se nos pide
for (foo = 1; foo<=4*(nummat-1); foo++)
{
do
{
fscanf(archivo, "%c" , &c );
} while( c != '\n' && c != '\r' );
}

//Inicializamos la salida, suponemos que los vectores apuntan a nulo en esta llamada.
//Los arreglos (leidos en ese orden)

//Guardamos como valor logico cada uno de los siguientes.
bar[0] = ( *ifil = (int *) malloc( sizeof(int) ) ) == NULL;
bar[1] = ( *icol = (int *) malloc( sizeof(int) ) ) == NULL;
bar[2] = ( *xval = (double *) malloc( sizeof(double) ) ) == NULL;
bar[3] = ( *b = (double *) malloc( sizeof(double) ) ) == NULL;

//Si alguno de los malloc retorna NULL, imprimimos un mensaje alusivo y salimos
if(bar[0] || bar[1] || bar[2] || bar[3])
{
printf("No hubo memoria disponible al inicializar los pares de vectores asociados a la matriz\n");
exit(EXIT_FAILURE);
}
//El contador de ifil, icol, y del arreglo b


//Nota importante: los punteros son pasados por referencia, o sea punteros a punteros de tipo
//Las lineas de arriba dicen el contenido del puntero a puntero (que es un puntero) es la direccion
//de memoria que retorna la función malloc en el tamaño especificado

*nfil = *ncol = nval = nb = 0;

//Leemos la entrada que se nos pide. La lectura completa consta de 4 líneas
//Con este for lo que pretendo es ahorrarnos líneas repetitivas


for(foo = 1; foo<=4; foo++)
{
do
{
switch(foo)
{

//Lectura de ifil

case 1:
fscanf(archivo, "%d" , (*ifil+(*nfil)) );
*ifil = (int *) realloc( *ifil, (++(*nfil)+1)*sizeof(int) );
break;


//Lectura de icol
case 2:
fscanf(archivo, "%d" , (*icol+(*ncol)) );
*icol = (int *) realloc( *icol, (++(*ncol)+1)*sizeof(int) );
break;

//Lectura de xval
case 3:
fscanf(archivo, "%lf" , *xval+nval);
*xval = (double *) realloc( *xval, (++nval+1)*sizeof(double) );
break;

//Lectura de b
case 4:
fscanf(archivo, "%lf" , *b+nb);
*b = (double *) realloc( *b, (++nb+1)*sizeof(double) );
break;


}

} while ( (nbytes = fscanf(archivo, "%c", &c)) != EOF && c != '\n' && c != '\r' );
//Cada iteración del for lee una línea
}

fclose(archivo);

//Dado que, cuando llega al fin de línea igualmente se le agrega
//un espacio a los vectores, quitamos este espacio extra

//Damos esta solución para el problema de validar que se ha hecho la reserva
//de manera correcta, en lugar de uno por uno, se haen todos, si alguno
//no ha salido como se esperaba, el programa no tiene sentido que continue
bar[0] = ( *ifil = (int *) realloc( *ifil, (*nfil)*sizeof(int) ) ) == NULL;
bar[1] = ( *icol = (int *) realloc( *icol, (*ncol)*sizeof(int) ) ) == NULL;
bar[2] = ( *xval = (double *) realloc( *xval, (nval)*sizeof(double) ) ) == NULL;
bar[3] = ( *b = (double *) realloc( *b, (nb)*sizeof(double) ) ) == NULL;

//Si alguno de los realloc retorna NULL, imprimimos un mensaje alusivo y salimos
if(bar[0] || bar[1] || bar[2] || bar[3])
{
printf("No hubo memoria suficiente para hacer realloc a los vectores de una matriz\n");
exit(EXIT_FAILURE);
}

if ( *nfil != *ncol || (*nfil != nval || *ncol != nval)  ) return -2;


//Se toma el mayor valor de los indices de filas/columnas como orden de la matriz
aux = **ifil;
for(foo=1; foo<*nfil; foo++)
if( *(*ifil+foo) > aux )
aux = *(*ifil+foo);
*nfil = aux;

aux = **icol;
for(foo=1; foo<*ncol; foo++)
if( *(*icol+foo) > aux )
aux = *(*icol+foo);
*ncol = aux;
//Para este punto, *nfil tiene el orden de la matriz en cuanto
//a filas y *ncol en cuanto a las columnas y por tanto podemos
//comparar con el segundo if para retornar

//Por defecto se retorna el numero de elementos no nulos encontrados en xval
a = nval;
//De encontrarse un problema de consistencia, cambiamos el valor de a

if ( *ncol != nb ) a = -1;

//Retornamos
return a;
}

Saludos, quizá pronto traiga más de este tema

Nota quiero agradecer a mi profesor quien hizo esfuerzo a lo largo del curso trayendo laboratorios interesantes.

PD: como siempre, acepto sugerencias y consejos