Multiplicación de matrices por el algoritmo Strassen

Multiplicación de matrices por el algoritmo Strassen

El algoritmos de Strassen, llamado así en honor del matemático Volker Strassen, es un algoritmo diseñado para la multiplicación de matrices. Este código que os voy a presentar es óptimo para matrices grandes.

Os voy a dejar, en éste artículo, una propuesta de código realizada en mis tiempos de universidad. Está programado en lenguaje c.

#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <process.h>
#include <dos.h>
#include <alloc.h>
#include <math.h>

#include "array2d.h"
#include "lib_cron.h"

/*-------------------------------------------------------------------------*/
/*                       CABECERAS DE FUNCIONES                            */
/*-------------------------------------------------------------------------*/

void lee (imatriz2d matrix, imatriz2d matrix2, int filas1, int columnas1,
					int filas2, int columnas2, FILE *fp);
int maximo (int x1, int j1, int x2, int j2);
int pot_dos (int num);
void strassen (int *p1a, int dim1, int *p1b, int dim2, int *p1c, int dim3,
							 int semidimen);
void suma (int *p1a, int dim1, int *p1b, int dim2, int *p1c, int dim3,
					 int dimen);
void resta (int *p1a, int dim1, int *p1b, int dim2, int *p1c, int dim3,
						int dimen);
void cambiasigno (imatriz2d matriz, int dimen);

/*-------------------------------------------------------------------------*/
/*                           PROGRAMA PRINCIPAL                            */
/*-------------------------------------------------------------------------*/

int main (void) {

	int i, j;
	unsigned int max;
	int tama;
	int filas1, columnas1, filas2, columnas2;
	imatriz2d matriz1, matriz2, matriz3;
	int largo;
	FILE *fp;
	struct  time t1, t2;
	ldiv_t segundo;
	double segundo2;
	unsigned long int centesimas;


	clrscr ();

	if ((fp = fopen ("datos.txt", "r")) == NULL){
		printf ("\n Error en la apertura del fichero");
		exit (1);
	}

	fscanf (fp, "%d %d %d %d ", &amp;filas1, &amp;columnas1, &amp;filas2, &amp;columnas2);

	if (columnas1 != filas2){
		printf ("\nLas dos matrices no son multiplicables");
		exit (1);
	}
	max = maximo (filas1, columnas1, filas2, columnas2);

	if ((tama = pot_dos(max)) == -1){
		printf ("\nUna de las matrices sobrepasa el tama¤o permitido");
		exit (1);
	}

	if (tama != 0) { /*no es potencia de 2*/
		matriz1 = icreamatriz2d (tama, tama);
		matriz2 = icreamatriz2d (tama, tama);
		matriz3 = icreamatriz2d (tama, tama);
		largo = tama;
	}
	else {  /* es potencia de dos*/
		matriz1 = icreamatriz2d (max, max);
		matriz2 = icreamatriz2d (max, max);
		matriz3 = icreamatriz2d (max, max);
		largo = max;
	}

	/* Llegados a este punto ya tenemos en memoria las dos matrices a
		multiplicar, s¢lo falta rellenarlas con las matrices originales y
		con una banda de 0's en caso de ser necesario */


	for (i=0; i<largo; i++){
		for (j=0; j<largo; j++){
			matriz1[i][j] = 0;
			matriz2[i][j] = 0;
			matriz3[i][j] = 0;
		}
	}

	lee (matriz1, matriz2, filas1, columnas1, filas2, columnas2, fp);

	for (i=0; i<filas1; i++){
		printf ("\n");
		for (j=0; j<columnas1; j++){
			printf("%d ", matriz1[i][j]);
		}
	}

	printf ("\n\n");

	for (i=0; i<filas2; i++){
		printf ("\n");
		for (j=0; j<columnas2; j++){
			printf("%d ", matriz2[i][j]);
		}
	}

	cronostart (&amp;t1);

	strassen (*matriz1,largo, *matriz2, largo, *matriz3, largo, largo/2);

	cronostop (&amp;t2);
	printf ("\a");
	centesimas = cronodifcent (&amp;t1, &amp;t2);
	segundo = ldiv (centesimas, 100L);
	printf ("\n\nEl n£mero de cent‚simas es: %lu", centesimas);
	segundo2 = segundo.quot+segundo.rem*0.01;
	printf ("\nEl n£mero de segundos es: %.2lf", segundo2);

	printf ("\n\n");

	for (i=0; i<filas1; i++){
		printf ("\n");
		for (j=0; j<columnas2; j++){
			printf("%d ", matriz3[i][j]);
		}
	}

	ifreematriz2d (&amp;matriz1);
	ifreematriz2d (&amp;matriz2);
	ifreematriz2d (&amp;matriz3);

	fclose (fp);
	fp = NULL;

	return 0;
} /*del main*/

