Commit 251e344f authored by Martin Perdacher's avatar Martin Perdacher
Browse files

performance improvement and bugfix

parent 740ae767

cmake/.gitignore

+2 −0
Original line number Diff line number Diff line
@@ -2,3 +2,5 @@ build/ 
debug/

nohup.out

*.csv

*.sh

bmatrix

cmake/CMakeLists.txt

+1 −0
Original line number Diff line number Diff line
@@ -14,6 +14,7 @@ elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel") 
  set(CMAKE_CXX_FLAGS  "${CMAKE_CXX_FLAGS} -xmic-avx512 -qopenmp -qopt-report=2 -DNDEBUG -O3 -liomp5 -lpthread -lmkl_core -lmkl_intel_lp64 -lmkl_intel_thread")

  # xeon phi (knl) specific:

  set(CMAKE_CXX_FLAGS  "${CMAKE_CXX_FLAGS} -lmemkind")

  add_definitions(-DDEBUG)

endif()



add_definitions(-DNUM_THREADS=64)

cmake/blasJoin/blasJoin.cpp

+28 −33
Original line number Diff line number Diff line 


#include "blasJoin.h"



double blasJoin(const double *x, const size_t N, const size_t D, const double EPS, const unsigned int THREADS, const size_t BLOCKSIZE){

double blasJoin(const double *x, const size_t N, const size_t D, const double EPS, const unsigned int THREADS, const size_t BLOCKSIZE, size_t *joinCounts){



    assert( BLOCKSIZE < N && BLOCKSIZE < 21000 );

    assert( BLOCKSIZE < N && BLOCKSIZE < 21000 && BLOCKSIZE > 1);



    double elapsed=0.0;

    CUtilTimer timer;

    double *iresult = NULL;

    // bool result[(N * N )/ 2 - N];

    int joinCounts = 0;



    *joinCounts = 0;



    if ( NUM_THREADS != 0 ){

         omp_set_num_threads(THREADS);
@@ -17,10 +17,10 @@ double blasJoin(const double *x, const size_t N, const size_t D, const double EP 
    }



    double *p = (double*) ddr_alloc(sizeof(double) * N );

    const double EPS_SQUARED = (EPS * EPS) / 4;

    const double EPS_SQUARED = (EPS * EPS) / 4.0 ;



    // pre-compute the values P(i) and P(j)

    #pragma simd

    // #pragma simd

    for ( int i = 0 ; i < N ; i++ ){

        p[i] = EPS_SQUARED - cblas_ddot(D, &x[i*D], 1, &x[i*D], 1) / 2.0;

    }
@@ -33,35 +33,36 @@ double blasJoin(const double *x, const size_t N, const size_t D, const double EP 
    size_t blockRow = 0, blockCol=0;



    timer.start();

    // lower triangle excluding the main diagonal



    #pragma ivdep

    #pragma vector always

    // #pragma ivdep

    // #pragma vector always

    for ( int i = 1 ; i < NRC ; i++ ){

        blockRow = ( i + 1 == NRC ) ? N - (BLOCKSIZE * i ) : BLOCKSIZE;

        blockRow = ( i + 1 == NRC ) ? N - (BLOCKSIZE * (NRC - 1) ) : BLOCKSIZE;

        for ( int j = 0 ; j < i ; j++ ){

            blockCol = ( j + 1 == NRC ) ? N - (BLOCKSIZE * (j + 1) ) : BLOCKSIZE;

            blockCol = ( j + 1 == NRC ) ? N - (BLOCKSIZE * (NRC - 1) ) : BLOCKSIZE;



            // perform regular matrix multiplication

            // C := alpha*A*B' + beta*C



            cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, blockRow, blockCol, D, 1.0, x, D, x, D, 0.0, iresult, blockCol);

            int accum=0;

            cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, blockRow, blockCol, D, 1.0, &x[i*BLOCKSIZE*D], D, &x[j*BLOCKSIZE*D], D, 0.0, iresult, blockCol);



            // count join partners on whole matrix

            #pragma omp parallel for reduction(+:accum) collapse(2)

            /*

               -count join partners on whole matrix

               -collapse(2) won't work here, since it is triangular

            */

            size_t accum=0;

            #pragma omp parallel for reduction(+:accum)

            for ( int k = 0 ; k < blockRow ; k++ ){

                for ( int l =0 ; l < blockCol ; l++ ){

                for ( int l =0 ; l < k ; l++ ){

                    /* with the 63-rshift we find out whether we are +/- and sum it up to efficiently determine count of join-partners

                     * prevent RAW-dependency with one-liner instead of the following:

                        accum +=  ( (long long) iresult[k * NRC + l] ) >> 63;

                        iresult[k*blockCol + l] += ( p[i*BLOCKSIZE + l] + p[i*BLOCKSIZE+k] );

                     * prevent RAW-dependency with one-liner

                     * count only join-partners (positive ones)

                    */

                    accum += ( (long long) (iresult[k * blockCol + l] + p[i*BLOCKSIZE+k] + p[j*BLOCKSIZE+l]) ) >> 63;

                    accum +=  - ( (long long) floor( - ( iresult[k * BLOCKSIZE + l] + p[i*BLOCKSIZE+k] + p[j*BLOCKSIZE+l] ) ) >> 63);

                }

            }



            joinCounts += accum;

            *joinCounts += accum;

        }

    }
@@ -71,28 +72,22 @@ double blasJoin(const double *x, const size_t N, const size_t D, const double EP 
        // perform symmetric matrix multiplication on lower triangle only

        // C := alpha*A*A' + beta*C



        // printf("%d %d %d\n", i, blockRow, NRC);fflush(stdout);

        cblas_dsyrk(CblasRowMajor, CblasLower, CblasNoTrans, blockRow, D, 1.0, x + BLOCKSIZE * D * i, D, 0.0, iresult, blockRow);



        int accum=0;

        // count join partners on lower triangle

        #pragma omp parallel for reduction(+:accum) collapse(2)

        size_t accum=0;

        #pragma omp parallel for reduction(+:accum)

        for ( int k = 0 ; k < blockRow ; k++ ){

            for ( int l = 0 ; l < k ; l++ ){

                accum += ( (long long) (iresult[k * blockCol + l] + p[i*BLOCKSIZE+k] + p[i*BLOCKSIZE+l]) ) >> 63;

                accum += - ( (long long) floor( - ( iresult[k * BLOCKSIZE + l] + p[i*BLOCKSIZE+k] + p[i*BLOCKSIZE+l] ) ) >> 63);

            }

        }



        joinCounts += accum;

        *joinCounts += accum;

    }



    timer.stop();



    joinCounts = ((N * N) / 2) - N + joinCounts;

#ifdef DEBUG

    printf("joinCounts: %d\n", joinCounts);

#endif

    ddr_free(iresult);

    ddr_free(p);



    return timer.get_time();

}

cmake/blasJoin/blasJoin.h

+1 −1
Original line number Diff line number Diff line
@@ -13,6 +13,6 @@ 
#include "../util/allocation.h"

#include "../util/timer.h"



double blasJoin(const double *x, const size_t N, const size_t D, const double EPS, const unsigned int THREADS, const size_t BLOCKSIZE);

double blasJoin(const double *x, const size_t N, const size_t D, const double EPS, const unsigned int THREADS, const size_t BLOCKSIZE, size_t *joinCounts);



#endif

cmake/main.cpp

+17 −7
Original line number Diff line number Diff line
@@ -8,20 +8,20 @@ 
#include "blasJoin/blasJoin.h"



int main(int argc, char** argv) {

    double *x = NULL;



    CUtilTimer timer;

    char filename[] = "";

    double elapsed=0.0;

    bool isBinary = false;



    double *x = NULL;

    size_t threads = 64;

    size_t N=200, D=20;

    size_t blocksize=4000;

    size_t blocksize=4000, joinCounts;

    double EPS=0.2;



    parsing_args(argc,argv, &N, &EPS, &D, &threads, &blocksize, filename, isBinary);

    parsing_args(argc, argv, &N, &EPS, &D, &threads, &blocksize, filename, &isBinary);



    N = N * 1000;

    // N = N * 1000;



    if ( threads != 0 ){

        omp_set_num_threads(threads);
@@ -29,18 +29,28 @@ int main(int argc, char** argv) { 


    x = (double*) ddr_alloc(sizeof (double)*N * D);



    // size_t threads = 1;

    // size_t N=4, D=3;

    // size_t blocksize=2;

    // double EPS=0.3;

    // double x[12] = { 1.0000, 2.0000, 3.0000,

    //                  1.1000, 2.2000, 3.0000,

    //                  5.0000, 6.0000, 7.0000,

    //                  5.1000, 6.2000, 7.1000

    //              };



    if ( strcmp(filename,"" ) == 0) {

        random_init(x,N,D);

    }else{

        read_file(x, N, D, filename, isBinary);

    }



    elapsed = blasJoin( x, N, D, EPS, threads, blocksize);

    elapsed = blasJoin( x, N, D, EPS, threads, blocksize, &joinCounts);



#pragma omp parallel

{

    if (  omp_get_thread_num() == 0 ){

        printf("%ld; %ld; %f; %ld; %d; %f \n", N, D, EPS, blocksize, omp_get_num_threads(),elapsed);

        printf("%ld; %ld; %f; %ld; %d; %f ; %lu\n", N, D, EPS, blocksize, omp_get_num_threads(),elapsed, joinCounts);

    }

}