Re: [ublas] Matrix multiplication performance

Subject: Re: [ublas] Matrix multiplication performance
From: Joaquim Duran Comas (jdurancomas_at_[hidden])
Date: 2016-02-02 18:56:36

Michael,

I've a comment about the function malloc_: it is possible to reduce the
memory allocated, exactly 1 byte per call. The change to implement is:

void *
malloc_(std::size_t alignment, std::size_t size)
{
    alignment = std::max(alignment, alignof(void *));
    size     += alignment - 1;  // this is the change: introduce -1

    void *ptr  = std::malloc(size);
    void *ptr2 = (void *)(((uintptr_t)ptr + alignment) & ~(alignment-1));
    void **vp  = (void**) ptr2 - 1;
    *vp        = ptr;
    return ptr2;
}

The reason for this change is that to make sure that an aligned
address is allocated it is not needed to add alignment to size,
because the allocation of 'size' bytes is the first address.

As a example, let's assume that we want allocate 1 byte of data
aligned at 8 (ignore the space to allocate ptr address for simplicity)
and the first memory address is 0. Based on the original code, 9 bytes
should be allocated (1 + 8), addresses from 0 until 8, so it includes
two 8 bytes aligned addresses (addresses 0 and 8). However, based on
the change proposed, 8 bytes should be allocated (1 + 8 - 1),
addresses from 0 till 7, and only includes one 8 bytes aligned address
(address 0).

Joaquim Duran

2016-02-02 23:26 GMT+01:00 Michael Lehn <michael.lehn_at_[hidden]>:

> I also started to write down some notes on the rational behind the
> micro-kernel:
>
> http://www.mathematik.uni-ulm.de/~lehn/test_ublas/session6/page01.html
>
>
> On 02 Feb 2016, at 21:43, Michael Lehn <michael.lehn_at_[hidden]> wrote:
>
> (So the pages are re-done …)
>
>
>
> The bug is because NR is BlockSize<T>::NR divided by the vector length.
> So this
> should work
>
>
> //---------------------------------------------------------------------------------
> template <typename Index, typename T>
> typename std::enable_if<std::is_floating_point<T>::value,
> void>::type
> ugemm(Index kc, T alpha,
>                const T *A, const T *B,
>      T beta,
>      T *C, Index incRowC, Index incColC)
> {
>    A = (const T*) __builtin_assume_aligned (A, 128);
>    B = (const T*) __builtin_assume_aligned (B, 128);
>    static const unsigned div = 32/sizeof(T);
>    static const Index MR = BlockSize<T>::MR;
>    static const Index NR = BlockSize<T>::NR/div;
>
>    typedef T vx __attribute__((vector_size (32)));
>
>    vx P[MR*NR] __attribute__ ((aligned (128))) = {};
>    const vx *B_ = (vx *)B;
>    for (Index l=0; l<kc; ++l) {
>        for (Index i=0; i<MR; ++i) {
>            for (Index j=0; j<(NR); ++j) {
>                P[i * NR + j] += A[i]*B_[j];
>            }
>        }
>        A += MR;
>        B_ += NR;
>    }
>
>    T *P_ = (T *)P;
>    for (Index j=0; j<NR*div; ++j) {
>        for (Index i=0; i<MR; ++i) {
>            C[i*incRowC+j*incColC] *= beta;
>            C[i*incRowC+j*incColC] += alpha*P_[i * NR * div + j];
>        }
>    }
> }
>
> //---------------------------------------------------------------------------------
>
> However, as Karl pointed out it also should treat the case beta==0
> as a special case.
>
> But I don’t get any speedup.  For the performance relevant is
> mainly the loop
>
>    for (Index l=0; l<kc; ++l) {
>        for (Index i=0; i<MR; ++i) {
>            for (Index j=0; j<(NR); ++j) {
>                P[i * NR + j] += A[i]*B_[j];
>            }
>        }
>        A += MR;
>        B_ += NR;
>    }
>
> From my experience, only if this is optimized to the extrem a little bit
> extra performance can be achieved by treating some special cases for
> the update
>
>    T *P_ = (T *)P;
>    for (Index j=0; j<NR*div; ++j) {
>        for (Index i=0; i<MR; ++i) {
>            C[i*incRowC+j*incColC] *= beta;
>            C[i*incRowC+j*incColC] += alpha*P_[i * NR * div + j];
>        }
>    }
>
> The main reason my code looks weird is that I also want to use it
> for teaching.  The way it is coded you easily can map it forth and
> back with the assembly code.
>
>
> new kernel:
>
> template <typename Index, typename T>
> typename std::enable_if<std::is_floating_point<T>::value,
> void>::type
> ugemm(Index kc, T alpha,
>                const T *A, const T *B,
>      T beta,
>      T *C, Index incRowC, Index incColC)
> {
>    A = (const T*) __builtin_assume_aligned (A, 128);
>    B = (const T*) __builtin_assume_aligned (B, 128);
>    static const unsigned div = 32/sizeof(T);
>    static const Index MR = BlockSize<T>::MR;
>    static const Index NR = BlockSize<T>::NR/div;
>
>    typedef T vx __attribute__((vector_size (32)));
>
>    vx P[MR*NR] __attribute__ ((aligned (128))) = {};
>    const vx *B_ = (vx *)B;
>    for (Index l=0; l<kc; ++l) {
>        for (Index i=0; i<MR; ++i) {
>            for (Index j=0; j<(NR); ++j) {
>                P[i * NR + j] += A[i]*B_[j];
>            }
>        }
>        A += MR;
>        B_ += NR;
>    }
>
>    T *P_ = (T *)P;
>    for (Index j=0; j<NR; ++j) {
>        for (Index i=0; i<MR; ++i) {
>            C[i*incRowC+j*incColC] *= beta;
>            C[i*incRowC+j*incColC] += alpha*P_[i * NR + j];
>        }
>    }
> }
>
>
> _______________________________________________
> ublas mailing list
> ublas_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/ublas
> Sent to: michael.lehn_at_[hidden]
>
>
>
> _______________________________________________
> ublas mailing list
> ublas_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/ublas
> Sent to: jdurancomas_at_[hidden]
>