On 2 Lut, 17:47, Fil <f...@p...onet.pl> wrote:
> Witam!
>
> Są dwie procedury:
> CPP:
> void MulTab(int N, int Q, int M, double** A, double** B, double** C)
> {
>   for (int i = 0; i < N; ++i)
>     for (int j = 0; j < M; ++j) {
>       C[i][j] = 0;
>       for (int k = 0; k < Q; ++k) C[i][j] += A[i][k] * B[k][j];
>     }
>   }
>

Troche na boku dyskusji o wyzszosci siąt bozego narodzenia.
Zle to robisz;) Przynajmniej zamien petle meijscami, lepsza
kolejnosc to ikj. Przy takich zadaniach RAM jest wolny,
wiec trzeba pomagać danym trzymac cie w cache.

Odrwocenie petli pomaga. Jeszcze lepiej dziala mnozenie blokami.


test 1: n=384 najlepszy z 10 prob zwykly piorytet
test 2 n=960 najlepszy z 3 prob, real time.

void dgemm_ijk(double *A, double *B, double *C, int N)
{
int i,j,k;
for (i = 0; i < N; i++)
for (j = 0; j < N; j++)
for (k = 0; k < N; k++)
C[i*N+j] += A[i*N+k]*B[k*N+j];
}


1) 3.738s
2) ad mortem defacatum (nie bral udzialu w drugim tescie)

Odwrocmy kolejnosc. Obie macierze wejsciowe czytane sa
w prawidlowej kolejnosci.

void dgemm_ikj(double *A, double *B, double *C, int N)
{
int i,j,k;
for (i = 0; i < N; i++)
for (k = 0; k < N; k++)
for (j = 0; j < N; j++)
C[i*N+j] += A[i*N+k]*B[k*N+j];
}

1)1.011
2) 10.299 (a jest 3 razy szybszy niz wersja trywialna!)

Blkokowo, pewna podmacierz miesci sie w cache, to ja tam
zostawmy i wykorzystajmy ile sie da.

template <class T,int SM> void dgemm_bikj(T *A, T *B, T *C, int N)
{
int i,j,k,ii,kk,jj;
for (i = 0; i < N; i+=SM)
for (k = 0; k < N; k+=SM)
for (j = 0; j < N; j+=SM)
for (ii = i; ii < i+SM; ii++)
for (kk = k; kk < k+SM; kk++)
for (jj = j; jj < j+SM; jj++)
C[ii*N+jj] += A[ii*N+kk]*B[kk*N+jj];
}

SM=30

1) 0.174
2) 2.639

Kompilator nie wpadl na wszytkie mozliwe zabawy ze wskaznikami,
dodajemy brzydkie chakierstwa.

template <class T,int SM> void dgemm_bikj2(T *A, T *B, T *C, int N)
{
int i,j,k,ii,kk,jj;
double *AA,*BB,*CC;
for (i = 0; i < N; i+=SM)
for (k = 0; k < N; k+=SM)
for (j = 0; j < N; j+=SM)
for (AA=A+N*i,CC=C+N*i,ii = i; ii < i+SM; ii++,AA+=N,CC+=N)
for (BB=B+N*k, kk = k ; kk < k+SM; kk++,BB+=N)
for (jj = j; jj < j+SM; jj++)
CC[jj] = AA[kk]*BB[jj];

}

SM=32
1) 0.122
2) 1.961 1.841 (SM=64)
Przy 30s w wersji trywialnej.

Matlab (intel math kernel najpewniej), bez sse2
>> A=rand(N);B=rand(N);
>> tic, A*B; toc

1) 0.111
2) 1.560

Podsumowując, kod trywialny jest jest 20 razy gorszy niz gotowce
(atlas powinien byc niewiele gorszy niz matalb). Mimo zawziecia
i nawet zastosowania hakerstwa (te wskazniki, cos ich VC++ nie
zoptymalizowal w przypadku dgemm_bikj jak w dgemm_bikj2)
nie udalo sie dobić do wynikow pakietow (a moje wyniki sa naciagniete,
zakaldam, ze N jest wieloktornoscia szerokosci bloku!)

Kod powstal do podobnej dyskusji tu albo na *.c
Warto przejrzec zalinkowany wtedy material (to raczej podstawy)
http://wazniak.mimuw.edu.pl/index.php?title=MN06#Jak
_napisa.C4.87_kod_.C5.BAle_wykorzystuj.C4.85cy_pami.
C4.99.C4.87_podr.C4.99czn.C4.85.3F


A, wracajac do pierwotnehgo problemu, warto zerknac, czy C i pascal w
tej samej
kolejnosci odkladaja macierze (dwuwymairowe tablice) w pamieci.
Czesciej w programach numerycznych uzywa sie jednej tablicy
(tak jest w blas/lapack, tak jest w powyzszym kodzie) a nie tablic
2D.



pozdrawiam
bartekltg