1 /*
2 Copyright (C) 2021 The Blosc Developers <blosc@blosc.org>
3 https://blosc.org
4 License: BSD 3-Clause (see LICENSE.txt)
5
6 Example program demonstrating how the different compression params affects
7 the performance of root finding.
8
9 To compile this program:
10
11 $ gcc -O3 find_roots.c -o find_roots -lblosc2
12
13 To run:
14
15 $ ./find_roots
16 Blosc version info: 2.0.0a6.dev ($Date:: 2018-05-18 #$)
17 Creation time for X values: 0.178 s, 4274.5 MB/s
18 Compression for X values: 762.9 MB -> 27.3 MB (28.0x)
19 Computing Y polynomial: 0.342 s, 4463.3 MB/s
20 Compression for Y values: 762.9 MB -> 54.0 MB (14.1x)
21 Roots found at: 1.350000023841858, 4.450000286102295, 8.500000953674316,
22 Find root time: 0.401 s, 3806.8 MB/s
23
24 */
25
26 #include <stdio.h>
27 #include "blosc2.h"
28
29 #define KB 1024.
30 #define MB (1024*KB)
31 #define GB (1024*MB)
32
33
34 #define NCHUNKS 500
35 #define CHUNKSIZE (200 * 1000) // fits well in modern L3 caches
36 #define NTHREADS 4
37
38
fill_buffer(double * x,int nchunk)39 void fill_buffer(double *x, int nchunk) {
40 double incx = 10. / (NCHUNKS * CHUNKSIZE);
41
42 for (int i = 0; i < CHUNKSIZE; i++) {
43 x[i] = incx * (nchunk * CHUNKSIZE + i);
44 }
45 }
46
process_data(const double * x,double * y)47 void process_data(const double *x, double *y) {
48
49 for (int i = 0; i < CHUNKSIZE; i++) {
50 double xi = x[i];
51 //y[i] = ((.25 * xi + .75) * xi - 1.5) * xi - 2;
52 y[i] = (xi - 1.35) * (xi - 4.45) * (xi - 8.5);
53 }
54 }
55
find_root(const double * x,const double * y,const double prev_value)56 void find_root(const double *x, const double *y,
57 const double prev_value) {
58 double pv = prev_value;
59 int last_root_idx = -1;
60
61 for (int i = 0; i < CHUNKSIZE; i++) {
62 double yi = y[i];
63 if (((yi > 0) - (yi < 0)) != ((pv > 0) - (pv < 0))) {
64 if (last_root_idx != (i - 1)) {
65 printf("%.16g, ", x[i]);
66 last_root_idx = i; // avoid the last point (ULP effects)
67 }
68 }
69 pv = yi;
70 }
71 }
72
73
compute_vectors(void)74 int compute_vectors(void) {
75 static double buffer_x[CHUNKSIZE];
76 static double buffer_y[CHUNKSIZE];
77 const int32_t isize = CHUNKSIZE * sizeof(double);
78 int dsize;
79 long nbytes = 0;
80 blosc2_schunk *sc_x, *sc_y;
81 int nchunk;
82 blosc_timestamp_t last, current;
83 double ttotal;
84 double prev_value;
85
86 /* Create a super-chunk container for input (X values) */
87 blosc2_cparams cparams = BLOSC2_CPARAMS_DEFAULTS;
88 cparams.typesize = sizeof(double);
89 cparams.compcode = BLOSC_LZ4;
90 cparams.clevel = 9;
91 cparams.filters[0] = BLOSC_TRUNC_PREC;
92 cparams.filters_meta[0] = 23; // treat doubles as floats
93 cparams.nthreads = NTHREADS;
94 blosc2_dparams dparams = BLOSC2_DPARAMS_DEFAULTS;
95 dparams.nthreads = NTHREADS;
96 blosc2_storage storage = {.cparams=&cparams, .dparams=&dparams};
97 sc_x = blosc2_schunk_new(&storage);
98
99 /* Create a super-chunk container for output (Y values) */
100 sc_y = blosc2_schunk_new(&storage);
101
102 /* Now fill the buffer with even values between 0 and 10 */
103 blosc_set_timestamp(&last);
104 for (nchunk = 0; nchunk < NCHUNKS; nchunk++) {
105 fill_buffer(buffer_x, nchunk);
106 blosc2_schunk_append_buffer(sc_x, buffer_x, isize);
107 nbytes += (long) isize;
108 }
109 blosc_set_timestamp(¤t);
110 ttotal = blosc_elapsed_secs(last, current);
111 printf("Creation time for X values: %.3g s, %.1f MB/s\n",
112 ttotal, (double) nbytes / (ttotal * MB));
113 printf("Compression for X values: %.1f MB -> %.1f MB (%.1fx)\n",
114 sc_x->nbytes / MB, sc_x->cbytes / MB,
115 (1. * sc_x->nbytes) / sc_x->cbytes);
116
117 /* Retrieve the chunks and compute the polynomial in another super-chunk */
118 blosc_set_timestamp(&last);
119 for (nchunk = 0; nchunk < NCHUNKS; nchunk++) {
120 dsize = blosc2_schunk_decompress_chunk(sc_x, nchunk, buffer_x, isize);
121 if (dsize < 0) {
122 printf("Decompression error. Error code: %d\n", dsize);
123 return dsize;
124 }
125 process_data(buffer_x, buffer_y);
126 blosc2_schunk_append_buffer(sc_y, buffer_y, isize);
127 }
128 blosc_set_timestamp(¤t);
129 ttotal = blosc_elapsed_secs(last, current);
130 printf("Computing Y polynomial: %.3g s, %.1f MB/s\n",
131 ttotal, 2. * (double) nbytes / (ttotal * MB)); // 2 super-chunks involved
132 printf("Compression for Y values: %.1f MB -> %.1f MB (%.1fx)\n",
133 sc_y->nbytes / MB, sc_y->cbytes / MB,
134 (1. * sc_y->nbytes) / sc_y->cbytes);
135
136 /* Find the roots of the polynomial */
137 printf("Roots found at: ");
138 blosc_set_timestamp(&last);
139 prev_value = buffer_y[0];
140 for (nchunk = 0; nchunk < NCHUNKS; nchunk++) {
141 dsize = blosc2_schunk_decompress_chunk(sc_y, nchunk, (void *) buffer_y, isize);
142 if (dsize < 0) {
143 printf("Decompression error. Error code: %d\n", dsize);
144 return dsize;
145 }
146 dsize = blosc2_schunk_decompress_chunk(sc_x, nchunk, (void *) buffer_x, isize);
147 if (dsize < 0) {
148 printf("Decompression error. Error code: %d\n", dsize);
149 return dsize;
150 }
151 find_root(buffer_x, buffer_y, prev_value);
152 prev_value = buffer_y[CHUNKSIZE - 1];
153 }
154 blosc_set_timestamp(¤t);
155 ttotal = blosc_elapsed_secs(last, current);
156 printf("\n");
157 printf("Find root time: %.3g s, %.1f MB/s\n",
158 ttotal, 2. * (double) nbytes / (ttotal * MB)); // 2 super-chunks involved
159
160 /* Free resources */
161 /* Destroy the super-chunk */
162 blosc2_schunk_free(sc_x);
163 blosc2_schunk_free(sc_y);
164 return 0;
165 }
166
167
main(void)168 int main(void) {
169 printf("Blosc version info: %s (%s)\n",
170 BLOSC_VERSION_STRING, BLOSC_VERSION_DATE);
171
172 /* Initialize the Blosc compressor */
173 blosc_init();
174
175 compute_vectors();
176
177 /* Destroy the Blosc environment */
178 blosc_destroy();
179
180 return 0;
181 }
182