Rev 108 | Blame | Compare with Previous | Last modification | View Log | RSS feed
/* Ira Snyder* 04-20-2005* License: Public Domain** CS331: Analysis of Algorithms* Project #1** Due: 05-11-2005** Changelog follows:** 04-20-2005:* Added code for Classical Method* Added code for Divide-and-Conquer Method* Switched to char (1 byte) from int (4 byte) to increase possible matrix size** 04-21-2005:* Changed creatematrix() to use malloc() instead of calloc() for speed* Added createzeromatrix() which initializes the returned matrix to 0* Changed the base case of multiply_dnq(). This gives a 4x speed up.* Added Strassen's Method** CODE IS FEATURE COMPLETE AT THIS POINT** 08-08-2005* Changed to using a (hopefully) faster way of copying memory, using memcpy()** 08-11-2005* Testing shows that the memcpy() way of doing things is a little faster, but* not by any significant amount. I'm leaving it here and enabled by default* just for fun.*/#include <stdio.h>#include <stdlib.h>#include <time.h>/* function prototypes */void printmatrix (char **matrix, int size);char** creatematrix (int size);char** createzeromatrix (int size);char** freematrix (char **matrix, int size);void fillmatrix (char **matrix, int size);void multiply_classical (char **matrix1, char **matrix2, char **result, int size);void copymatrix (char **source, char **dest, int startrow, int startcol, int size);void multiply_dnq (char **matrix1, char **matrix2, char **result, int size);void multiply_strassen (char **matrix1, char **matrix2, char **result, int size);void copymatrix (char **source, char **dest, int startrow, int startcol, int size);void copymatrixback (char **source, char **dest, int startrow, int startcol, int size);void addmatrices (char **matrix1, char **matrix2, char **result, int size);void submatrices (char **matrix1, char **matrix2, char **result, int size);void startclock (void);void stopclock (void);void printtimetaken (char *funcname);/* global timeval for timing the function calls */struct timeval start_time, end_time, lapsed;int main (void){int size;char **matrix1;char **matrix2;char **result;/* for loop runs in powers of 2 (2, 4, 8, 16, 32 ...) */for (size=2; size<=4096; size*=2){/* create matrices */matrix1 = creatematrix (size);matrix2 = creatematrix (size);result = createzeromatrix (size);/* fill the matrices to be multiplied (with 1's) */fillmatrix (matrix1, size);fillmatrix (matrix2, size);/* print the current iteration's title */printf("Using %d x %d matrix\n", size, size);/* run the classical method */startclock();multiply_classical (matrix1, matrix2, result, size);stopclock();printtimetaken("classical");/* run the divide-and-conquer method */startclock();multiply_dnq (matrix1, matrix2, result, size);stopclock();printtimetaken("div-n-con");/* run the strassen method */startclock();multiply_strassen (matrix1, matrix2, result, size);stopclock();printtimetaken("strassens");printf("\n");/* free all the memory used for the matrices */freematrix (matrix1, size);freematrix (matrix2, size);freematrix (result, size);}return 0;}/* start the timer */void startclock (void){/* gettimeofday(&start_time, (struct timeval*)0); */gettimeofday(&start_time, NULL);}/* stop the timer */void stopclock (void){/* gettimeofday(&end_time, (struct timeval*)0); */gettimeofday(&end_time, NULL);}/* print the time taken to run, in microseconds */void printtimetaken (char *funcname){double total_usecs = (end_time.tv_sec-start_time.tv_sec) * 1000000.0+ (end_time.tv_usec-start_time.tv_usec);printf("%s : %.3lf milliseconds\n", funcname, total_usecs/1000.0);}/* print out the contents of a matrix */void printmatrix (char **matrix, int size){int i,j;for (i=0; i<size; i++){for (j=0; j<size; j++)printf ("%i ", matrix[i][j]);printf("\n");}}/* allocate a matrix with dimensions size x size** NOTE: The values are NOT INITIALIZED** WARNING: These get big fast!* size = 4096 is 16MB* size = 8192 is 64MB* size = 16384 is 256MB* size = 32768 is 1024MB (1GB)*/char** creatematrix (int size){char **matrix;int i;/* try to allocate first part of memory */matrix = (char**) malloc (size * sizeof(char*));/* check to make sure it allocated okay, exit if it didn't */if (matrix == NULL){fprintf (stderr,"out of memory allocating size=%i matrix\n", size);exit (1);}for (i=0; i<size; i++){/* try to allocate the array cells */matrix[i] = (char*) malloc (size * sizeof(char));/* check to make sure it allocated okay, exit if it didn't */if (matrix[i] == NULL){fprintf (stderr, "out of memory allocating size=%i matrix\n", size);exit (1);}}return matrix;}/* allocate a matrix with dimensions size x size** NOTE: The values are initialized to 0** WARNING: These get big fast!* size = 4096 is 16MB* size = 8192 is 64MB* size = 16384 is 256MB* size = 32768 is 1024MB (1GB)*/char** createzeromatrix (int size){char **matrix;int i;/* try to allocate first part of memory */matrix = (char**) calloc (size, sizeof(char*));/* check to make sure it allocated okay, exit if it didn't */if (matrix == NULL){fprintf (stderr,"out of memory allocating size=%i matrix\n", size);exit (1);}for (i=0; i<size; i++){/* try to allocate the array cells */matrix[i] = (char*) calloc (size, sizeof(char));/* check to make sure it allocated okay, exit if it didn't */if (matrix[i] == NULL){fprintf (stderr, "out of memory allocating size=%i matrix\n", size);exit (1);}}return matrix;}/* free all the memory associated with the matrix */char** freematrix (char **matrix, int size){int i;for (i=0; i<size; i++)free (matrix[i]);free (matrix);return NULL;}/* fill a matrix with 1's */void fillmatrix (char **matrix, int size){int i,j;for (i=0; i<size; i++)for (j=0; j<size; j++)matrix[i][j] = 1;/* if you want numbers 0-9 to fill the matrix, use:* matrix[i][j] = (i+j)%10;*/}/* multiply with the classical method (row * column) */void multiply_classical (char **matrix1, char **matrix2, char **result, int size){int i, j, k;for (i=0; i<size; i++)for (j=0; j<size; j++)for (k=0; k<size; k++){#ifdef DEBUGprintf("i=%d j=%d k=%d\n", i,j,k);#endifresult[i][j] += matrix1[i][k] * matrix2[k][j];}}/* Divide-and-Conquer method of matrix multiplication** The matrix is broken up as follows:** |* quad1 | quad2* -----------------* quad3 | quad4* |** NOTE: Only works on matrices of size = 2^n*/void multiply_dnq (char **matrix1, char **matrix2, char **result, int size){char **m1quad1, **m1quad2, **m1quad3, **m1quad4;char **m2quad1, **m2quad2, **m2quad3, **m2quad4;char **resquad1, **resquad2, **resquad3, **resquad4;char **temp1, **temp2, **temp3, **temp4;char **temp5, **temp6, **temp7, **temp8;int newsize = size/2;/* base case */if (size == 2){result[0][0] = matrix1[0][0] * matrix2[0][0] + matrix1[0][1] * matrix2[1][0];result[0][1] = matrix1[0][0] * matrix2[0][1] + matrix1[0][1] * matrix2[1][1];result[1][0] = matrix1[1][0] * matrix2[0][0] + matrix1[1][1] * matrix2[1][0];result[1][1] = matrix1[1][0] * matrix2[0][1] + matrix1[1][1] * matrix2[1][1];return;}/* create the matrices to hold matrix1 */m1quad1 = creatematrix(newsize);m1quad2 = creatematrix(newsize);m1quad3 = creatematrix(newsize);m1quad4 = creatematrix(newsize);/* create the matrices to hold matrix2 */m2quad1 = creatematrix(newsize);m2quad2 = creatematrix(newsize);m2quad3 = creatematrix(newsize);m2quad4 = creatematrix(newsize);/* create the matrices to hold the result */resquad1 = creatematrix(newsize);resquad2 = creatematrix(newsize);resquad3 = creatematrix(newsize);resquad4 = creatematrix(newsize);/* create the temporary matrices for addition */temp1 = creatematrix(newsize);temp2 = creatematrix(newsize);temp3 = creatematrix(newsize);temp4 = creatematrix(newsize);temp5 = creatematrix(newsize);temp6 = creatematrix(newsize);temp7 = creatematrix(newsize);temp8 = creatematrix(newsize);/* split matrix1 into it's 4 quadrants */copymatrix (matrix1, m1quad1, 0, 0, newsize);copymatrix (matrix1, m1quad2, 0, newsize, newsize);copymatrix (matrix1, m1quad3, newsize, 0, newsize);copymatrix (matrix1, m1quad4, newsize, newsize, newsize);/* split matrix2 into it's 4 quadrants */copymatrix (matrix2, m2quad1, 0, 0, newsize);copymatrix (matrix2, m2quad2, 0, newsize, newsize);copymatrix (matrix2, m2quad3, newsize, 0, newsize);copymatrix (matrix2, m2quad4, newsize, newsize, newsize);/* calculate each of the 8 multiplications */multiply_dnq (m1quad1, m2quad1, temp1, newsize);multiply_dnq (m1quad2, m2quad3, temp2, newsize);multiply_dnq (m1quad1, m2quad2, temp3, newsize);multiply_dnq (m1quad1, m2quad4, temp4, newsize);multiply_dnq (m1quad3, m2quad1, temp5, newsize);multiply_dnq (m1quad4, m2quad3, temp6, newsize);multiply_dnq (m1quad3, m2quad2, temp7, newsize);multiply_dnq (m1quad4, m2quad4, temp8, newsize);/* add the 8 multiplications to get the 4 quadrants */addmatrices (temp1, temp2, resquad1, newsize);addmatrices (temp3, temp4, resquad2, newsize);addmatrices (temp5, temp6, resquad3, newsize);addmatrices (temp7, temp8, resquad4, newsize);/* copy the quadrants back into result */copymatrixback (resquad1, result, 0, 0, newsize);copymatrixback (resquad2, result, 0, newsize, newsize);copymatrixback (resquad3, result, newsize, 0, newsize);copymatrixback (resquad4, result, newsize, newsize, newsize);/* free all the temporary matrices */freematrix (m1quad1, newsize);freematrix (m1quad2, newsize);freematrix (m1quad3, newsize);freematrix (m1quad4, newsize);freematrix (m2quad1, newsize);freematrix (m2quad2, newsize);freematrix (m2quad3, newsize);freematrix (m2quad4, newsize);freematrix (resquad1, newsize);freematrix (resquad2, newsize);freematrix (resquad3, newsize);freematrix (resquad4, newsize);freematrix (temp1, newsize);freematrix (temp2, newsize);freematrix (temp3, newsize);freematrix (temp4, newsize);freematrix (temp5, newsize);freematrix (temp6, newsize);freematrix (temp7, newsize);freematrix (temp8, newsize);}/* copies a portion of the source to the destination */void copymatrix (char **source, char **dest, int startrow, int startcol, int size){int i;#ifndef FAST_MEMCPYint j;#endiffor (i=0; i<size; i++)#ifdef FAST_MEMCPYmemcpy (dest[i], source[i+startrow], sizeof(char)*size);#elsefor (j=0; j<size; j++)dest[i][j] = source[i+startrow][j+startcol];#endif}/* copies from the source to a portion of the destination */void copymatrixback (char **source, char **dest, int startrow, int startcol, int size){int i;#ifndef FAST_MEMCPYint j;#endiffor (i=0; i<size; i++)#ifdef FAST_MEMCPYmemcpy (dest[i+startrow], source[i], sizeof(char)*size);#elsefor (j=0; j<size; j++)dest[i+startrow][j+startcol] = source[i][j];#endif}/* adds two matrices together */void addmatrices (char **matrix1, char **matrix2, char **result, int size){int i, j;for (i=0; i<size; i++)for (j=0; j<size; j++)result[i][j] = matrix1[i][j] + matrix2[i][j];}void multiply_strassen (char **matrix1, char **matrix2, char **result, int size){char **m1quad1, **m1quad2, **m1quad3, **m1quad4;char **m2quad1, **m2quad2, **m2quad3, **m2quad4;char **resquad1, **resquad2, **resquad3, **resquad4;char **p, **q, **r, **s, **t, **u, **v;char **temp1, **temp2, **temp3, **temp4, **temp5, **temp6, **temp7;char **temp8, **temp9, **temp10, **temp11, **temp12, **temp13, **temp14;int newsize = size/2;/* base case */if (size == 2){result[0][0] = matrix1[0][0] * matrix2[0][0] + matrix1[0][1] * matrix2[1][0];result[0][1] = matrix1[0][0] * matrix2[0][1] + matrix1[0][1] * matrix2[1][1];result[1][0] = matrix1[1][0] * matrix2[0][0] + matrix1[1][1] * matrix2[1][0];result[1][1] = matrix1[1][0] * matrix2[0][1] + matrix1[1][1] * matrix2[1][1];return;}/* create the matrices to hold matrix1 */m1quad1 = creatematrix(newsize);m1quad2 = creatematrix(newsize);m1quad3 = creatematrix(newsize);m1quad4 = creatematrix(newsize);/* create the matrices to hold matrix2 */m2quad1 = creatematrix(newsize);m2quad2 = creatematrix(newsize);m2quad3 = creatematrix(newsize);m2quad4 = creatematrix(newsize);/* create the matrices to hold the result */resquad1 = creatematrix(newsize);resquad2 = creatematrix(newsize);resquad3 = creatematrix(newsize);resquad4 = creatematrix(newsize);/* create the temporary matrices for addition */temp1 = creatematrix(newsize);temp2 = creatematrix(newsize);temp3 = creatematrix(newsize);temp4 = creatematrix(newsize);temp5 = creatematrix(newsize);temp6 = creatematrix(newsize);temp7 = creatematrix(newsize);temp8 = creatematrix(newsize);temp9 = creatematrix(newsize);temp10 = creatematrix(newsize);temp11 = creatematrix(newsize);temp12 = creatematrix(newsize);temp13 = creatematrix(newsize);temp14 = creatematrix(newsize);p = creatematrix(newsize);q = creatematrix(newsize);r = creatematrix(newsize);s = creatematrix(newsize);t = creatematrix(newsize);u = creatematrix(newsize);v = creatematrix(newsize);/* split matrix1 into it's 4 quadrants */copymatrix (matrix1, m1quad1, 0, 0, newsize);copymatrix (matrix1, m1quad2, 0, newsize, newsize);copymatrix (matrix1, m1quad3, newsize, 0, newsize);copymatrix (matrix1, m1quad4, newsize, newsize, newsize);/* split matrix2 into it's 4 quadrants */copymatrix (matrix2, m2quad1, 0, 0, newsize);copymatrix (matrix2, m2quad2, 0, newsize, newsize);copymatrix (matrix2, m2quad3, newsize, 0, newsize);copymatrix (matrix2, m2quad4, newsize, newsize, newsize);/* calculate all the additions necessary for pqrstuv */addmatrices (m1quad1, m1quad4, temp1, newsize);addmatrices (m2quad1, m2quad4, temp2, newsize);addmatrices (m1quad3, m1quad4, temp3, newsize);submatrices (m2quad2, m2quad4, temp4, newsize);submatrices (m2quad3, m2quad1, temp5, newsize);addmatrices (m1quad1, m1quad2, temp6, newsize);submatrices (m1quad3, m1quad1, temp7, newsize);addmatrices (m2quad1, m2quad2, temp8, newsize);submatrices (m1quad2, m1quad4, temp9, newsize);addmatrices (m2quad3, m2quad4, temp10, newsize);/* call multiply_strassen() */multiply_strassen (temp1, temp2, p, newsize);multiply_strassen (temp3, m2quad1, q, newsize);multiply_strassen (m1quad1, temp4, r, newsize);multiply_strassen (m1quad4, temp5, s, newsize);multiply_strassen (temp6, m2quad4, t, newsize);multiply_strassen (temp7, temp8, u, newsize);multiply_strassen (temp9, temp10, v, newsize);/* find first result quadrant */addmatrices (p, s, temp11, newsize);addmatrices (t, v, temp12, newsize);submatrices (temp11, temp12, resquad1, newsize);/* find second result quadrant */addmatrices (r, t, resquad2, newsize);/* find third result quadrant */addmatrices (q, s, resquad3, newsize);/* find fourth result quadrant */addmatrices (p, r, temp13, newsize);addmatrices (q, u, temp14, newsize);submatrices (temp13, temp14, resquad4, newsize);/* copy the quadrants back into result */copymatrixback (resquad1, result, 0, 0, newsize);copymatrixback (resquad2, result, 0, newsize, newsize);copymatrixback (resquad3, result, newsize, 0, newsize);copymatrixback (resquad4, result, newsize, newsize, newsize);/* free all the temporary matrices */freematrix (m1quad1, newsize);freematrix (m1quad2, newsize);freematrix (m1quad3, newsize);freematrix (m1quad4, newsize);freematrix (m2quad1, newsize);freematrix (m2quad2, newsize);freematrix (m2quad3, newsize);freematrix (m2quad4, newsize);freematrix (resquad1, newsize);freematrix (resquad2, newsize);freematrix (resquad3, newsize);freematrix (resquad4, newsize);freematrix (p, newsize);freematrix (q, newsize);freematrix (r, newsize);freematrix (s, newsize);freematrix (t, newsize);freematrix (u, newsize);freematrix (v, newsize);freematrix (temp1, newsize);freematrix (temp2, newsize);freematrix (temp3, newsize);freematrix (temp4, newsize);freematrix (temp5, newsize);freematrix (temp6, newsize);freematrix (temp7, newsize);freematrix (temp8, newsize);freematrix (temp9, newsize);freematrix (temp10, newsize);freematrix (temp11, newsize);freematrix (temp12, newsize);freematrix (temp13, newsize);freematrix (temp14, newsize);}/* subtracts two matrices (matrix1 - matrix2) */void submatrices (char **matrix1, char **matrix2, char **result, int size){int i, j;for (i=0; i<size; i++)for (j=0; j<size; j++)result[i][j] = matrix1[i][j] - matrix2[i][j];}