/*
 *  =====================================================================
 *  Ordering of loops in matrix multiplication
 *  =====================================================================
 *  Copyright (C) 1993-96 by Mark Austin and David Mazzoni.
 *                                                                     
 *  This software is provided "as is" without express or implied warranty.
 *  Permission is granted to use this software on any computer system,
 *  and to redistribute it freely, subject to the following restrictions:
 * 
 *  1. The authors are not responsible for the consequences of use of
 *     this software, even if they arise from defects in the software.
 *  2. The origin of this software must not be misrepresented, either
 *     by explicit claim or by omission.
 *  3. Altered versions must be plainly marked as such, and must not
 *     be misrepresented as being the original software.
 *  4. This notice is to remain intact.
 *
 *  Written By: Mark Austin                                  January 1993
 *  =====================================================================
 */

#include <stdio.h>
#include "miscellaneous.h"
#include "defs.h"

enum { MatrixSize = 50};

double a[ MatrixSize ][ MatrixSize ];
double b[ MatrixSize ][ MatrixSize ];
double c[ MatrixSize ][ MatrixSize ];

int main( void ) {
int ii, ij;

    /* [a] : Initialize Matrices */

    for(ii = 1; ii <= MatrixSize; ii = ii + 1) 
        for(ij = 1; ij <= MatrixSize; ij = ij + 1)  {
            a[ii-1][ij-1] = (double) ii + 1;
            b[ii-1][ij-1] = (double) ii + 1;
    }

    /* [b] : Multiply Matrices */

    mult_slow();
    mult_fast();
}

mult_slow() {
int ii, ij, ik;

    for(ij = 1; ij <= MatrixSize; ij = ij + 1) 
    for(ik = 1; ik <= MatrixSize; ik = ik + 1) 
    for(ii = 1; ii <= MatrixSize; ii = ii + 1)
        c[ii-1][ij-1] = c[ii-1][ij-1] + a[ii-1][ik-1]*b[ik-1][ij-1];

}

mult_fast() {
int ii, ij, ik;

    for(ii = 1; ii <= MatrixSize; ii = ii + 1)
    for(ik = 1; ik <= MatrixSize; ik = ik + 1) 
    for(ij = 1; ij <= MatrixSize; ij = ij + 1) 
        c[ii-1][ij-1] += a[ii-1][ik-1]*b[ik-1][ij-1];

}