Dnia Fri, 03 May 2013 04:02:49 -0700 po głębokim namyśle firr kenobi
rzekł:

> drugie dla mnie najciekawze zdanie
>
> "Systems programming is definitely in no way a goal. If that's the
> target, C is still the language of choice, of course. However, C is
> still the golden standard for general performance -- although, of
> course, Fortran beats it in some cases."
>
> jest to tajemnica ktora sie przejawia od jakiegos czasu w moim zyciu - w
> jakich to cases fortran beats c ??? (mozliwe ze tak jest ale to raczej z
> jakichs 'pzypadkowych' powodow typu pewnie cos ze zmiennymi lokalnymi
> albo conwencja wywolan funkcji w c albo cos takiego - pozatym raczej
> chyba nie ma powodu wiec nie powinno tak byc ale chetnie bym sie
> dowiedzial)

Fortran przez lata królował w obliczeniach numerycznych, teraz już
bardziej używane jest C++.

Na pytanie 'dlaczego' nie ma trywialnej odpowiedzi: Fortran numerycznie
jest szybszy od C dość często, Python jest wolniejszy tak ze 20x minimum
praktycznie zawsze - tak po prostu jest, zaakceptuj to jak pogodę.

> z innych ciekawostek: CO TO ROBI ?
>
> function randmatstat(t)
> n = 5 v = zeros(t)
> w = zeros(t)
> for i = 1:t
> a = randn(n,n)
> b = randn(n,n)
> c = randn(n,n)
> d = randn(n,n)
> P = [a b c d]
> Q = [a b; c d]
> v[i] = trace((P.'*P)^4)
> w[i] = trace((Q.'*Q)^4)
> end std(v)/mean(v), std(w)/mean(w)
> end

Wklejam wersję c. Od "co" ważniejsze dla mnie jest "jak".

Edek

PS. Kod w c:

struct double_pair randmatstat(int t) {
int n = 5;
struct double_pair r;
double *v = (double*)calloc(t,sizeof(double));
double *w = (double*)calloc(t,sizeof(double));
double *a = (double*)malloc((n)*(n)*sizeof(double));
double *b = (double*)malloc((n)*(n)*sizeof(double));
double *c = (double*)malloc((n)*(n)*sizeof(double));
double *d = (double*)malloc((n)*(n)*sizeof(double));
double *P = (double*)malloc((n)*(4*n)*sizeof(double));
double *Q = (double*)malloc((2*n)*(2*n)*sizeof(double));
double *PtP1 = (double*)malloc((4*n)*(4*n)*sizeof(double));
double *PtP2 = (double*)malloc((4*n)*(4*n)*sizeof(double));
double *QtQ1 = (double*)malloc((2*n)*(2*n)*sizeof(double));
double *QtQ2 = (double*)malloc((2*n)*(2*n)*sizeof(double));
for (int i=0; i < t; i++) {
randmtzig_fill_randn(a, n*n);
randmtzig_fill_randn(b, n*n);
randmtzig_fill_randn(c, n*n);
randmtzig_fill_randn(d, n*n);
memcpy(P+0*n*n, a, n*n*sizeof(double));
memcpy(P+1*n*n, b, n*n*sizeof(double));
memcpy(P+2*n*n, c, n*n*sizeof(double));
memcpy(P+3*n*n, d, n*n*sizeof(double));
for (int j=0; j < n; j++) {
for (int k=0; k < n; k++) {
Q[2*n*j+k] = a[k];
Q[2*n*j+n+k] = b[k];
Q[2*n*(n+j)+k] = c[k];
Q[2*n*(n+j)+n+k] = d[k];
}
}
cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans,
n, n, 4*n, 1.0, P, 4*n, P, 4*n, 0.0, PtP1, 4*n);
cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
4*n, 4*n, 4*n, 1.0, PtP1, 4*n, PtP1, 4*n, 0.0, PtP2,
4*n);
cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
4*n, 4*n, 4*n, 1.0, PtP2, 4*n, PtP2, 4*n, 0.0, PtP1,
4*n);
for (int j=0; j < n; j++) {
v[i] += PtP1[(n+1)*j];
}
cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans,
2*n, 2*n, 2*n, 1.0, Q, 2*n, Q, 2*n, 0.0, QtQ1, 2*n);
cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
2*n, 2*n, 2*n, 1.0, QtQ1, 2*n, QtQ1, 2*n, 0.0, QtQ2,
2*n);
cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
2*n, 2*n, 2*n, 1.0, QtQ2, 2*n, QtQ2, 2*n, 0.0, QtQ1,
2*n);
for (int j=0; j < 2*n; j++) {
w[i] += QtQ1[(2*n+1)*j];
}
}
free(PtP1);
free(PtP2);
free(QtQ1);
free(QtQ2);
free(P);
free(Q);
free(a);
free(b);
free(c);
free(d);
double v1=0.0, v2=0.0, w1=0.0, w2=0.0;
for (int i=0; i < t; i++) {
v1 += v[i]; v2 += v[i]*v[i];
w1 += w[i]; w2 += w[i]*w[i];
}
free(v);
free(w);
r.s1 = sqrt((t*(t*v2-v1*v1))/((t-1)*v1*v1));
r.s2 = sqrt((t*(t*w2-w1*w1))/((t-1)*w1*w1));
return r;
}