← Back to branch summary

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

« back to all changes in this revision

Viewing changes to test/perf/perf.cpp

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