test.cpp

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <iostream>
#include <cmath>
#include <omp.h>
#include <fstream>
#define leave 32
using namespace std;
fstream ff;

static void mm_generate(double* matA, double* matB, double* matC, const int M, const int N, const int K,
                        const int strideA, const int strideB, const int strideC) {
	#pragma omp parallel for
	for (int i = 0; i < M; i++) {
		for (int j = 0; j < N; j++) {
			double sum = 0.0f;
			for (int k = 0; k < K; k++) {
				sum += matA[i*strideA + k] * matB[k*strideB + j];
			}
			matC[i*strideC + j] = sum;
		}
	}
}


static void mm_CoppersmithWinograd(double* matA, double* matB, double* matC, const int M, const int N, const int K,
                             const int strideA, const int strideB, const int strideC) {
	if((M <= leave) || (M%2 != 0 || N%2 != 0 || K%2 != 0)) {
		return mm_generate(matA, matB, matC, M, N, K, strideA, strideB, strideC);
	}
	double* S1 = (double*) malloc((M/2) * (K/2) * sizeof(double));
	double* S2 = (double*) malloc((M/2) * (K/2) * sizeof(double));
	double* S3 = (double*) malloc((M/2) * (K/2) * sizeof(double));
	double* S4 = (double*) malloc((M/2) * (K/2) * sizeof(double));
	#pragma omp task shared(S1) shared(S2) shared(S3) shared(S4)
	// #pragma omp task shared(S2)
	// #pragma omp task shared(S3)
	// #pragma omp task shared(S4)
	// shared(S2) shared(S3) shared(S4)
	{	
		for (int i = 0; i < M/2; i++) {
			for (int j = 0; j < K/2; j++) {
				int idxA, offset, idxS = i * (K/2) + j;
				//S1     = A21 + A22
				idxA     = (i + (M/2)) * strideA + j;
				offset   = K/2;
				S1[idxS] = matA[idxA] + matA[idxA + offset];
				//S2     = S1 - A11
				idxA     = i * strideA + j;
				S2[idxS] = S1[idxS] - matA[idxA];
				//S3     = A11 - A21
				offset   = (M/2) * strideA;
				S3[idxS] = matA[idxA] - matA[idxA + offset];
				//S4     = A12 - S2
				idxA     = i * strideA + (K/2) + j;
				S4[idxS] = matA[idxA] - S2[idxS];
			}
		}
	}


	double* T1 = (double*) malloc((K/2) * (N/2) * sizeof(double));
	double* T2 = (double*) malloc((K/2) * (N/2) * sizeof(double));
	double* T3 = (double*) malloc((K/2) * (N/2) * sizeof(double));
	double* T4 = (double*) malloc((K/2) * (N/2) * sizeof(double));
	// #pragma omp task shared(T1) shared(T2) shared(T3) shared(T4)
	// #pragma omp task shared(T1)
	// #pragma omp task shared(T2)
	// #pragma omp task shared(T3)
	// #pragma omp task shared(T4)
	{
		#pragma omp parallel for
		for (int i = 0; i < K/2; i++) {
			for (int j = 0; j < N/2; j++) {
				int idxB, offset, idxT = i * (N/2) + j;
				//T1     = B12 - B11
				idxB     = i * strideB + j;
				offset   = (N/2);
				T1[idxT] = matB[idxB + offset] - matB[idxB];
				//T2     = B22 - T1
				idxB     = (i + (K/2)) * strideB + (N/2) + j;
				T2[idxT] = matB[idxB] - T1[idxT];
				//T3     = B22 - B12
				idxB     = i * strideB + (N/2) + j;
				offset   = ((K/2)) * strideB;
				T3[idxT] = matB[idxB + offset] - matB[idxB];
				//T4     = T2 - B21
				idxB     = (i + (K/2)) * strideB + j;
				T4[idxT] = T2[idxT] - matB[idxB];
			}
		}
	}

	// cout << "Hello" << endl;
	//M1 = A11 * B11
	double* M1 = (double*) malloc((M/2) * (N/2) * sizeof(double));
	{
		mm_CoppersmithWinograd(matA, matB, &M1[0], M/2, N/2, K/2, strideA, strideB, N/2);
	}
	//M2 = A12 * B21
	double* M2 = (double*) malloc((M/2) * (N/2) * sizeof(double));
	mm_CoppersmithWinograd(&matA[K/2], &matB[(K/2)*strideB], &M2[0], M/2, N/2, K/2, strideA, strideB, N/2);
	//M3 = S4 * B22
	double* M3 = (double*) malloc((M/2) * (N/2) * sizeof(double));
	mm_CoppersmithWinograd(&S4[0], &matB[(K/2) * strideB + (N/2)], &M3[0], M/2, N/2, K/2, K/2, strideB, N/2);
	//M4 = A22 * T4
	double* M4 = (double*) malloc((M/2) * (N/2) * sizeof(double));
	mm_CoppersmithWinograd(&matA[(M/2) * strideA + (K/2)], &T4[0], &M4[0], M/2, N/2, K/2, strideA, N/2, N/2);
	//M5 = S1 * T1
	double* M5 = (double*) malloc((M/2) * (N/2) * sizeof(double));
	mm_CoppersmithWinograd(&S1[0], &T1[0], &M5[0], M/2, N/2, K/2, K/2, N/2, N/2);
	//M6 = S2 * T2
	double* M6 = (double*) malloc((M/2) * (N/2) * sizeof(double));
	mm_CoppersmithWinograd(&S2[0], &T2[0], &M6[0], M/2, N/2, K/2, K/2, N/2, N/2);
	//M7 = S3 * T3
	double* M7 = (double*) malloc((M/2) * (N/2) * sizeof(double));
	mm_CoppersmithWinograd(&S3[0], &T3[0], &M7[0], M/2, N/2, K/2, K/2, N/2, N/2);
	//C11 = U1 = M1 + M2
	//C12 = U5 = U4 + M3 = U2 + M5 + M3 = M1 + M6 + M5 + M3
	//C21 = U6 = U3 - M4 = U2 + M7 - M4 = M1 + M6 + M7 - M4
	//C22 = U7 = U3 + M5 = U2 + M7 + M5 = M1 + M6 + M7 + M5
	for (int i = 0; i < M/2; i++) {
		for (int j = 0; j < N/2; j++) {
			int idx = i * (N/2) + j;
			matC[i*strideC + j] = M1[idx] + M2[idx];
			matC[i*strideC + j + (N/2)] = M1[idx] + M6[idx] + M5[idx] + M3[idx];
			matC[(i+(M/2))*strideC + j] = M1[idx] + M6[idx] + M7[idx] - M4[idx];
			matC[(i+(M/2))*strideC + j + (N/2)] = M1[idx] + M6[idx] + M7[idx] + M5[idx];
		}
	}

    free(S1);           S1=NULL;
    free(S2);           S2=NULL;
    free(S3);           S3=NULL;
    free(S4);           S4=NULL;
    free(T1);           T1=NULL;
    free(T2);           T2=NULL;
    free(T3);           T3=NULL;
    free(T4);           T4=NULL;
    free(M1);           M1=NULL;
    free(M2);           M2=NULL;
    free(M3);           M3=NULL;
    free(M4);           M4=NULL;
    free(M5);           M5=NULL;
    free(M6);           M6=NULL;
    free(M7);           M7=NULL;
}
// static void mm_CoppersmithWinograd(double* matA, double* matB, double* matC, const int M, const int N, const int K) {
// 	mm_CoppersmithWinograd(matA, matB, matC, M, N, K, K, N, N);
// }


