← Back to branch summary

~ubuntu-branches/ubuntu/wily/julia/wily

« back to all changes in this revision

Viewing changes to test/perf/perf.c

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-01-16 12:29:42 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130116122942-x86e42akjq31repw
Tags: 0.0.0+20130107.gitd9656f41-1
* New upstream snashot
* No longer try to rebuild helpdb.jl.
   + debian/rules: remove helpdb.jl from build-arch rule
   + debian/control: move back python-sphinx to Build-Depends-Indep
* debian/copyright: reflect upstream changes
* Add Build-Conflicts on libatlas3-base (makes linalg tests fail)
* debian/rules: replace obsolete USE_DEBIAN makeflag by a list of
  USE_SYSTEM_* flags
* debian/rules: on non-x86 systems, use libm instead of openlibm
* dpkg-buildflags.patch: remove patch, applied upstream
* Refreshed other patches

Show diffs side-by-side

added added

removed removed

Lines of Context:
test/perf/perf.c
 
1
#include <complex.h>
 
2
 
 
3
#define DSFMT_MEXP 19937
 
4
#include "perf.h"
 
5
#include "../../deps/random/randmtzig.c"
 
6
 
 
7
double *myrand(int n) {
 
8
    double *d = (double *)malloc(n*sizeof(double));
 
9
    dsfmt_gv_fill_array_close_open(d, n);
 
10
    return d;
 
11
}
 
12
 
 
13
#define NITER 5
 
14
 
 
15
double clock_now()
 
16
{
 
17
    struct timeval now;
 
18
 
 
19
    gettimeofday(&now, NULL);
 
20
    return (double)now.tv_sec + (double)now.tv_usec/1.0e6;
 
21
}
 
22
 
 
23
int fib(int n) {
 
24
    return n < 2 ? n : fib(n-1) + fib(n-2);
 
25
}
 
26
 
 
27
long parse_int(const char *s, long base) {
 
28
    long n = 0;
 
29
 
 
30
    for (unsigned int i=0; i<strlen(s); ++i) {
 
31
        char c = s[i];
 
32
        long d = 0;
 
33
        if (c >= '0' && c <= '9') d = c-'0';
 
34
        else if (c >= 'A' && c <= 'Z') d = c-'A' + (int) 10;
 
35
        else if (c >= 'a' && c <= 'z') d = c-'a' + (int) 10;
 
36
        else exit(-1);
 
37
 
 
38
        if (base <= d) exit(-1);
 
39
        n = n*base + d;
 
40
    }
 
41
    return n;
 
42
}
 
43
 
 
44
double *ones(int m, int n) {
 
45
    double *a = (double *) malloc(m*n*sizeof(double));
 
46
    for (int k=0; k<m*n; ++k) {
 
47
        a[k] = 1.0;
 
48
    }
 
49
    return a;
 
50
}
 
51
 
 
52
double *matmul_aat(int n, double *b) {
 
53
    double *c = (double *) malloc(n*n*sizeof(double));
 
54
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans, n, n, n, 1.0, b, n, b, n, 0.0, c, n);
 
55
    return c;
 
56
}
 
57
 
 
58
int mandel(double complex z) {
 
59
    int n = 0;
 
60
    double complex c = z;
 
61
    for (n=0; n<=79; ++n) {
 
62
        if (cabs(z) > 2.0) {
 
63
            n -= 1;
 
64
            break;
 
65
        }
 
66
        z = cpow(z,2)+c;
 
67
    }
 
68
    return n+1;
 
69
}
 
70
 
 
71
int mandelperf() {
 
72
    int mandel_sum = 0;
 
73
    for (double re=-2.0; re<=0.5; re+=0.1) {
 
74
        for (double im=-1.0; im<=1.0; im+=0.1) {
 
75
            int m = mandel(re+im*I);
 
76
            mandel_sum += m;
 
77
        }
 
78
    }
 
79
    return mandel_sum;
 
80
}
 
