#include "stdio.h"

////////////////////////////////////
// Image parameters

#define PIXEL_SIZE	2
#define NL		1024
#define NP		1024
#define BLOCK_SIZE	1024

#define PRINTF		if(lid==0) tty_printf

#define TA(c,l,p)  (A[c][((NP)*(l))+(p)])
#define TB(c,p,l)  (B[c][((NL)*(p))+(l)])
#define TC(c,l,p)  (C[c][((NP)*(l))+(p)])
#define TD(c,l,p)  (D[c][((NP)*(l))+(p)])
#define TZ(c,l,p)  (Z[c][((NP)*(l))+(p)])

#define max(x,y) ((x) > (y) ? (x) : (y))
#define min(x,y) ((x) < (y) ? (x) : (y))

///////////////////////////////////////////
// tricks to read parameters from ldscript
///////////////////////////////////////////

struct plaf;

extern struct plaf seg_heap_base;
extern struct plaf NB_PROCS;
extern struct plaf NB_CLUSTERS;

/////////////
void main()
{

//////////////////////////////////
// convolution kernel parameters
// The content of this section is
// Philips proprietary information.
///////////////////////////////////

    int	vnorm  = 115;
    int	vf[35];
    vf[0]  = 1;
    vf[1]  = 1;
    vf[2]  = 2;
    vf[3]  = 2;
    vf[4]  = 2;
    vf[5]  = 2;
    vf[6]  = 3;
    vf[7]  = 3;
    vf[8]  = 3;
    vf[9]  = 4;
    vf[10] = 4;
    vf[11] = 4;
    vf[12] = 4;
    vf[13] = 5;
    vf[14] = 5;
    vf[15] = 5;
    vf[16] = 5;
    vf[17] = 5;
    vf[18] = 5;
    vf[19] = 5;
    vf[20] = 5;
    vf[21] = 5;
    vf[22] = 4;
    vf[23] = 4;
    vf[24] = 4;
    vf[25] = 4;
    vf[26] = 3;
    vf[27] = 3;
    vf[28] = 3;
    vf[29] = 2;
    vf[30] = 2;
    vf[31] = 2;
    vf[32] = 2;
    vf[33] = 1;
    vf[34] = 1;

    int hrange = 100;
    int hnorm  = 201;

    unsigned int date      = 0;

    int c;                                              	// cluster index for loops
    int l;                                              	// line index for loops
    int p;                                              	// pixel index for loops
    int x;                                              	// filter index for loops

    int pid                 = procid();                         // processor id
    int nprocs              = (int)&NB_PROCS;          		// number of processors per cluster
    int nclusters           = (int)&NB_CLUSTERS;       		// number of clusters
    int lid                 = pid%nprocs;                       // local task id
    int cid                 = pid/nprocs;                       // cluster task id
    int base                = (unsigned int)&seg_heap_base;     // base address for shared buffers
    int increment           = (0x80000000 / nclusters) * 2;     // cluster increment
    int ntasks              = nclusters * nprocs;               // number of tasks
    int nblocks             = (NP*NL*PIXEL_SIZE)/BLOCK_SIZE;  	// number of blocks per image 

    int lines_per_task      = NL/ntasks;			// number of lines per task
    int lines_per_cluster   = NL/nclusters;			// number of lines per cluster
    int pixels_per_task     = NP/ntasks;			// number of columns per task
    int pixels_per_cluster  = NP/nclusters;			// number of columns per cluster

    int first, last;

    PRINTF("\n*** Processor %d entering main at cycle %d ***\n\n", pid, proctime());
    
    //////////////////////////
    //  parameters checking
    if( (nprocs != 1) && (nprocs != 2) && (nprocs != 4) )
    {
        PRINTF("NB_PROCS must be 1, 2 or 4\n");
        while(1);
    }
    if( (nclusters !=  4) && (nclusters !=  8) && (nclusters != 16) && 
        (nclusters != 32) && (nclusters != 64) && (nclusters !=128) && (nclusters != 256) )
    {
        PRINTF("NB_CLUSTERS must be a power of 2 between 4 and 256\n");
        while(1);
    }
    if( pid >= ntasks )
    {
        PRINTF("processor id %d larger than NB_CLUSTERS*NB_PROCS\n", pid);
        while(1);
    }
    if ( NL % nclusters != 0 )
    {
        PRINTF("NB_CLUSTERS must be a divider of NL");
        while(1);
    }
    if( NP % nclusters != 0 )
    {
        PRINTF("NB_CLUSTERS must be a divider of NP");
        while(1);
    }

    //////////////////////////////////////////////////////////////////
    // Arrays of pointers on the shared, distributed buffers  
    // containing the images (sized for the worst case : 256 clusters)
    unsigned short*	A[256];
    int*		B[256];
    int*		C[256];
    int*		D[256];
    unsigned char*	Z[256];
    
    // Arrays of pointers on the instrumentation arrays
    // These arrays are indexed by the cluster index (sized for the worst case : 256 clusters)
    // each pointer points on the base adress of an array of 4 (NPROCS max) unsigned int
    unsigned int*       LOAD_START[256];
    unsigned int*       LOAD_ENDED[256];
    unsigned int*       VERT_START[256];
    unsigned int*       VERT_ENDED[256];
    unsigned int*       HORI_START[256];
    unsigned int*       HORI_ENDED[256];
    unsigned int*       DISP_START[256];
    unsigned int*       DISP_ENDED[256];

    // The shared, distributed buffers addresses are computed
    // from the seg_heap_base value defined in the ldscript file
    // and from the cluster increment = 4Gbytes/nclusters.
    // These arrays of pointers are identical and
    // replicated in the stack of each task 
    for( c=0 ; c<nclusters ; c++)
    {
        A[c] = (unsigned short*)	(base 				+ increment*c);
        Z[c] = (unsigned char*)		(base + 2*NP*NL/nclusters 	+ increment*c);
        B[c] = (int*)			(base + 4*NP*NL/nclusters 	+ increment*c);
        C[c] = (int*)			(base + 8*NP*NL/nclusters 	+ increment*c);
        D[c] = (int*)			(base + 12*NP*NL/nclusters 	+ increment*c);

        LOAD_START[c] = (unsigned int*) (base + 3*NL*NP/nclusters 	+ increment*c);
        LOAD_ENDED[c] = (unsigned int*) (base + 3*NL*NP/nclusters + 16	+ increment*c);
        VERT_START[c] = (unsigned int*) (base + 3*NL*NP/nclusters + 32	+ increment*c);
        VERT_ENDED[c] = (unsigned int*) (base + 3*NL*NP/nclusters + 48	+ increment*c);
        HORI_START[c] = (unsigned int*) (base + 3*NL*NP/nclusters + 64	+ increment*c);
        HORI_ENDED[c] = (unsigned int*) (base + 3*NL*NP/nclusters + 80	+ increment*c);
        DISP_START[c] = (unsigned int*) (base + 3*NL*NP/nclusters + 96	+ increment*c);
        DISP_ENDED[c] = (unsigned int*) (base + 3*NL*NP/nclusters + 112	+ increment*c);
    }

    PRINTF("NCLUSTERS = %d\n", nclusters); 
    PRINTF("NPROCS    = %d\n\n", nprocs); 

    PRINTF("*** Starting barrier init at cycle %d ***\n", proctime());

    //  barriers initialization
    barrier_init(0, ntasks);
    barrier_init(1, ntasks);
    barrier_init(2, ntasks);
    barrier_init(3, ntasks);

    PRINTF("*** Completing barrier init at cycle %d ***\n", proctime());

    ////////////////////////////////////////////////////////
    // pseudo parallel load from disk to A[c] buffers
    // only task running on processor with (lid==0) does it
    // nblocks/nclusters are loaded in each cluster

    if ( lid == 0 )
    {
        int p;
        date  = proctime();
        PRINTF("\n*** Starting load at cycle %d\n", date);
        for ( p=0 ; p<nprocs ; p++ ) LOAD_START[cid][p] = date;

        if( ioc_read(nblocks*cid/nclusters, 
                     A[cid] , 
                     nblocks/nclusters) )
        {
            PRINTF("echec ioc_read\n");
            while(1);
        }
        if ( ioc_completed() )
        {
            PRINTF("echec ioc_completed\n");
            while(1);
        }

        date  = proctime();
        PRINTF("*** Completing load at cycle %d\n", date);
        for ( p=0 ; p<nprocs ; p++ ) LOAD_ENDED[cid][p] = date;
    }

    barrier_wait(0);

    //////////////////////////////////////////////////////////
    // parallel horizontal filter : 
    // B <= transpose(FH(A))
    // D <= A - FH(A)
    // Each task computes (NL/ntasks) lines 
    // The image must be extended :
    // if (z<0) 	TA(cid,l,z) == TA(cid,l,0)
    // if (z>NP-1)	TA(cid,l,z) == TA(cid,l,NP-1)

    date  = proctime();
    PRINTF("\n*** Starting horizontal filter at cycle %d\n", date);
    HORI_START[cid][lid] = date;

    // l = absolute line index / p = absolute pixel index  
    // first & last define which lines are handled by a given task(cid,lid)

    first = (cid*nprocs + lid)*lines_per_task;
    last  = first + lines_per_task;

    for ( l=first ; l<last ; l++)
    {
        // src_c and src_l are the cluster index and the line index for A & D
        int src_c = l/lines_per_cluster;
        int src_l = l%lines_per_cluster;

        // We use the spécific values of the horizontal ep-filter for optimisation:
        // sum(p) = sum(p-1) + TA[p+hrange] - TA[p-hrange-1]
        // To minimize the number of tests, the loop on pixels is split in three domains 

        int sum_p = (hrange+2)*TA(src_c, src_l, 0);
        for ( x = 1 ; x < hrange ; x++) sum_p = sum_p + TA(src_c, src_l, x);

        // first domain : from 0 to hrange
        for ( p=0 ; p<hrange+1 ; p++)
        {
            // dst_c and dst_p are the cluster index and the pixel index for B
            int dst_c = p/pixels_per_cluster;
            int dst_p = p%pixels_per_cluster;
            sum_p = sum_p + (int)TA(src_c, src_l, p+hrange) - (int)TA(src_c, src_l, 0);
            TB(dst_c, dst_p, l) = sum_p/hnorm;
            TD(src_c, src_l, p) = (int)TA(src_c, src_l, p) - sum_p/hnorm;
        }
        // second domain : from (hrange+1) to (NP-hrange-1)
        for ( p = hrange+1 ; p < NP-hrange ; p++)
        {
            // dst_c and dst_p are the cluster index and the pixel index for B
            int dst_c = p/pixels_per_cluster;
            int dst_p = p%pixels_per_cluster;
            sum_p = sum_p + (int)TA(src_c, src_l, p+hrange) - (int)TA(src_c, src_l, p-hrange-1);
            TB(dst_c, dst_p, l) = sum_p/hnorm;
            TD(src_c, src_l, p) = (int)TA(src_c, src_l, p) - sum_p/hnorm;
        }
        // third domain : from (NP-hrange) to (NP-1)
        for ( p = NP-hrange ; p < NP ; p++)
        {
            // dst_c and dst_p are the cluster index and the pixel index for B
            int dst_c = p/pixels_per_cluster;
            int dst_p = p%pixels_per_cluster;
            sum_p = sum_p + (int)TA(src_c, src_l, NP-1) - (int)TA(src_c, src_l, p-hrange-1);
            TB(dst_c, dst_p, l) = sum_p/hnorm;
            TD(src_c, src_l, p) = (int)TA(src_c, src_l, p) - sum_p/hnorm;
        }

        PRINTF(" - line %d computed at cycle %d\n", l, proctime());
    }

    date  = proctime();
    PRINTF("*** Completing horizontal filter at cycle %d\n", date);
    HORI_ENDED[cid][lid] = date;

    barrier_wait(1);

    //////////////////////////////////////////////////////////
    // parallel vertical filter : 
    // C <= transpose(FV(B))
    // Each task computes (NP/ntasks) columns
    // The image must be extended :
    // if (l<0) 	TB(cid,p,x) == TB(cid,p,0)
    // if (l>NL-1)	TB(cid,p,x) == TB(cid,p,NL-1)

    date  = proctime();
    PRINTF("\n*** starting vertical filter at cycle %d\n", date);
    VERT_START[cid][lid] = date;

    // l = absolute line index / p = absolute pixel index
    // first & last define which pixels are handled by a given task(cid,lid)

    first = (cid*nprocs + lid)*pixels_per_task;
    last  = first + pixels_per_task;

    for ( p=first ; p<last ; p++)
    {
        // src_c and src_p are the cluster index and the pixel index for B
        int src_c = p/pixels_per_cluster;
        int src_p = p%pixels_per_cluster;

        int sum_l;

        // We use the specific values of the vertical ep-filter
        // To minimize the number of tests, the NL lines are split in three domains 

        // first domain : explicit computation for the first 18 values
        for ( l=0 ; l<18 ; l++)
        {
            // dst_c and dst_l are the cluster index and the line index for C
            int dst_c = l/lines_per_cluster;
            int dst_l = l%lines_per_cluster;

            for ( x=0, sum_l=0 ; x<35 ; x++ )
            {
                sum_l = sum_l + vf[x] * TB(src_c, src_p, max(l-17+x,0) );
            }
            TC(dst_c, dst_l, p) = sum_l/vnorm;
        }
        // second domain
        for ( l = 18 ; l < NL-17 ; l++ )
        {
            // dst_c and dst_l are the cluster index and the line index for C
            int dst_c = l/lines_per_cluster;
            int dst_l = l%lines_per_cluster;

            sum_l = sum_l + TB(src_c, src_p, l+4)
                          + TB(src_c, src_p, l+8)
                          + TB(src_c, src_p, l+11)
                          + TB(src_c, src_p, l+15)
                          + TB(src_c, src_p, l+17)
                          - TB(src_c, src_p, l-5)
                          - TB(src_c, src_p, l-9)
                          - TB(src_c, src_p, l-12)
                          - TB(src_c, src_p, l-16)
                          - TB(src_c, src_p, l-18);
            TC(dst_c, dst_l, p) = sum_l/vnorm;
        }
        // third domain
        for ( l = NL-17 ; l < NL ; l++ )
        {
            // dst_c and dst_l are the cluster index and the line index for C
            int dst_c = l/lines_per_cluster;
            int dst_l = l%lines_per_cluster;

            sum_l = sum_l + TB(src_c, src_p, min(l+4,NL-1))
                          + TB(src_c, src_p, min(l+8,NL-1))
                          + TB(src_c, src_p, min(l+11,NL-1))
                          + TB(src_c, src_p, min(l+15,NL-1))
                          + TB(src_c, src_p, min(l+17,NL-1))
                          - TB(src_c, src_p, l-5)
                          - TB(src_c, src_p, l-9)
                          - TB(src_c, src_p, l-12)
                          - TB(src_c, src_p, l-16)
                          - TB(src_c, src_p, l-18);
            TC(dst_c, dst_l, p) = sum_l/vnorm;
        }
        PRINTF(" - column %d computed at cycle %d\n", p, proctime());
    }

    date  = proctime();
    PRINTF("*** Completing vertical filter at cycle %d\n", date);
    VERT_ENDED[cid][lid] = date;

    barrier_wait(2);

    ////////////////////////////////////////////////////////////////
    // final computation and parallel display 
    // Z <= D + C
    // Each processor use its private DMA channel to display 
    // the resulting image, line  per line (one byte per pixel).
    // Eah processor computes & displays (NL/ntasks) lines. 

    date  = proctime();
    PRINTF("\n*** Starting display at cycle %d\n", date);
    DISP_START[cid][lid] = date;

    first = lid*lines_per_task;
    last  = first + lines_per_task;

    for ( l=first ; l<last ; l++)
    {
        for ( p=0 ; p<NP ; p++)
        {
           TZ(cid,l,p) = (unsigned char)(((TD(cid,l,p) + TC(cid,l,p))>>8) & 0xFF);
        }
        fb_write(NP*(cid*lines_per_cluster+l), &TZ(cid,l,0), NP);
    }

    date  = proctime();
    PRINTF("*** Completing display at cycle %d\n", date);
    DISP_ENDED[cid][lid] = date;

    barrier_wait(3);

    /////////////////////////////////////////////////////////
    // Instrumentation (done by processor 0 in cluster 0)   

    if ( pid == 0 )
    {
        date  = proctime();
        PRINTF("\n*** Starting Instrumentation at cycle %d\n\n", date);

        int cc, pp;
        unsigned int min_load_start = 1000000000;
        unsigned int max_load_start = 0;
        unsigned int min_load_ended = 1000000000;
        unsigned int max_load_ended = 0;

        unsigned int min_hori_start = 1000000000;
        unsigned int max_hori_start = 0;
        unsigned int min_hori_ended = 1000000000;
        unsigned int max_hori_ended = 0;

        unsigned int min_vert_start = 1000000000;
        unsigned int max_vert_start = 0;
        unsigned int min_vert_ended = 1000000000;
        unsigned int max_vert_ended = 0;

        unsigned int min_disp_start = 1000000000;
        unsigned int max_disp_start = 0;
        unsigned int min_disp_ended = 1000000000;
        unsigned int max_disp_ended = 0;

        for ( cc=0 ; cc<nclusters ; cc++ )
        {
            for ( pp=0 ; pp<nprocs ; pp++ )
            {
                if ( LOAD_START[cc][pp] < min_load_start ) min_load_start = LOAD_START[cc][pp];
                if ( LOAD_START[cc][pp] > max_load_start ) max_load_start = LOAD_START[cc][pp];
                if ( LOAD_ENDED[cc][pp] < min_load_ended ) min_load_ended = LOAD_ENDED[cc][pp];
                if ( LOAD_ENDED[cc][pp] > max_load_ended ) max_load_ended = LOAD_ENDED[cc][pp];

                if ( HORI_START[cc][pp] < min_hori_start ) min_hori_start = HORI_START[cc][pp];
                if ( HORI_START[cc][pp] > max_hori_start ) max_hori_start = HORI_START[cc][pp];
                if ( HORI_ENDED[cc][pp] < min_hori_ended ) min_hori_ended = HORI_ENDED[cc][pp];
                if ( HORI_ENDED[cc][pp] > max_hori_ended ) max_hori_ended = HORI_ENDED[cc][pp];

                if ( VERT_START[cc][pp] < min_vert_start ) min_vert_start = VERT_START[cc][pp];
                if ( VERT_START[cc][pp] > max_vert_start ) max_vert_start = VERT_START[cc][pp];
                if ( VERT_ENDED[cc][pp] < min_vert_ended ) min_vert_ended = VERT_ENDED[cc][pp];
                if ( VERT_ENDED[cc][pp] > max_vert_ended ) max_vert_ended = VERT_ENDED[cc][pp];

                if ( DISP_START[cc][pp] < min_disp_start ) min_disp_start = DISP_START[cc][pp];
                if ( DISP_START[cc][pp] > max_disp_start ) max_disp_start = DISP_START[cc][pp];
                if ( DISP_ENDED[cc][pp] < min_disp_ended ) min_disp_ended = DISP_ENDED[cc][pp];
                if ( DISP_ENDED[cc][pp] > max_disp_ended ) max_disp_ended = DISP_ENDED[cc][pp];
            }
        }
        PRINTF(" - LOAD_START : min = %d / max = %d / med = %d / delta = %d\n",
        min_load_start, max_load_start, (min_load_start+max_load_start)/2, max_load_start-min_load_start);
        PRINTF(" - LOAD_END   : min = %d / max = %d / med = %d / delta = %d\n",
        min_load_ended, max_load_ended, (min_load_ended+max_load_ended)/2, max_load_ended-min_load_ended);

        PRINTF(" - HORI_START : min = %d / max = %d / med = %d / delta = %d\n",
        min_hori_start, max_hori_start, (min_hori_start+max_hori_start)/2, max_hori_start-min_hori_start);
        PRINTF(" - HORI_END   : min = %d / max = %d / med = %d / delta = %d\n",
        min_hori_ended, max_hori_ended, (min_hori_ended+max_hori_ended)/2, max_hori_ended-min_hori_ended);

        PRINTF(" - VERT_START : min = %d / max = %d / med = %d / delta = %d\n",
        min_vert_start, max_vert_start, (min_vert_start+max_vert_start)/2, max_vert_start-min_vert_start);
        PRINTF(" - VERT_END   : min = %d / max = %d / med = %d / delta = %d\n",
        min_vert_ended, max_vert_ended, (min_vert_ended+max_vert_ended)/2, max_vert_ended-min_vert_ended);

        PRINTF(" - DISP_START : min = %d / max = %d / med = %d / delta = %d\n",
        min_disp_start, max_disp_start, (min_disp_start+max_disp_start)/2, max_disp_start-min_disp_start);
        PRINTF(" - DISP_END   : min = %d / max = %d / med = %d / delta = %d\n",
        min_disp_ended, max_disp_ended, (min_disp_ended+max_disp_ended)/2, max_disp_ended-min_disp_ended);

        PRINTF(" - BARRIER LOAD/HORI = %d\n", min_hori_start - max_load_ended);
        PRINTF(" - BARRIER HORI/VERT = %d\n", min_vert_start - max_hori_ended);
        PRINTF(" - BARRIER VERT/DISP = %d\n", min_disp_start - max_vert_ended);

        PRINTF(" - LOAD              = %d\n", max_load_ended);
        PRINTF(" - FILTER            = %d\n", max_vert_ended - max_load_ended);
        PRINTF(" - DISPLAY           = %d\n", max_disp_ended - max_vert_ended);

        PRINTF("\nBEGIN LOAD_START\n");
        for ( cc=0 ; cc<nclusters ; cc++ )
        {
            for ( pp=0 ; pp<nprocs ; pp++ )
            {
                PRINTF("cluster= %d proc= %d date= %d\n", cc, pp, LOAD_START[cc][pp]);  
            }
        }
        PRINTF("END\n");
        PRINTF("\nBEGIN LOAD_ENDED\n");
        for ( cc=0 ; cc<nclusters ; cc++ )
        {
            for ( pp=0 ; pp<nprocs ; pp++ )
            {
                PRINTF("cluster= %d proc= %d date= %d\n", cc, pp, LOAD_ENDED[cc][pp]);  
            }
        }
        PRINTF("END\n");
        PRINTF("\nBEGIN HORI_START\n");
        for ( cc=0 ; cc<nclusters ; cc++ )
        {
            for ( pp=0 ; pp<nprocs ; pp++ )
            {
                PRINTF("cluster= %d proc= %d date= %d\n", cc, pp, HORI_START[cc][pp]);  
            }
        }
        PRINTF("END\n");
        PRINTF("\nBEGIN HORI_ENDED\n");
        for ( cc=0 ; cc<nclusters ; cc++ )
        {
            for ( pp=0 ; pp<nprocs ; pp++ )
            {
                PRINTF("cluster= %d proc= %d date= %d\n", cc, pp, HORI_ENDED[cc][pp]);  
            }
        }
        PRINTF("END\n");
        PRINTF("\nBEGIN VERT_START\n");
        for ( cc=0 ; cc<nclusters ; cc++ )
        {
            for ( pp=0 ; pp<nprocs ; pp++ )
            {
                PRINTF("cluster= %d proc= %d date= %d\n", cc, pp, VERT_START[cc][pp]);  
            }
        }
        PRINTF("END\n");
        PRINTF("\nBEGIN VERT_ENDED\n");
        for ( cc=0 ; cc<nclusters ; cc++ )
        {
            for ( pp=0 ; pp<nprocs ; pp++ )
            {
                PRINTF("cluster= %d proc= %d date= %d\n", cc, pp, VERT_ENDED[cc][pp]);  
            }
        }
        PRINTF("END\n");
        PRINTF("\nBEGIN DISP_START\n");
        for ( cc=0 ; cc<nclusters ; cc++ )
        {
            for ( pp=0 ; pp<nprocs ; pp++ )
            {
                PRINTF("cluster= %d proc= %d date= %d\n", cc, pp, DISP_START[cc][pp]);  
            }
        }
        PRINTF("END\n");
        PRINTF("\nBEGIN DISP_ENDED\n");
        for ( cc=0 ; cc<nclusters ; cc++ )
        {
            for ( pp=0 ; pp<nprocs ; pp++ )
            {
                PRINTF("cluster= %d proc= %d date= %d\n", cc, pp, DISP_ENDED[cc][pp]);  
            }
        }
        PRINTF("END\n");
    }

    while(1);

} // end main()