#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>
#include <omp.h>

double gettime(void) {
   struct timeval tv;
   gettimeofday(&tv, NULL);
   return tv.tv_sec + 1e-6 * tv.tv_usec;
}

void matfill(long N, double *mat, double val) {
   long i, j;
   
   for(i = 0; i < N; i ++){
      for(j = 0; j < N; j ++){
                mat[i * N + j] = val;
      }
   }
}

void matmul(long N, double *a, double *b, double *c) {
/*************************************************************************
 ***                  Matrix-Matrix Multiplication                     ***
 ***  (cache blocking,loop unrolling,OpenMP tasks,Strassen algorithm)  ***
 ***                                                                   ***
 ***                  HP-SEE Computing Challenge                       ***
 ***      http://www.hp-see.eu/hp-see-computing-challenge              ***
 ***                                                                   ***
 ***            (C) Alexandros Sokratis Papadakis 2012                 ***
 ***                   papadakis@ceid.upatras.gr                       ***
 *************************************************************************/
/*******************************************************
*                   internal constants                 *
*******************************************************/
        #define NUM_THREADS      8
        #define BLOCK_SIZE       50 //IMPORTANT: if BLOCK_SIZE!=50, correct the loop unrolling of core_mul,core_add,core_sub function
	#define THRESHOLD_MATRIX 7000
	int THRESHOLD_MULT;
/*******************************************************
*                   internal functions                 *
*******************************************************/
        void matmul2(long N, double *X, double *Y, double *Z) {
                // S=X*Y classic 
                // X,Y square matrices of size N
                int i,j,k;
                int tol1=N/NUM_THREADS;
                int thread;

                double *b_t=(double *)malloc(N*N*sizeof(double));
                if(b_t==NULL){
                        printf("\nnot enough memory\n");
                        exit(1);
                } 
              
                // save matrix b as b-transpose to perform better dot products (C language: row wise memory)
                for(thread=0;thread<NUM_THREADS;thread++)
                        #pragma omp task firstprivate(thread) private(i,j)
                        {
                                for (i=thread*tol1;i<(thread+1)*tol1;i++)
                                        for(j=0;j<N;j++)
                                                b_t[j*N+i]=Y[i*N+j];
                        }

                #pragma omp taskwait

                for(i=NUM_THREADS*tol1;i<N;i++)
                        for(j=0;j<N;j++)
                                 b_t[j*N+i]=Y[i*N+j];


                for(thread=0;thread<NUM_THREADS;thread++)
                        #pragma omp task firstprivate(thread) private(i,j)
                        {
                                for (i=thread*tol1;i<(thread+1)*tol1;i++){
                        		for (j = 0; j < N; j++) {
                                		double sum = 0.0;
                                		for ( k = 0; k < N; k++)
                                        		sum += X[i*N + k]*b_t[j*N + k];
                                		Z[i*N + j] = sum;
                        		}
				}
                        }

                #pragma omp taskwait

                for(i=NUM_THREADS*tol1;i<N;i++){
                        for (j = 0; j < N; j++) {
                                double sum = 0.0;
                                for ( k = 0; k < N; k++)
                                        sum += X[i*N + k]*b_t[j*N + k];
                                Z[i*N + j] = sum;
                        }
		}

                free(b_t);
        }

	void core_mul(int i,int num,int N,double *a,double *b_t,double *c){
		int j, k, m, C, fi, fj, fii, fjj,  el, t1, t2, kk, u, l, pp;
                                
		fi=i*BLOCK_SIZE*N;
                for (j = 0; j < num; j++) {
                	fj=j*BLOCK_SIZE;
                        for (k = 0; k < BLOCK_SIZE; k++) {
                        	fii=fi+k*N;
                                for (m = 0; m < BLOCK_SIZE; m++) {
                                	fjj=fj+m;
                                        el=fii+fjj;
					double sum=0.0;
                                        for (C = 0; C < num; C++) {
                                                kk=C*BLOCK_SIZE;
                                                t1=fii+kk;
                                                t2=fjj*N+kk;
						//for (pp=0; pp < BLOCK_SIZE; pp++)
                                                //	sum+=a[t1+pp]*b_t[t2+pp];
							
						//
						// loop unrolling the above for-loop
						// IMPORTANT: if BLOCK_SIZE!=50 you need to do the appropriate changes
                                                sum+=a[t1]*b_t[t2];
                                                sum+=a[t1+1]*b_t[t2+1];
                                                sum+=a[t1+2]*b_t[t2+2];
                                                sum+=a[t1+3]*b_t[t2+3];
                                                sum+=a[t1+4]*b_t[t2+4];
                                                sum+=a[t1+5]*b_t[t2+5];
                                                sum+=a[t1+6]*b_t[t2+6];
                                                sum+=a[t1+7]*b_t[t2+7];
                                                sum+=a[t1+8]*b_t[t2+8];
                                                sum+=a[t1+9]*b_t[t2+9];
                                                sum+=a[t1+10]*b_t[t2+10];
                                                sum+=a[t1+11]*b_t[t2+11];
                                                sum+=a[t1+12]*b_t[t2+12];
                                                sum+=a[t1+13]*b_t[t2+13];
                                                sum+=a[t1+14]*b_t[t2+14];
                                                sum+=a[t1+15]*b_t[t2+15];
                                                sum+=a[t1+16]*b_t[t2+16];
                                                sum+=a[t1+17]*b_t[t2+17];
                                                sum+=a[t1+18]*b_t[t2+18];
                                                sum+=a[t1+19]*b_t[t2+19];
                                                sum+=a[t1+20]*b_t[t2+20];
                                                sum+=a[t1+21]*b_t[t2+21];
                                                sum+=a[t1+22]*b_t[t2+22];
                                                sum+=a[t1+23]*b_t[t2+23];
                                                sum+=a[t1+24]*b_t[t2+24];
                                                sum+=a[t1+25]*b_t[t2+25];
                                                sum+=a[t1+26]*b_t[t2+26];
                                                sum+=a[t1+27]*b_t[t2+27];
                                                sum+=a[t1+28]*b_t[t2+28];
                                                sum+=a[t1+29]*b_t[t2+29];
                                                sum+=a[t1+30]*b_t[t2+30];
                                                sum+=a[t1+31]*b_t[t2+31];
                                                sum+=a[t1+32]*b_t[t2+32];
                                                sum+=a[t1+33]*b_t[t2+33];
                                                sum+=a[t1+34]*b_t[t2+34];
                                                sum+=a[t1+35]*b_t[t2+35];
                                                sum+=a[t1+36]*b_t[t2+36];
                                                sum+=a[t1+37]*b_t[t2+37];
                                                sum+=a[t1+38]*b_t[t2+38];
                                                sum+=a[t1+39]*b_t[t2+39];
                                                sum+=a[t1+40]*b_t[t2+40];
                                                sum+=a[t1+41]*b_t[t2+41];
                                                sum+=a[t1+42]*b_t[t2+42];
                                                sum+=a[t1+43]*b_t[t2+43];
                                                sum+=a[t1+44]*b_t[t2+44];
                                                sum+=a[t1+45]*b_t[t2+45];
                                                sum+=a[t1+46]*b_t[t2+46];
                                                sum+=a[t1+47]*b_t[t2+47];
                                                sum+=a[t1+48]*b_t[t2+48];
                                                sum+=a[t1+49]*b_t[t2+49];
					}
					c[el]=sum;
				}
			}
		}
	}
        void matmul1(long N, double *a, double *b, double *c) {
                // c=a*b classic way enhanced with cache blocking,loop unrolling
                // a,b square matrices size N 
                int i,j;
                int num=N/BLOCK_SIZE;
                int tol1=N/NUM_THREADS;
                int thread;

                double *b_t = (double *) malloc(N * N * sizeof(double));
                if(b_t==NULL){
                        printf("\nnot enough memory\n");
                        exit(1);
                }

                if(N%BLOCK_SIZE!=0){
                        matmul2(N,a,b,c);
                        return;
                }

                // save matrix b as b-transpose to perform better dot products (C language: row wise memory)
                for(thread=0;thread<NUM_THREADS;thread++)
                        #pragma omp task firstprivate(thread) private(i,j)
                        {
                                for (i=thread*tol1;i<(thread+1)*tol1;i++)
                                        for(j=0;j<N;j++)
                                                b_t[j*N+i]=b[i*N+j];
                        }

                #pragma omp taskwait

                for(i=NUM_THREADS*tol1;i<N;i++)
                        for(j=0;j<N;j++)
                                 b_t[j*N+i]=b[i*N+j];

                tol1=num/NUM_THREADS;

                for(thread=0;thread<NUM_THREADS;thread++)
                        #pragma omp task firstprivate(thread) private(i)
                        {
                                for (i=thread*tol1;i<(thread+1)*tol1;i++)
                                        core_mul(i,num,N,a,b_t,c);
                        }


                #pragma omp taskwait

                for(i=NUM_THREADS*tol1;i<num;i++)
                        core_mul(i,num,N,a,b_t,c);

                free(b_t);
        }

	void core_add(int i,int num,long N,double *X, double *Y, double *S){
		int j,fj,fi,k,m,el;

        	fi=i*BLOCK_SIZE*N;
                for (j = 0; j < num; j++) {
                	fj=j*BLOCK_SIZE;
                        for (k = 0; k < BLOCK_SIZE; k++) {
                      		el=fi+k*N+fj;
                                //for (m=0; m < BLOCK_SIZE; m++)
                                //      S[el+m] = X[el+m] + Y[el+m];

                                //
                                // loop unrolling the above for-loop
                                // IMPORTANT: if BLOCK_SIZE!=50 you need to do the appropriate changes

				S[el] = X[el] + Y[el];
				S[el+1] = X[el+1] + Y[el+1];
                                S[el+2] = X[el+2] + Y[el+2];
                                S[el+3] = X[el+3] + Y[el+3];
                                S[el+4] = X[el+4] + Y[el+4];
                                S[el+5] = X[el+5] + Y[el+5];
                                S[el+6] = X[el+6] + Y[el+6];
                                S[el+7] = X[el+7] + Y[el+7];
                                S[el+8] = X[el+8] + Y[el+8];
                                S[el+9] = X[el+9] + Y[el+9];
                                S[el+10] = X[el+10] + Y[el+10];
                                S[el+11] = X[el+11] + Y[el+11];
                                S[el+12] = X[el+12] + Y[el+12];
                                S[el+13] = X[el+13] + Y[el+13];
                                S[el+14] = X[el+14] + Y[el+14];
                                S[el+15] = X[el+15] + Y[el+15];
                                S[el+16] = X[el+16] + Y[el+16];
                                S[el+17] = X[el+17] + Y[el+17];
                                S[el+18] = X[el+18] + Y[el+18];
                                S[el+19] = X[el+19] + Y[el+19];
                                S[el+20] = X[el+20] + Y[el+20];
                                S[el+21] = X[el+21] + Y[el+21];
                                S[el+22] = X[el+22] + Y[el+22];
                                S[el+23] = X[el+23] + Y[el+23];
                                S[el+24] = X[el+24] + Y[el+24];
                                S[el+25] = X[el+25] + Y[el+25];
                                S[el+26] = X[el+26] + Y[el+26];
                                S[el+27] = X[el+27] + Y[el+27];
                                S[el+28] = X[el+28] + Y[el+28];
                                S[el+29] = X[el+29] + Y[el+29];
                                S[el+30] = X[el+30] + Y[el+30];
                                S[el+31] = X[el+31] + Y[el+31];
                                S[el+32] = X[el+32] + Y[el+32];
                                S[el+33] = X[el+33] + Y[el+33];
                                S[el+34] = X[el+34] + Y[el+34];
                                S[el+35] = X[el+35] + Y[el+35];
                                S[el+36] = X[el+36] + Y[el+36];
                                S[el+37] = X[el+37] + Y[el+37];
                                S[el+38] = X[el+38] + Y[el+38];
                                S[el+39] = X[el+39] + Y[el+39];
                                S[el+40] = X[el+40] + Y[el+40];
                                S[el+41] = X[el+41] + Y[el+41];
                                S[el+42] = X[el+42] + Y[el+42];
                                S[el+43] = X[el+43] + Y[el+43];
                                S[el+44] = X[el+44] + Y[el+44];
                                S[el+45] = X[el+45] + Y[el+45];
                                S[el+46] = X[el+46] + Y[el+46];
                                S[el+47] = X[el+47] + Y[el+47];
                                S[el+48] = X[el+48] + Y[el+48];
                                S[el+49] = X[el+49] + Y[el+49];
                         }
                }
	}

        void core_sub(int i,int num,long N,double *X, double *Y, double *S){
                int j,fj,fi,k,m,el;

                fi=i*BLOCK_SIZE*N;
                for (j = 0; j < num; j++) {
                        fj=j*BLOCK_SIZE;
                        for (k = 0; k < BLOCK_SIZE; k++) {
                                el=fi+k*N+fj;
                                //for (m=0; m < BLOCK_SIZE; m++)
                                //      S[el+m] = X[el+m] + Y[el+m];

                                //
                                // loop unrolling the above for-loop
                                // IMPORTANT: if BLOCK_SIZE!=50 you need to do the appropriate changes

                                S[el] = X[el] - Y[el];
                                S[el+1] = X[el+1] - Y[el+1];
                                S[el+2] = X[el+2] - Y[el+2];
                                S[el+3] = X[el+3] - Y[el+3];
                                S[el+4] = X[el+4] - Y[el+4];
                                S[el+5] = X[el+5] - Y[el+5];
                                S[el+6] = X[el+6] - Y[el+6];
                                S[el+7] = X[el+7] - Y[el+7];
                                S[el+8] = X[el+8] - Y[el+8];
                                S[el+9] = X[el+9] - Y[el+9];
                                S[el+10] = X[el+10] - Y[el+10];
                                S[el+11] = X[el+11] - Y[el+11];
                                S[el+12] = X[el+12] - Y[el+12];
                                S[el+13] = X[el+13] - Y[el+13];
                                S[el+14] = X[el+14] - Y[el+14];
                                S[el+15] = X[el+15] - Y[el+15];
                                S[el+16] = X[el+16] - Y[el+16];
                                S[el+17] = X[el+17] - Y[el+17];
                                S[el+18] = X[el+18] - Y[el+18];
                                S[el+19] = X[el+19] - Y[el+19];
                                S[el+20] = X[el+20] - Y[el+20];
                                S[el+21] = X[el+21] - Y[el+21];
                                S[el+22] = X[el+22] - Y[el+22];
                                S[el+23] = X[el+23] - Y[el+23];
                                S[el+24] = X[el+24] - Y[el+24];
                                S[el+25] = X[el+25] - Y[el+25];
                                S[el+26] = X[el+26] - Y[el+26];
                                S[el+27] = X[el+27] - Y[el+27];
                                S[el+28] = X[el+28] - Y[el+28];
                                S[el+29] = X[el+29] - Y[el+29];
                                S[el+30] = X[el+30] - Y[el+30];
                                S[el+31] = X[el+31] - Y[el+31];
                                S[el+32] = X[el+32] - Y[el+32];
                                S[el+33] = X[el+33] - Y[el+33];
                                S[el+34] = X[el+34] - Y[el+34];
                                S[el+35] = X[el+35] - Y[el+35];
                                S[el+36] = X[el+36] - Y[el+36];
                                S[el+37] = X[el+37] - Y[el+37];
                                S[el+38] = X[el+38] - Y[el+38];
                                S[el+39] = X[el+39] - Y[el+39];
                                S[el+40] = X[el+40] - Y[el+40];
                                S[el+41] = X[el+41] - Y[el+41];
                                S[el+42] = X[el+42] - Y[el+42];
                                S[el+43] = X[el+43] - Y[el+43];
                                S[el+44] = X[el+44] - Y[el+44];
                                S[el+45] = X[el+45] - Y[el+45];
                                S[el+46] = X[el+46] - Y[el+46];
                                S[el+47] = X[el+47] - Y[el+47];
                                S[el+48] = X[el+48] - Y[el+48];
                                S[el+49] = X[el+49] - Y[el+49];
                         }
                }
        }
               
	void madd(long N, double *X, double *Y, double *S) {
		// S=X+Y
		// X,Y square matrices of size N
		int i,j,k,m,fi,fj,fii,fjj,el;
		int num=N/BLOCK_SIZE;
                int tol1=N/NUM_THREADS;
                int thread;


		if(N%BLOCK_SIZE){
                	for(thread=0;thread<NUM_THREADS;thread++)
                        	#pragma omp task firstprivate(thread) private(i,j)
                        	{
                                	for (i=thread*tol1;i<(thread+1)*tol1;i++)
                                        	for(j=0;j<N;j++)
                                                	S[i*N + j] = X[i*N + j] + Y[i*N + j];
                        	}

                	#pragma omp taskwait

		
                	for(i=NUM_THREADS*tol1;i<N;i++)
                        	for(j=0;j<N;j++)
                                	 S[i*N + j] = X[i*N + j] + Y[i*N + j];
		}
		else{	
			tol1=num/NUM_THREADS;
			
			for(thread=0;thread<NUM_THREADS;thread++)
				#pragma omp task firstprivate(thread) private(i)
				{
                        		for (i=thread*tol1;i<(thread+1)*tol1;i++) {
						core_add(i,num,N,X,Y,S);
					}
				}
			
			#pragma omp taskwait

			for(i=NUM_THREADS*tol1;i<num;i++)
				core_add(i,num,N,X,Y,S);
		}

	}

	void msub(long N, double *X, double *Y, double *S) {
		// S=X-Y
                // X,Y square matrices of size N
                int i,j,k,m,fi,fj,fii,fjj,el;
                int num=N/BLOCK_SIZE;
                int tol1=N/NUM_THREADS;
                int thread;


                if(N%BLOCK_SIZE){
                        for(thread=0;thread<NUM_THREADS;thread++)
                                #pragma omp task firstprivate(thread) private(i,j)
                                {
                                        for (i=thread*tol1;i<(thread+1)*tol1;i++)
                                                for(j=0;j<N;j++)
                                                        S[i*N + j] = X[i*N + j] - Y[i*N + j];
                                }

                        #pragma omp taskwait


                        for(i=NUM_THREADS*tol1;i<N;i++)
                                for(j=0;j<N;j++)
                                         S[i*N + j] = X[i*N + j] - Y[i*N + j];
                }
                else{
                        tol1=num/NUM_THREADS;

                        for(thread=0;thread<NUM_THREADS;thread++)
                                #pragma omp task firstprivate(thread) private(i)
                                {
                                        for (i=thread*tol1;i<(thread+1)*tol1;i++) {
                                                core_sub(i,num,N,X,Y,S);
                                        }
                                }

                        #pragma omp taskwait

                        for(i=NUM_THREADS*tol1;i<num;i++)
                                core_sub(i,num,N,X,Y,S);
                } 

	
	}

	void matmul0(long N, double *X, double *Y, double *Z) {
		//   Z=X*Y strassen way, theoretical runtime is O(N^2.807) 
		//   X,Y matrices of size N 
		//   Matrices X and Y are split into four smaller
		//   (N/2)x(N/2) matrices as follows:
		//          _    _          _   _
		//     X = | A  B |    Y = | E F |
		//         | C  D |        | G H |
		//          -    -          -   -
		//   Then we build the following 7 matrices (requiring
		//   seven (N/2)x(N/2) matrix multiplications -- this is
		//   where the 2.807 = log2(7) improvement comes from):
		//     P0 = A*(F - H);
		//     P1 = (A + B)*H
		//     P2 = (C + D)*E
		//     P3 = D*(G - E);
		//     P4 = (A + D)*(E + H)
		//     P5 = (B - D)*(G + H)
		//     P6 = (A - C)*(E + F)
		//   The final result is
		//        _                                            _
		//   Z = | (P3 + P4) + (P5 - P1)   P0 + P1              |
		//       | P2 + P3                 (P0 + P4) - (P2 + P6)|
		//        -                                            -
		//  
  		if (N%2!=0 || N<=THRESHOLD_MULT) {
			if(N%BLOCK_SIZE==0)
				matmul1(N,X,Y,Z); // matrix-matrix multiplication  classic way enhanced with cache blocking,loop unrolling	
			else
   		 		matmul2(N, X, Y, Z); // matrix-matrix multiplication classic way
    			return;
  		}

		int n = N/2;      // size of sub-matrices
  		double *P[7];     // allocate temp matrices off heap
  		const int sz = n*n*sizeof(double);
  		int i,j;

  		double *A= (double *) malloc(sz);    // A-D matrices embedded in X
  		double *B= (double *) malloc(sz);
  		double *C= (double *) malloc(sz);
  		double *D= (double *) malloc(sz);

  		double *E= (double *) malloc(sz);    // E-H matrices embeded in Y
  		double *F= (double *) malloc(sz);
  		double *G= (double *) malloc(sz);
  		double *H= (double *) malloc(sz);

  		double *Z11 = (double *) malloc(sz);
  		double *Z12 = (double *) malloc(sz);
  		double *Z21 = (double *) malloc(sz);
  		double *Z22 = (double *) malloc(sz);


                int tol1=n/NUM_THREADS;
                int tol2=n%NUM_THREADS;
                int thread;

                for(thread=0;thread<NUM_THREADS;thread++)
                        #pragma omp task firstprivate(thread) private(i,j)
                        {
                                for (i=thread*tol1;i<(thread+1)*tol1;i++){
                                        for(j=0;j<n;j++){
        		                        A[i*n+j]=X[i*N+j];
	                	                B[i*n+j]=X[i*N+j+n];
                                		C[i*n+j]=X[(i+n)*N+j];
                                		D[i*n+j]=X[(i+n)*N+j+n];

                                		E[i*n+j]=Y[i*N+j];
                                		F[i*n+j]=Y[i*N+j+n];
                                		G[i*n+j]=Y[(i+n)*N+j];
                                		H[i*n+j]=Y[(i+n)*N+j+n];
					}
				}

                        }

                #pragma omp taskwait

		
                for(i=NUM_THREADS*tol1;i<n;i++){
                        for(j=0;j<n;j++){
                                A[i*n+j]=X[i*N+j];
                                B[i*n+j]=X[i*N+j+n];
                                C[i*n+j]=X[(i+n)*N+j];
                                D[i*n+j]=X[(i+n)*N+j+n];

                                E[i*n+j]=Y[i*N+j];
                                F[i*n+j]=Y[i*N+j+n];
                                G[i*n+j]=Y[(i+n)*N+j];
                                H[i*n+j]=Y[(i+n)*N+j+n];
			}
		}
                                 

  		for ( i = 0; i < 7; i++)
    			P[i] = (double *) malloc(sz);

		#pragma omp task 
		{
  			// P0 = A*(F - H);
  			double *T = (double *) malloc(sz);
  			msub(n, F, H, T);
  			matmul0(n, A, T, P[0]);
			free(T);
		}
		#pragma omp task 
		{   
  			// P1 = (A + B)*H
  			double *T = (double *) malloc(sz);
  			madd(n, A, B, T);
  			matmul0(n, T, H, P[1]);
			free(T);
		}
		#pragma omp task
		{
  			// P2 = (C + D)*E
			double *T = (double *) malloc(sz);
  			madd(n, C, D, T);
  			matmul0(n, T, E, P[2]);
			free(T);
		}
		#pragma omp  task 
		{
  			// P3 = D*(G - E);
			double *T = (double *) malloc(sz);
  			msub(n, G, E, T);
  			matmul0(n, D, T, P[3]);
			free(T);
		}
		#pragma omp  task 
		{
  			// P4 = (A + D)*(E + H)
			double *T = (double *) malloc(sz);
                        double *U = (double *) malloc(sz);
			madd(n, A, D, T);
  			madd(n, E, H, U);
  			matmul0(n, T, U, P[4]);
			free(T);
			free(U);
		}
		#pragma omp  task 
		{
  			// P5 = (B - D)*(G + H)
			double *T = (double *) malloc(sz);
                        double *U = (double *) malloc(sz);
  			msub(n, B, D, T);
  			madd(n, G, H, U);
  			matmul0(n, T, U, P[5]);
			free(T);
			free(U);
		}
		#pragma omp  task 
		{
  			// P6 = (A - C)*(E + F)
			double *T = (double *) malloc(sz);
                        double *U = (double *) malloc(sz);
  			msub(n, A, C, T);
  			madd(n, E, F, U);
  			matmul0(n, T, U, P[6]);
			free(T);
			free(U);
		}

		#pragma omp taskwait

		#pragma omp task
		{
  			// Z upper left = (P3 + P4) + (P5 - P1)
                        double *T = (double *) malloc(sz);
                        double *U = (double *) malloc(sz);
  			madd(n, P[4], P[3], T);
  			msub(n, P[5], P[1], U);
  			madd(n, T, U, Z11);
			free(T);
			free(U);
		}

		#pragma omp task 
		{
  			// Z lower left = P2 + P3
  			madd(n, P[2], P[3], Z21);

  			// Z upper right = P0 + P1
  			madd(n, P[0], P[1], Z12);
		}

		#pragma omp task
		{
  			// Z lower right = (P0 + P4) - (P2 + P6)
                        double *T = (double *) malloc(sz);
                        double *U = (double *) malloc(sz);
  			madd(n, P[0], P[4], T);
  			madd(n, P[2], P[6], U);
  			msub(n, T, U, Z22);
			free(T);
			free(U);
		}

		#pragma omp taskwait

                tol1=n/NUM_THREADS;
              
                for(thread=0;thread<NUM_THREADS;thread++)
                        #pragma omp task firstprivate(thread) private(i,j)
                        {
                                for (i=thread*tol1;i<(thread+1)*tol1;i++){
		                        for(j=0;j<n;j++){
                		                Z[i*N+j]=Z11[i*n+j];
                                		Z[i*N+j+n]=Z12[i*n+j];
                                		Z[(i+n)*N+j]=Z21[i*n+j];
                                		Z[(i+n)*N+j+n]=Z22[i*n+j];
                        		}
				}
                        }

                #pragma omp taskwait


                for(i=NUM_THREADS*tol1;i<n;i++){
                        for(j=0;j<n;j++){
                                Z[i*N+j]=Z11[i*n+j];
                                Z[i*N+j+n]=Z12[i*n+j];
                                Z[(i+n)*N+j]=Z21[i*n+j];
                                Z[(i+n)*N+j+n]=Z22[i*n+j];
                        }
		}


  		// deallocate temp matrices
  		free(A);
  		free(B);
  		free(C);
  		free(D);
  		free(E);
  		free(F);
  		free(G);
  		free(H); 
  		for (i = 6; i >= 0; i--)
    			free(P[i]);

	}

/*******************************************************
*                   function's code                    *
*******************************************************/
        #pragma omp parallel num_threads(NUM_THREADS)
        {
               	#pragma omp single
               	{
			if(N<=THRESHOLD_MATRIX)
				THRESHOLD_MULT=180;
			else
				THRESHOLD_MULT=1200;

			matmul0(N, a, b, c); 
               	}
        }
}


int main(int argc, char **argv) {
   long N;
   double *A, *B, *C, t;

   N = atoi(argv[1]);
   A = (double *) malloc(N * N * sizeof(double));
   B = (double *) malloc(N * N * sizeof(double));
   C = (double *) malloc(N * N * sizeof(double));
   matfill(N, A, 1.0);
   matfill(N, B, 2.0);
   matfill(N, C, 0.0);

   t = gettime();
   matmul(N, A, B, C);
   t = gettime() - t;

   fprintf(stdout, "%ld\t%le\t%le\n", N, t, (2 * N - 1) * N * N / t);
   fflush(stdout);

   free(A);
   free(B);
   free(C);

   return EXIT_SUCCESS;
}