Subject: Re: [ublas] Matrix multiplication performance
From: Michael Lehn (michael.lehn_at_[hidden])
Date: 2016-02-10 09:50:00
Hi Palik,
sorry for the late response. It’s the last week of the semester so it is pretty busy … It is pretty
cool that prefetching improves the results so significantly. I will add this to the examples on my
page.
Today I added an example on how-to implement a fast LU factorization with ublas. The
implementation still needs polishing but the algorithms are pretty good:
http://www.mathematik.uni-ulm.de/~lehn/test_ublas/session7/page01.html
Cheers,
Michael
On 03 Feb 2016, at 20:53, palik imre <imre_palik_at_[hidden]> wrote:
> Hi Michael,
>
> first of all, the link on your page to test_ublas.tgz is still pointing to session4.tgz
>
> I also realised, that we start to become memory bound.
> I am getting approximately 10% improvement with the following setup on my old AMD box:
>
>
> 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_ = (const vx *)B;
> for (Index l=0; l<kc; ++l) {
> if (!(kc&3))
> __builtin_prefetch(A + 16, 0);
> 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];
> }
> }
> }
>
> template <typename Index, typename T, typename Beta, typename TC>
> void
> mgemm(Index mc, Index nc, Index kc,
> T alpha,
> const T *A, const T *B,
> Beta beta,
> TC *C, Index incRowC, Index incColC)
> {
> const Index MR = BlockSize<T>::MR;
> const Index NR = BlockSize<T>::NR;
> const Index mp = (mc+MR-1) / MR;
> const Index np = (nc+NR-1) / NR;
> const Index mr_ = mc % MR;
> const Index nr_ = nc % NR;
>
> #if defined(_OPENMP)
> #pragma omp parallel for
> #endif
> for (Index j=0; j<np; ++j) {
> const Index nr = (j!=np-1 || nr_==0) ? NR : nr_;
> T C_[BlockSize<T>::MR*BlockSize<T>::NR];
> __builtin_prefetch(B + j * kc * NR, 0);
> __builtin_prefetch(B + j * kc * NR + 16, 0);
> for (Index i=0; i<mp; ++i) {
> const Index mr = (i!=mp-1 || mr_==0) ? MR : mr_;
> __builtin_prefetch(A + i * kc * MR, 0);
> if (mr==MR && nr==NR) {
> ugemm(kc, alpha,
> &A[i*kc*MR], &B[j*kc*NR],
> beta,
> &C[i*MR*incRowC+j*NR*incColC],
> incRowC, incColC);
> } else {
> std::fill_n(C_, MR*NR, T(0));
> ugemm(kc, alpha,
> &A[i*kc*MR], &B[j*kc*NR],
> T(0),
> C_, Index(1), MR);
> gescal(mr, nr, beta,
> &C[i*MR*incRowC+j*NR*incColC],
> incRowC, incColC);
> geaxpy(mr, nr, T(1), C_, Index(1), MR,
> &C[i*MR*incRowC+j*NR*incColC],
> incRowC, incColC);
> }
> }
> }
> }
>