const int blockSize = 12;

void readBlock(double outputMatrix[], double inputMatrix[], int beginX, int endX, int beginY, int endY, int n) {

        int i = 0;
        for (int x = beginX; x < endX; x++) {
                for (int y = beginY; y < endY; y++) {
                        outputMatrix[i] = inputMatrix[x*n + y];
                        i++;
                }
                for (int y = endY-beginY; y < blockSize; y++) {
                        outputMatrix[i] = 0;
                        i++;
                }


        }

        for (int x = endX - beginX; x < blockSize; x++) {
                for (int y = 0; y < blockSize; y++) {
                        outputMatrix[i] = 0;
                        i++;
                }

        }

}


int min(int a, int b) {
        if (a < b) return a;
        else return b;

}

void multiplyMatrices(double inputMatrix1[], double inputMatrix2[], double outputMatrix[], int endx, int endy) {
        for (int i = 0; i < endx; i++) {
                for (int j = 0; j < endy; j++) {
                        double cvalue = outputMatrix[i*blockSize + j];
                        for (int l = 0; l < blockSize; l+=2) {
                                double a1 = inputMatrix1[l*blockSize + j] * inputMatrix2[i*blockSize + l];;
                                double a2 = inputMatrix1[(l+1)*blockSize + j] * inputMatrix2[i*blockSize + l + 1];
                                cvalue += a1 + a2;

                        }
                        outputMatrix[i*blockSize + j] = cvalue;
                }
        }

}


void writeMatrix(double* inputMatrix, double* outputMatrix, int beginX, int endX, int beginY, int endY, int n) {
        int k,l=0;
        for (int i = beginX; i < endX; i++) {
                k = 0;
                for (int j = beginY; j < endY; j++) {
                        outputMatrix[i*n + j] = inputMatrix[l*blockSize+k];
                        k++;
                }
                l++;
        }

}


void square_dgemm(int n, double* A, double* B, double* C) {



        for (int x = 0; x < n; x += blockSize) {
                for (int y = 0; y < n; y += blockSize) {
                        double CC[blockSize*blockSize];

                        readBlock(CC, C, x, min(x + blockSize, n), y, min(y + blockSize, n), n);

                        for (int z = 0; z < n; z += blockSize) {
                                double BC[blockSize*blockSize];
                                double AC[blockSize*blockSize];
                                int zAdj = min(z + blockSize, n);
                                readBlock(AC, A, z, zAdj,y, min(y+blockSize, n), n);

                                readBlock(BC, B, x, min(x + blockSize, n), z,zAdj, n);

                                multiplyMatrices(AC, BC, CC, min(x+blockSize,n)-x,min(y+blockSize,n)-y);

                        }
                        writeMatrix(CC, C, x, min(x + blockSize, n), y, min(y + blockSize, n), n);




                }

        }
}