81
 
 
82
void quicksort(double *a, int lo, int hi) {
 
83
    int i = lo;
 
84
    int j = hi;
 
85
    while (i < hi) {
 
86
        double pivot = a[(lo+hi)/2];
 
87
        // Partition
 
88
        while (i <= j) {
 
89
            while (a[i] < pivot) {
 
90
                i = i + 1;
 
91
            }
 
92
            while (a[j] > pivot) {
 
93
                j = j - 1;
 
94
            }
 
95
            if (i <= j) {
 
96
                double t = a[i];
 
97
                a[i] = a[j];
 
98
                a[j] = t;
 
99
                i = i + 1;
 
100
                j = j - 1;
 
101
            }
 
102
        }
 
103
 
 
104
        // Recursion for quicksort
 
105
        if (lo < j) {
 
106
            quicksort(a, lo, j);
 
107
        }
 
108
        lo = i;
 
109
        j = hi;
 
110
    }
 
111
}
 
112
 
 
113
double pisum() {
 
114
    double sum = 0.0;
 
115
    for (int j=0; j<500; ++j) {
 
116
        sum = 0.0;
 
117
        for (int k=1; k<=10000; ++k) {
 
118
            sum += 1.0/(k*k);
 
119
        }
 
120
    }
 
121
    return sum;
 
122
}
 
123
 
 
124
struct double_pair { double s1, s2; };
 
125
 
 
126
struct double_pair randmatstat(int t) {
 
127
    int n = 5;
 
128
    struct double_pair r;
 
129
    double *v = (double*)calloc(t,sizeof(double));
 
130
    double *w = (double*)calloc(t,sizeof(double));
 
131
    double *a = (double*)malloc((n)*(n)*sizeof(double));
 
132
    double *b = (double*)malloc((n)*(n)*sizeof(double));
 
133
    double *c = (double*)malloc((n)*(n)*sizeof(double));
 
134
    double *d = (double*)malloc((n)*(n)*sizeof(double));
 
135
    double *P = (double*)malloc((n)*(4*n)*sizeof(double));
 
136
    double *Q = (double*)malloc((2*n)*(2*n)*sizeof(double));
 
137
    double *PtP1 = (double*)malloc((4*n)*(4*n)*sizeof(double));
 
138
    double *PtP2 = (double*)malloc((4*n)*(4*n)*sizeof(double));
 
139
    double *QtQ1 = (double*)malloc((2*n)*(2*n)*sizeof(double));
 
140
    double *QtQ2 = (double*)malloc((2*n)*(2*n)*sizeof(double));
 
141
    for (int i=0; i < t; i++) {
 
142
        randmtzig_fill_randn(a, n*n);
 
143
        randmtzig_fill_randn(b, n*n);
 
144
        randmtzig_fill_randn(c, n*n);
 
145
        randmtzig_fill_randn(d, n*n);
 
146
        memcpy(P+0*n*n, a, n*n*sizeof(double));
 
147
        memcpy(P+1*n*n, b, n*n*sizeof(double));
 
148
        memcpy(P+2*n*n, c, n*n*sizeof(double));
 
149
        memcpy(P+3*n*n, d, n*n*sizeof(double));
 
150
        for (int j=0; j < n; j++) {
 
151
            for (int k=0; k < n; k++) {
 
152
                Q[2*n*j+k]       = a[k];
 
153
                Q[2*n*j+n+k]     = b[k];
 
154
                Q[2*n*(n+j)+k]   = c[k];
 
155
                Q[2*n*(n+j)+n+k] = d[k];
 
156
            }
 
157
        }
 
158
        cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans,
 
159
                    n, n, 4*n, 1.0, P, 4*n, P, 4*n, 0.0, PtP1, 4*n);
 
160
        cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
 
161
                    4*n, 4*n, 4*n, 1.0, PtP1, 4*n, PtP1, 4*n, 0.0, PtP2, 4*n);
 
162
        cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
 
