lunes, julio 26, 2021

Multiplicación de matrices por el algoritmo Strassen

Deberías leer

Antonio García Pratshttp://antoniogarciaprats.eu/
Estudios de Licenciatura en CC Quimicas e Ingeniería Técnica en Informática de Gestión por la Universidad de Jaén. Social Media autodidacta. Siempre en búsqueda de nuevos retos que proporcionen algún sentido a mi existencia.

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 ", &filas1, &columnas1, &filas2, &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 (&t1);
	strassen (*matriz1,largo, *matriz2, largo, *matriz3, largo, largo/2);
	cronostop (&t2);
	printf ("\a");
	centesimas = cronodifcent (&t1, &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 (&matriz1);
	ifreematriz2d (&matriz2);
	ifreematriz2d (&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 ", &matrix1[i][j]);
		}
	}
	for (i=0; i<filas2; i++){
		for (j=0; j<columnas2; j++){
			fscanf(fp, "%d ", &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 (&maus1);
		ifreematriz2d (&maus2);
		ifreematriz2d (&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.

- Advertisement -spot_img

Más publicaciones

Dejar respuesta

Please enter your comment!
Please enter your name here

- Advertisement -spot_img

Últimas publicaciones