/*
 * It is the test function
 */
static void mm_test(int M, int N, int K, int rangeTop, int rangeButtom) {
	// unsigned seed = time(0);
	// srand(1017);
	int randomNum[5] = {1001, 8808, 9909, 1242, 1293};
	double start,end;
	double time1 = 0.0;
	double time2 = 0.0;
	double time3 = 0.0;
	// for (int i = 0; i < 1; i++) {
	double * mA = (double*) malloc(M*K*sizeof(double));
	double * mB = (double*) malloc(K*N*sizeof(double));
	double * mC = (double*) malloc(M*N*sizeof(double));
	double * mD = (double*) malloc(M*N*sizeof(double));
	double * mE = (double*) malloc(M*N*sizeof(double));
	srand(1017);
	for (int j = 0; j < M*K; j++) {
		int t = rand() % (rangeTop - rangeButtom) + rangeButtom + 1;
		mA[j] = pow(-1, t) * rand() / double(RAND_MAX);
	}
	srand(1001);
	for (int j = 0; j < K*N; j++) {
		int t = rand() % (rangeTop - rangeButtom) + rangeButtom + 1;
		mB[j] = pow(-1, t) * rand() / double(RAND_MAX);
	}
	start = omp_get_wtime();
	omp_set_num_threads(4);
	#pragma omp parallel
	{
		#pragma omp single
		{
			mm_CoppersmithWinograd(mA, mB, mE, M, N, K, K, N, N);
		}
	}
	end = omp_get_wtime();
	double endtime3 = end-start;
	// time3 = time3 + endtime3;
	// printf("Winograd%d time: %lfms\n", i, endtime3*1000);
	cout << endtime3*1000 << endl;
	ff << endtime3*1000 << endl;
	// }
}
int main() {
	int M, N, K, rangeTop, rangeButtom;
	int num[57] = {32, 64, 128, 256, 512, 1024, 2048, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000, 2100, 2200, 2300, 2400, 2500, 2600, 2700, 2800, 2900, 3000, 3100, 3200, 3300, 3400, 3500, 3600, 3700, 3800, 3900, 4000, 4100, 4200, 4300, 4400, 4500, 4600, 4700, 4800, 4900, 5000};
	ff.open("test.txt", ios::out);
	ff << "OMP -- CW" << endl;
	ff.close();
	ff.open("test.txt", ios::out|ios::app);
	rangeTop = 2;
	rangeButtom = 0;
	for (int i = 0; i < 57; i++) {
		// cout << num[i] << endl;
		mm_test(num[i], num[i], num[i], rangeTop, rangeButtom);
		// cout << endl;
	}
	return 0;
}