Subversion Repositories programming

/* 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 DEBUG

                printf("i=%d j=%d k=%d\n", i,j,k);

                #endif

                result[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_MEMCPY

    int j;

#endif

    for (i=0; i<size; i++)

#ifdef FAST_MEMCPY

        memcpy (dest[i], source[i+startrow], sizeof(char)*size);

#else

        for (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_MEMCPY

    int j;

#endif

    for (i=0; i<size; i++)

#ifdef FAST_MEMCPY

        memcpy (dest[i+startrow], source[i], sizeof(char)*size);

#else

        for (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];

}