GEMM_chunked_wrong.jdf

extern "C" %{

/**
 * This second example shows how to create a simple jdf that has only one single task.
 *    JDF syntax
 *    parsec_JDFNAME_New()
 *    parsec_context_add_taskpool()
 *    parsec_data_collection_init()
 *
 * Can play with the HelloWorld bounds to show embarissingly parallel algorithm.
 *
 * @version 3.0
 * @email parsec-users@icl.utk.edu
 *
 */

/* Paste code to allocate a matrix in desc if cond_init is true */
#define PASTE_CODE_ALLOCATE_MATRIX(DC, COND, TYPE, INIT_PARAMS)      \
    TYPE##_t DC;                                                     \
    if(COND) {                                                          \
        TYPE##_init INIT_PARAMS;                                        \
        DC.mat = parsec_data_allocate((size_t)DC.super.nb_local_tiles * \
                                        (size_t)DC.super.bsiz *      \
                                        (size_t)parsec_datadist_getsizeoftype(DC.super.mtype)); \
        parsec_data_collection_set_key((parsec_data_collection_t*)&DC, #DC);          \
    }

#include <math.h>
#include "parsec.h"
#include "parsec/data_dist/matrix/matrix.h"
#include "parsec/data_dist/matrix/two_dim_rectangle_cyclic.h"

#define CHUNK_SIZE_X 2
#define CHUNK_SIZE_Y 2
#define CHUNK_SIZE_Z 2

#define MATRIX_M 4 // M rows
#define MATRIX_N 4 // N columns
#define MATRIX_K 8 // depth
#define TILE_SIZE 8

#define CHUNK_NUMBER_X (MATRIX_M/CHUNK_SIZE_X)
#define CHUNK_NUMBER_Y (MATRIX_N/CHUNK_SIZE_Y)
#define CHUNK_NUMBER_Z (MATRIX_K/CHUNK_SIZE_Z)

%}


descA  [ type="parsec_matrix_block_cyclic_t*" ]
descB  [ type="parsec_matrix_block_cyclic_t*" ]
descC  [ type="parsec_matrix_block_cyclic_t*" ]


rank   [ type="int" ]
nodes   [ type="int" ]
matrix_size_m   [ type="int" ] // row number
matrix_size_n   [ type="int" ] // column number
matrix_size_k   [ type="int" ]
chunk_number_x   [ type="int" ]
chunk_number_y   [ type="int" ]
chunk_number_z   [ type="int" ]
chunk_size_x   [ type="int" ]
chunk_size_y   [ type="int" ]
chunk_size_z   [ type="int" ]


FillA(m, k)

m = 0 .. matrix_size_m-1
k = 0 .. matrix_size_k-1

: descA(m, k)

RW A <- descA(m, k)
    -> A GEMM(m, 0 .. matrix_size_n-1, k, m/chunk_number_x, 0 .. chunk_number_y-1, 0)

BODY
    //int rank = this_task->parsec_context->my_rank;

    double *matA = A;

printf("process %d fills matrix A at coordinates %d, %d\n", rank, m, k);

    for(int i=0;i<TILE_SIZE*TILE_SIZE;++i)
    {
        //matA[i] = (double)((151+i+i*i*8)%256);
        matA[i] = 1;
    }
END

FillB(k, n)

n = 0 .. matrix_size_n-1
k = 0 .. matrix_size_k-1

: descB(k, n)

RW B <- descB(k, n)
    -> B GEMM(0 .. matrix_size_m-1, n, k, 0 .. chunk_number_x-1, n/chunk_number_y, 0)

BODY
    double *matB = B;

    for(int i=0;i<TILE_SIZE*TILE_SIZE;++i)
    {
        //matB[i] = (double)((151+i+i*i*8)%256);
        matB[i] = 1;
    }
END

/*
Sync(x, y, z)

x = 0 .. chunk_number_x-1
y = 0 .. chunk_number_y-1
z = 0 .. chunk_number_z-1

// TODO check
: descC(x, y)

CTL X <- X GEMM(
        x*chunk_size_x .. (x+1)*chunk_size_x-1,
        y*chunk_size_y .. (y+1)*chunk_size_y-1,
        (z+1)*chunk_size_z-1,
        x, y, z)
    -> X GEMM(
        x*chunk_size_x .. (x+1)*chunk_size_x-1,
        y*chunk_size_y .. (y+1)*chunk_size_y-1,
        (z+1)*chunk_size_z,
        x, y, z+1)

BODY
    //printf("Syncing ... (%d, %d, %d) (depends on (%d-%d, %d-%d, %d))\n", x, y, z,
    //        x*chunk_size_x, (x+1)*chunk_size_x-1,
    //        y*chunk_size_y, (y+1)*chunk_size_y-1,
    //        (z+1)*chunk_size_z-1);
    printf("Syncing ... (%d, %d, %d) (will liberate GEMM(%d-%d, %d-%d, %d, %d, %d, %d))\n", x, y, z,
            x*chunk_size_x, (x+1)*chunk_size_x-1,
            y*chunk_size_y, (y+1)*chunk_size_y-1,
            (z+1)*chunk_size_z,
            x, y, z+1);
END
*/

GEMM(m, n, k, x, y, z)

m = 0 .. matrix_size_m-1
n = 0 .. matrix_size_n-1
k = 0 .. matrix_size_k-1

x = 0 .. chunk_number_x-1
y = 0 .. chunk_number_y-1
z = 0 .. chunk_number_z-1

: descC(m, n)

/*
READ A <- A FillA(m, k)
READ B <- B FillB(k, n)
RW C <- (k==0) ? descC(m, n) : C GEMM(m, n, k-1)
    -> (k<matrix_size_k-1) ? C GEMM(m, n, k+1) : C PrintC(m, n)


// old

RW C <- (k==0) ? descC(m, n) : C GEMM(m, n, k-1, x, y, (k%chunk_size_z==0)?-1:z)
    -> (k<matrix_size_k-1) ? C GEMM(m, n, k+1, x, y, (k%chunk_size_z==chunk_size_z-1)?-1:z)// : C PrintC(m, n)
*/


READ A <- A FillA(m, k)
READ B <- B FillB(k, n)
/*
RW C <- (k==0) ? descC(m, n) : C GEMM(m, n, k-1, x, y, (k-1)/chunk_size_z)
    -> (k<matrix_size_k-1) ? C GEMM(m, n, k+1, x, y, (k+1)/chunk_size_z)
*/
RW C <- (k==0) ? descC(m, n) : C GEMM(m, n, k-1, x, y, z)
    -> (k<matrix_size_k-1) ? C GEMM(m, n, k+1, x, y, z) : C PrintC(m, n)

/*
CTL X <- (k%chunk_size_z==0 && k != 0) ? X Sync(x, y, z-1)
    -> (k%chunk_size_z==chunk_size_z-1 && k != matrix_size_k-1) ? X Sync(x, y, z)
*/

BODY
    //printf("[GEMM] process %d is at coordinates %d, %d, (k=%d), chunk=(%d, %d, %d)\n", rank, m, n, k, x, y, z);
    printf("[GEMM] GEMM(%d, %d, %d, %d, %d, %d) will liberate GEMM(%d, %d, %d, %d, %d, %d)\n",
            m, n, k, x, y, z,
            m, n, k+1, x, y, (k+1)/chunk_size_z);

    double *matA = A;
    double *matB = B;
    double *matC = C;

    for(int i=0;i<TILE_SIZE;++i)
    {
        for(int j=0;j<TILE_SIZE;++j)
        {
            for(int k=0;k<TILE_SIZE;++k)
            {
                matC[i+j*TILE_SIZE] += matA[i+k*TILE_SIZE]*matB[k+j*TILE_SIZE];
            }
        }
    }
END


PrintC(m, n)

m = 0 .. matrix_size_m-1
n = 0 .. matrix_size_n-1

: descC(m, n)

READ C <- C GEMM(m, n, matrix_size_k-1, m/chunk_number_x, n/chunk_number_y, chunk_number_z-1)

BODY
    double *matC = C;

    printf("process %d matC %d, %d = %f\n", rank, m, n, matC[0]);
END


extern "C" %{

int main(int argc, char *argv[])
{
    parsec_context_t* parsec;
    int rc;
    int rank, world;
    parsec_test_taskpool_t *tp;
    //int mycounter;

    parsec_arena_datatype_t adt;
    parsec_datatype_t otype;

#if defined(PARSEC_HAVE_MPI)
    {
        int provided;
        MPI_Init_thread(&argc, &argv, MPI_THREAD_SERIALIZED, &provided);
    }
    MPI_Comm_size(MPI_COMM_WORLD, &world);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
#else
    world = 1;
    rank = 0;
#endif

    parsec = parsec_init(-1, &argc, &argv);

    //mycounter = 300 + rank;
    int test_data = 0;

    int nodes = world;

    //int KP = sqrt(nodes);
    //int KQ = nodes/KP;
    int KP = 2;
    int KQ = 2;

    PASTE_CODE_ALLOCATE_MATRIX(descA, true, parsec_matrix_block_cyclic, (
                               &descA,
                               PARSEC_MATRIX_DOUBLE,
                               PARSEC_MATRIX_TILE,
                               rank,
                               TILE_SIZE,      TILE_SIZE,   /* Tile size */
                               MATRIX_M*TILE_SIZE,    MATRIX_K*TILE_SIZE,   /* Global matrix size (what is stored)*/
                               0,           0,    /* Staring point in the global matrix */
                               MATRIX_M,    MATRIX_K,    /* Submatrix size (the one concerned by the computation */
                               KP,          KQ,    /* process process grid*/
                               1,           1,   /* k-cyclicity */
                               0,           0)   /* starting point on the process grid*/
                            );

    PASTE_CODE_ALLOCATE_MATRIX(descB, true, parsec_matrix_block_cyclic, (
                               &descB,
                               PARSEC_MATRIX_DOUBLE,
                               PARSEC_MATRIX_TILE,
                               rank,
                               TILE_SIZE,      TILE_SIZE,   /* Tile size */
                               MATRIX_K*TILE_SIZE,    MATRIX_N*TILE_SIZE,   /* Global matrix size (what is stored)*/
                               0,           0,    /* Staring point in the global matrix */
                               MATRIX_K,    MATRIX_N,    /* Submatrix size (the one concerned by the computation */
                               KP,          KQ,    /* process process grid*/
                               1,           1,   /* k-cyclicity */
                               0,           0)   /* starting point on the process grid*/
                            );

    PASTE_CODE_ALLOCATE_MATRIX(descC, true, parsec_matrix_block_cyclic, (
                               &descC,
                               PARSEC_MATRIX_DOUBLE,
                               PARSEC_MATRIX_TILE,
                               rank,
                               TILE_SIZE,      TILE_SIZE,   /* Tile size */
                               MATRIX_M*TILE_SIZE,    MATRIX_N*TILE_SIZE,   /* Global matrix size (what is stored)*/
                               0,           0,    /* Staring point in the global matrix */
                               MATRIX_M,    MATRIX_N,    /* Submatrix size (the one concerned by the computation */
                               KP,          KQ,    /* process process grid*/
                               1,           1,   /* k-cyclicity */
                               0,           0)   /* starting point on the process grid*/
                            );

    parsec_translate_matrix_type(PARSEC_MATRIX_DOUBLE, &otype);
    parsec_add2arena_rect(&adt, otype,
                                 descC.super.mb, descC.super.nb, descC.super.mb);

    tp = (parsec_test_taskpool_t*)parsec_test_new(&descA, &descB, &descC, rank, world, MATRIX_M, MATRIX_N, MATRIX_K, CHUNK_NUMBER_X, CHUNK_NUMBER_Y, CHUNK_NUMBER_Z, CHUNK_SIZE_X, CHUNK_SIZE_Y, CHUNK_SIZE_Z);

    assert( NULL != tp );
    tp->arenas_datatypes[PARSEC_test_DEFAULT_ADT_IDX] = adt;
    PARSEC_OBJ_RETAIN(adt.arena);

    rc = parsec_context_add_taskpool( parsec, (parsec_taskpool_t*)tp );
    PARSEC_CHECK_ERROR(rc, "parsec_context_add_taskpool");
    rc = parsec_context_start(parsec);
    PARSEC_CHECK_ERROR(rc, "parsec_context_start");
    rc = parsec_context_wait(parsec);
    PARSEC_CHECK_ERROR(rc, "parsec_context_wait");

printf("final : %f\n", ((double*)(descA.mat))[0]);

    parsec_taskpool_free(&tp->super);

    parsec_data_free(descA.mat);
    parsec_tiled_matrix_destroy( (parsec_tiled_matrix_t*)&descA);

    parsec_data_free(descB.mat);
    parsec_tiled_matrix_destroy( (parsec_tiled_matrix_t*)&descB);

    parsec_data_free(descC.mat);
    parsec_tiled_matrix_destroy( (parsec_tiled_matrix_t*)&descC);

    parsec_fini(&parsec);
#if defined(PARSEC_HAVE_MPI)
    MPI_Finalize();
#endif

    return 0;
}

%}