163
                    4*n, 4*n, 4*n, 1.0, PtP2, 4*n, PtP2, 4*n, 0.0, PtP1, 4*n);
 
164
        for (int j=0; j < n; j++) {
 
165
            v[i] += PtP1[(n+1)*j];
 
166
        }
 
167
        cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans,
 
168
                    2*n, 2*n, 2*n, 1.0, Q, 2*n, Q, 2*n, 0.0, QtQ1, 2*n);
 
169
        cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
 
170
                    2*n, 2*n, 2*n, 1.0, QtQ1, 2*n, QtQ1, 2*n, 0.0, QtQ2, 2*n);
 
171
        cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
 
172
                    2*n, 2*n, 2*n, 1.0, QtQ2, 2*n, QtQ2, 2*n, 0.0, QtQ1, 2*n);
 
173
        for (int j=0; j < 2*n; j++) {
 
174
            w[i] += QtQ1[(2*n+1)*j];
 
175
        }
 
176
    }
 
177
    free(PtP1);
 
178
    free(PtP2);
 
179
    free(QtQ1);
 
180
    free(QtQ2);
 
181
    free(P);
 
182
    free(Q);
 
183
    free(a);
 
184
    free(b);
 
185
    free(c);
 
186
    free(d);
 
187
    double v1=0.0, v2=0.0, w1=0.0, w2=0.0;
 
188
    for (int i=0; i < t; i++) {
 
189
        v1 += v[i]; v2 += v[i]*v[i];
 
190
        w1 += w[i]; w2 += w[i]*w[i];
 
191
    }
 
192
    free(v);
 
193
    free(w);
 
194
    r.s1 = sqrt((t*(t*v2-v1*v1))/((t-1)*v1*v1));
 
195
    r.s2 = sqrt((t*(t*w2-w1*w1))/((t-1)*w1*w1));
 
196
    return r;
 
197
}
 
198
 
 
199
double *randmatmul(int n) {
 
200
    double *A = myrand(n*n);
 
201
    double *B = myrand(n*n);
 
202
    double *C = (double*)malloc(n*n*sizeof(double));
 
203
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
 
204
                n, n, n, 1.0, A, n, B, n, 0.0, C, n);
 
205
    free(A);
 
206
    free(B);
 
207
    return C;
 
208
}
 
209
 
 
210
void printfd(int n) {
 
211
    FILE *f = fopen("/dev/null", "w");
 
212
    long i = 0;
 
213
    for (i = 0; i < n; i++)
 
214
        fprintf(f, "%ld %ld\n", i, i);
 
215
    fclose(f);
 
216
}
 
217
 
 
218
void print_perf(const char *name, double t) {
 
219
    printf("c,%s,%.6f\n", name, t*1000);
 
220
}
 