/*-------------------------------------------------------------------------*/

int maximo (int x1, int j1, int x2, int j2){

	int aus, aus2;

	if (x1 >= j1) aus = x1;
	else aus = j1;
	if (x2 >= j2) aus2 = x2;
	else aus2 = j2;
	if (aus >= aus2) return (aus);
	else return (aus2);

}

/*-------------------------------------------------------------------------*/

/* La funci¢n pot_dos devuelve 0 en caso de que el argumento sea potencia de
	dos, -1 en caso de que num sea menor que todos los elementos del array y
	la potencia de dos m s grande que ‚l en caso de que no lo sea. */

int pot_dos (int num){

	int array [11] = {1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024};
	int i;


	for (i=0; i<11; i++){
		if (num == array[i]) return (0);
		else if (num < array[i]) return (array[i]);
	}
	return (-1);
}

/*-------------------------------------------------------------------------*/

void lee (imatriz2d matrix1, imatriz2d matrix2, int filas1, int columnas1,
					int filas2, int columnas2, FILE *fp){

	int i, j;

	for (i=0; i<filas1; i++){
		for (j=0; j<columnas1; j++){
			fscanf(fp, "%d ", &amp;matrix1[i][j]);
		}
	}

	for (i=0; i<filas2; i++){
		for (j=0; j<columnas2; j++){
			fscanf(fp, "%d ", &amp;matrix2[i][j]);
		}
	}
}

/*-------------------------------------------------------------------------*/