221
 
 
222
int main() {
 
223
    // Initialize RNG
 
224
    dsfmt_gv_init_gen_rand(0);
 
225
 
 
226
    double t, tmin;
 
227
 
 
228
    // fib(20)
 
229
    assert(fib(20) == 6765);
 
230
    int f = 0;
 
231
    tmin = INFINITY;
 
232
    for (int i=0; i<NITER; ++i) {
 
233
        t = clock_now();
 
234
        f += fib(20);
 
235
        t = clock_now()-t;
 
236
        if (t < tmin) tmin = t;
 
237
    }
 
238
    print_perf("fib", tmin);
 
239
 
 
240
    // parse_bin
 
241
    tmin = INFINITY;
 
242
    for (int i=0; i<NITER; ++i) {
 
243
        t = clock_now();
 
244
        char s[11];
 
245
        for (int k=0; k<1000; ++k) {
 
246
            uint32_t n = dsfmt_gv_genrand_uint32();
 
247
            sprintf(s, "%x", n);
 
248
            uint32_t m = (uint32_t)parse_int(s, 16);
 
249
            assert(m == n);
 
250
        }
 
251
        t = clock_now()-t;
 
252
        if (t < tmin) tmin = t;
 
253
    }
 
254
    print_perf("parse_int", tmin);
 
255
 
 
256
    // // array constructor
 
257
    // tmin = INFINITY;
 
258
    // for (int i=0; i<NITER; ++i) {
 
259
    //     t = clock_now();
 
260
    //     double *a = ones(200,200);
 
261
    //     free(a);
 
262
    //     t = clock_now()-t;
 
263
    //     if (t < tmin) tmin = t;
 
264
    // }
 
265
    // print_perf("ones", tmin);
 
266
    // 
 
267
    // // A*A'
 
268
    // //SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
 
269
    // double *b = ones(200, 200);
 
270
    // tmin = INFINITY;
 
271
    // for (int i=0; i<NITER; ++i) {
 
272
    //     t = clock_now();
 
273
    //     double *c = matmul_aat(200, b);
 
274
    //     free(c);
 
275
    //     t = clock_now()-t;
 
276
    //     if (t < tmin) tmin = t;
 
277
    // }
 
278
    // free(b);
 
279
    // print_perf("AtA", tmin);
 
280
 
 
281
    // mandel
 
282
    int mandel_sum;
 
283
    tmin = INFINITY;
 
284
    for (int i=0; i<NITER; ++i) {
 
285
        t = clock_now();
 
286
        mandel_sum = mandelperf();
 
287
        t = clock_now()-t;
 
288
        if (t < tmin) tmin = t;
 
289
    }
 
290
    assert(mandel_sum == 14719);
 
291
    print_perf("mandel", tmin);
 
292
 
 
293
    // sort
 
294
    tmin = INFINITY;
 
295
    for (int i=0; i<NITER; ++i) {
 
296
        t = clock_now();
 
297
        double *d = myrand(5000);
 
298
        quicksort(d, 0, 5000-1);
 
299
        free(d);
 
300
        t = clock_now()-t;
 
301
        if (t < tmin) tmin = t;
 
302
    }
 
303
    print_perf("quicksort", tmin);
 
304
 
 
305
    // pi sum
 
306
    double pi;
 
307
    tmin = INFINITY;
 
308
    for (int i=0; i<NITER; ++i) {
 
309
        t = clock_now();
 
310
        pi = pisum();
 
311
        t = clock_now()-t;
 
312
        if (t < tmin) tmin = t;
 
313
    }
 
314
    assert(fabs(pi-1.644834071848065) < 1e-12);
 
315
    print_perf("pi_sum", tmin);
 
316
 
 
317
    // rand mat stat
 
318
    struct double_pair r;
 
319
    tmin = INFINITY;
 
320
    for (int i=0; i<NITER; ++i) {
 
321
        t = clock_now();
 
322
        r = randmatstat(1000);
 
323
        t = clock_now()-t;
 
324
        if (t < tmin) tmin = t;
 
325
    }
 
326
    // assert(0.5 < r.s1 && r.s1 < 1.0 && 0.5 < r.s2 && r.s2 < 1.0);
 
327
    print_perf("rand_mat_stat", tmin);
 
328
 
 
329
    // rand mat mul
 
330
    tmin = INFINITY;
 
331
    for (int i=0; i<NITER; ++i) {
 
332
        t = clock_now();
 
333
        double *C = randmatmul(1000);
 
334
        assert(0 <= C[0]);
 
335
        free(C);
 
336
        t = clock_now()-t;
 
337
        if (t < tmin) tmin = t;
 
338
    }
 
339
    print_perf("rand_mat_mul", tmin);
 
340
 
 
341
    // printfd
 
342
    tmin = INFINITY;
 
343
    for (int i=0; i<NITER; ++i) {
 
344
        t = clock_now();
 
345
        printfd(100000);
 
346
        t = clock_now()-t;
 
347
        if (t < tmin) tmin = t;
 
348
    }
 
349
    print_perf("printfd", tmin);
 
350
 
 
351
    return 0;
 
352
}