void strassen (int *p1a, int dim1, int *p1b, int dim2, int *p1c, int dim3,
							 int semidimen){

	int i,j;

	int *p2a, *p3a, *p4a;
	int *p2b, *p3b, *p4b;
	int *p2c, *p3c, *p4c;


	int aus1, aus2, aus3, m;
	imatriz2d maus1, maus2, maus3;


	p2a = p1a+semidimen;
	p3a = p1a+dim1*semidimen;
	p4a = p3a+semidimen;

	p2b = p1b+semidimen;
	p3b = p1b+dim2*semidimen;
	p4b = p3b+semidimen;

	p2c = p1c+semidimen;
	p3c = p1c+dim3*semidimen;
	p4c = p3c+semidimen;


	if (semidimen == 1){
		m = p1a[0]*p1b[0];  /* m2 */
		p1c[0]+= m;
		p2c[0]+= m;
		p3c[0]+= m;
		p4c[0]+= m;

		m = p2a[0]*p3b[0];  /* m3 */
		p1c[0]+= m;

		aus1 = p1a[0]-p3a[0];
		aus2 = p4b[0]-p2b[0];
		m = aus1*aus2;            /* m4 */
		p3c[0]+= m;
		p4c[0]+= m;

		aus3 = -aus1+p4a[0];
		aus1 = aus2+p1b[0];
		m = aus3*aus1;            /* m1 */
		p2c[0]+= m;
		p3c[0]+= m;
		p4c[0]+= m;

		m = (aus3+p1a[0])*(-aus1+p4b[0]);    /* m5 */
		p2c[0]+= m;
		p4c[0]+= m;

		m = (-aus3+p2a[0])*p4b[0];     /* m6 */
		p2c[0]+= m;

		m = (aus1-p3b[0])*p4a[0];        /* m7 */
		p3c[0]-= m;

	}
	else{
		maus1 = icreamatriz2d (semidimen, semidimen);
		for (i=0; i<semidimen; i++)
			for (j=0; j<semidimen; j++)
				maus1[i][j] = 0;

								/* m2 */
		strassen (p1a, dim1, p1b, dim2, *maus1, semidimen, semidimen/2);

		suma (*maus1, semidimen, p1c, dim3, p1c, dim3, semidimen);
		suma (*maus1, semidimen, p2c, dim3, p2c, dim3, semidimen);
		suma (*maus1, semidimen, p3c, dim3, p3c, dim3, semidimen);
		suma (*maus1, semidimen, p4c, dim3, p4c, dim3, semidimen);

								/* m3 */
		for (i=0; i<semidimen; i++)
			for (j=0; j<semidimen; j++)
				maus1[i][j] = 0;

		strassen (p2a, dim1, p3b, dim2, *maus1, semidimen, semidimen/2);

		suma (*maus1, semidimen, p1c, dim3, p1c, dim3, semidimen);

								/* m4 */
		maus2 = icreamatriz2d (semidimen, semidimen);
							 //maus2 = a11-a21
		resta (p1a, dim1, p3a, dim1, *maus2, semidimen, semidimen);


		maus3 = icreamatriz2d (semidimen, semidimen);
							//maus3 = b22-b12
		resta (p4b, dim2, p2b, dim2, *maus3, semidimen, semidimen);

		for (i=0; i<semidimen; i++)
			for (j=0; j<semidimen; j++)
				maus1[i][j] = 0;
								//m4= maus2*maus3
		strassen (*maus2, semidimen, *maus3, semidimen, *maus1, semidimen,
							semidimen/2);

		suma (*maus1, semidimen, p3c, dim3, p3c, dim3, semidimen);
		suma (*maus1, semidimen, p4c, dim3, p4c, dim3, semidimen);


								/* m1  */

					 //maus2 = a22-aus1
		resta (p4a, dim1, *maus2, semidimen, *maus2, semidimen, semidimen);

				 //maus3 = aus2+b11
		suma (*maus3, semidimen, p1b, dim2, *maus3, semidimen, semidimen);

		for (i=0; i<semidimen; i++)
			for (j=0; j<semidimen; j++)
				maus1[i][j] = 0;
					//m1= maus2*maus3
		strassen (*maus2, semidimen, *maus3, semidimen,	*maus1, semidimen,
							semidimen/2);

		suma (*maus1, semidimen, p2c, dim3, p2c, dim3, semidimen);
		suma (*maus1, semidimen, p3c, dim3, p3c, dim3, semidimen);
		suma (*maus1, semidimen, p4c, dim3, p4c, dim3, semidimen);


							 /*  m5 */

							//maus2= aus3+a11
		suma  (*maus2, semidimen, p1a, dim1, *maus2, semidimen, semidimen);

					 //maus3= b22-aus1
		resta ( p4b, dim2, *maus3, semidimen, *maus3, semidimen, semidimen);

		for (i=0; i<semidimen; i++)
			for (j=0; j<semidimen; j++)
				maus1[i][j] = 0;

					//m5 = maus2*maus3
		strassen (*maus2, semidimen, *maus3, semidimen, *maus1, semidimen,
							semidimen/2);

		suma (*maus1, semidimen, p2c, dim3, p2c, dim3, semidimen);
		suma (*maus1, semidimen, p4c, dim3, p4c, dim3, semidimen);


							 /* m6 */
					//maus2=aus3
		resta (*maus2, semidimen, p1a, dim1, *maus2, semidimen, semidimen);

						 //maus2=a12-aus3
		resta (p2a, dim1, *maus2, semidimen, *maus2, semidimen, semidimen);

		for (i=0; i<semidimen; i++)
			for (j=0; j<semidimen; j++)
				maus1[i][j] = 0;
						// m6 = maus2*b22
		strassen (*maus2, semidimen, p4b, dim2, *maus1, semidimen,
							semidimen/2);

		suma (*maus1, semidimen, p2c, dim3, p2c, dim3, semidimen);

		//calcular de nuevo aus1; maus3=b22-aus1->aus1=b22-maus3
		resta (p4b, dim2, *maus3, semidimen, *maus3, semidimen, semidimen);
					//aus1-b21
		resta (*maus3, semidimen, p3b, dim3, *maus3, semidimen, semidimen);

		for (i=0; i<semidimen; i++)
			for (j=0; j<semidimen; j++)
				maus1[i][j] = 0;
				//m7=maus3*a22
		strassen (*maus3, semidimen, p4a, dim1, *maus1, semidimen,
							semidimen/2);

		cambiasigno (maus1, semidimen);
		suma (*maus1, semidimen, p3c, dim3, p3c, dim3, semidimen);
/* Liberar memoria */
		ifreematriz2d (&amp;maus1);
		ifreematriz2d (&amp;maus2);
		ifreematriz2d (&amp;maus3);

	}
} /* De la funci¢n strassen */

/*-------------------------------------------------------------------------*/
void suma (int *p1a, int dim1,int *p1b, int dim2, int *p1c, int dim3,
					 int dimen){

		int i, j;

		for (i=0; i<dimen; i++){
			for (j=0; j<dimen; j++){
				p1c[j] = p1a[j]+p1b[j];
			}
			p1a+= dim1;
			p1b+= dim2;
			p1c+= dim3;
		}
}

/*-------------------------------------------------------------------------*/

void resta (int *p1a, int dim1, int *p1b, int dim2, int *p1c, int dim3,
						int dimen){

	 int i, j;

	 for (i=0; i<dimen; i++){
		for (j=0; j<dimen; j++){
			p1c[j] = p1a[j]-p1b[j];
		}
		p1a+= dim1;
		p1b+= dim2;
		p1c+= dim3;
	 }
}

/*-------------------------------------------------------------------------*/

void cambiasigno (imatriz2d matriz, int dimen){

	int i, j;

	for  (i=0; i<dimen; i++){
		for (j=0; j<dimen; j++){
			matriz[i][j] = -matriz[i][j];
		}
	}
}

Termino dejando mi solicitud de que comentéis en la sección de comentarios vuestra opinión o, directamente, compartáis vuestra versión de códigos.

No Comments

Post a Comment