Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] Annotation of /benchmarks/programs/lic2d/bmark-hact.c
ViewVC logotype

Annotation of /benchmarks/programs/lic2d/bmark-hact.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1615 - (view) (download) (as text)

1 : glk 1615 /*! \file /lic2d/bmark-teem.c
2 :     *
3 :     * \author Nick Seltzer and Gordon Kindlmann
4 :     */
5 :    
6 :     /*
7 :     * COPYRIGHT (c) 2011 The Diderot Project (http://diderot-language.cs.uchicago.edu)
8 :     * All rights reserved.
9 :     */
10 :    
11 :     #include <stdio.h>
12 :     #include <math.h>
13 :     #include <stdbool.h>
14 :     #include "teem/nrrd.h"
15 :     #include "teem/ell.h"
16 :     #include "teem/gage.h"
17 :    
18 :     #include "teem-defs.h"
19 :    
20 :     STATIC_INLINE void scale2(double scl, double vec[2], double res[2])
21 :     {
22 :     res[0] = scl * vec[0];
23 :     res[1] = scl * vec[1];
24 :     }
25 :    
26 :     STATIC_INLINE void add2(double vec1[2], double vec2[2], double res[2])
27 :     {
28 :     res[0] = vec1[0] + vec2[0];
29 :     res[1] = vec1[1] + vec2[1];
30 :     }
31 :    
32 :     STATIC_INLINE double dot2(double vec1[2], double vec2[2])
33 :     {
34 :     return (vec1[0] * vec2[0] + vec1[1] * vec2[1]);
35 :     }
36 :    
37 :     STATIC_INLINE double mag2(double vec[2])
38 :     {
39 :     return sqrt(dot2(vec, vec));
40 :     }
41 :    
42 :     STATIC_INLINE void norm2(double vec[2], double res[2])
43 :     {
44 :     double mag = mag2(vec);
45 :     res[0] = vec[0] / mag;
46 :     res[1] = vec[1] / mag;
47 :     }
48 :    
49 :     STATIC_INLINE bool inside(double pos[2])
50 :     {
51 :     // XXX - Hack
52 :     if (pos[0] < 0.00 || pos[0] >= 6.78) return false;
53 :     if (pos[1] < 0.00 || pos[1] >= 3.72) return false;
54 :     return true;
55 :     }
56 :    
57 :     void atR(void *_data, double pos[2], double *out)
58 :     {
59 :     float *data = (float*)_data;
60 :     double temp[2][2];
61 :     int a = floor(pos[0] * 150) + 1;
62 :     int b = floor(pos[1] * 150) + 1;
63 :     for (int i = 0; i < 2; i++) {
64 :     for (int j = 0; j < 2; j++) {
65 :     temp[i][j] = *(data + (i + a) + (j + b) * 1020);
66 :     }
67 :     }
68 :     double temp2[2];
69 :     double t = pos[0]*150 - floor(pos[0]*150);
70 :     for (int i = 0; i < 2; i++) {
71 :     temp2[i] = (temp[0][i] * (1 - t) + temp[1][i] * t);
72 :     }
73 :     t = pos[1]*150 - floor(pos[1]*150);
74 :     *out = (temp2[0] * (1 - t) + temp2[1] * t);
75 :     }
76 :    
77 :     void at(void* _data, double pos[2], double out[2])
78 :     {
79 :     // XXX - Hack
80 :     float *data = (float*)_data;
81 :     double tempx[4][4];
82 :     double tempy[4][4];
83 :     int a = floor(pos[0] + 3);
84 :     int b = floor(pos[1] + 3);
85 :     for (int i = 0; i < 4; i++) {
86 :     for (int j = 0; j < 4; j++) {
87 :     tempx[i][j] = data[0 + 2*(i + a - 1 + 15*(j + b - 1))];
88 :     tempy[i][j] = data[1 + 2*(i + a - 1 + 15*(j + b - 1))];
89 :     }
90 :     }
91 :     double t = pos[0]+3 - a;
92 :     double tempxx[4];
93 :     double tempyx[4];
94 :     for (int i = 0; i < 4; i++) {
95 :     tempxx[i] = (tempx[0][i] + 4 * tempx[1][i] + tempx[2][i] +
96 :     (-3 * tempx[0][i] + 3 * tempx[2][i]) * t +
97 :     (3 * tempx[0][i] - 6 * tempx[1][i] + 3 * tempx[2][i]) * t * t +
98 :     (-1 * tempx[0][i] + 3 * tempx[1][i] - 3 * tempx[2][i] + tempx[3][i]) * t * t * t) / 6.0;
99 :     tempyx[i] = (tempy[0][i] + 4 * tempy[1][i] + tempy[2][i] +
100 :     (-3 * tempy[0][i] + 3 * tempy[2][i]) * t +
101 :     (3 * tempy[0][i] - 6 * tempy[1][i] + 3 * tempy[2][i]) * t * t +
102 :     (-1 * tempy[0][i] + 3 * tempy[1][i] - 3 * tempy[2][i] + tempy[3][i]) * t * t * t) / 6.0;
103 :     }
104 :     t = pos[1]+3 - b;
105 :     out[0] = (tempxx[0] + 4 * tempxx[1] + tempxx[2] +
106 :     (-3 * tempxx[0] + 3 * tempxx[2]) * t +
107 :     (3 * tempxx[0] - 6 * tempxx[1] + 3 * tempxx[2]) * t * t +
108 :     (-1 * tempxx[0] + 3 * tempxx[1] - 3 * tempxx[2] + tempxx[3]) * t * t * t) / 6.0;
109 :     out[1] = (tempyx[0] + 4 * tempyx[1] + tempyx[2] +
110 :     (-3 * tempyx[0] + 3 * tempyx[2]) * t +
111 :     (3 * tempyx[0] - 6 * tempyx[1] + 3 * tempyx[2]) * t * t +
112 :     (-1 * tempyx[0] + 3 * tempyx[1] - 3 * tempyx[2] + tempyx[3]) * t * t * t) / 6.0;
113 :     return;
114 :     }
115 :    
116 :     double crl_at(void* _data, double pos[2])
117 :     {
118 :     // XXX - Hack
119 :     float *data = (float*)_data;
120 :     double tempx[4][4];
121 :     double tempy[4][4];
122 :     double temp[2];
123 :     int a = floor(pos[0] * 50) + 5;
124 :     int b = floor(pos[1] * 50) + 5;
125 :     for (int i = 0; i < 4; i++) {
126 :     for (int j = 0; j < 4; j++) {
127 :     tempx[i][j] = *(data + (j + b - 1) * 700 + (i + a - 1) * 2 + 0);
128 :     tempy[i][j] = *(data + (j + b - 1) * 700 + (i + a - 1) * 2 + 1);
129 :     }
130 :     }
131 :     double t = pos[0]*50 - floor(pos[0]*50);
132 :     double tempxx[4];
133 :     double tempyx[4];
134 :     for (int i = 0; i < 4; i++) {
135 :     tempxx[i] = (tempx[0][i] + 4 * tempx[1][i] + tempx[2][i] +
136 :     (-3 * tempx[0][i] + 3 * tempx[2][i]) * t +
137 :     (3 * tempx[0][i] - 6 * tempx[1][i] + 3 * tempx[2][i]) * t * t +
138 :     (-1 * tempx[0][i] + 3 * tempx[1][i] - 3 * tempx[2][i] + tempx[3][i]) * t * t * t) / 6.0;
139 :     tempyx[i] = ((-3 * tempy[0][i] + 3 * tempy[2][i]) * 1 +
140 :     (3 * tempy[0][i] - 6 * tempy[1][i] + 3 * tempy[2][i]) * t * 2 +
141 :     (-1 * tempy[0][i] + 3 * tempy[1][i] - 3 * tempy[2][i] + tempy[3][i]) * t * t * 3) / 6.0;
142 :     }
143 :     t = pos[1]*50 - floor(pos[1]*50);
144 :     temp[0] = ((-3 * tempxx[0] + 3 * tempxx[2]) * 1 +
145 :     (3 * tempxx[0] - 6 * tempxx[1] + 3 * tempxx[2]) * t * 2 +
146 :     (-1 * tempxx[0] + 3 * tempxx[1] - 3 * tempxx[2] + tempxx[3]) * t * t * 3) / 6.0;
147 :     temp[1] = (tempyx[0] + 4 * tempyx[1] + tempyx[2] +
148 :     (-3 * tempyx[0] + 3 * tempyx[2]) * t +
149 :     (3 * tempyx[0] - 6 * tempyx[1] + 3 * tempyx[2]) * t * t +
150 :     (-1 * tempyx[0] + 3 * tempyx[1] - 3 * tempyx[2] + tempyx[3]) * t * t * t) / 6.0;
151 :     return temp[1] - temp[0];
152 :     }
153 :    
154 :     STATIC_INLINE double lerp(double out_min, double out_max, double in_min, double in, double in_max)
155 :     {
156 :     return (out_min * (in_max - in) + out_max * (in - in_min)) / (in_max - in_min);
157 :     }
158 :    
159 :     STATIC_INLINE unsigned char ulerp(double in_min, double in, double in_max)
160 :     {
161 :     return (unsigned char)((0.0f * (in_max - in) + 255.0f * (in - in_min)) / (in_max - in_min));
162 :     }
163 :    
164 :     int main (int argc, const char *argv[])
165 :     {
166 :     const char *me;
167 :     char *inFile = "imp-u.nrrd";
168 :     char *randFile = "../../data/R.nrrd";
169 :     char *outFile = "bmark-teem.nrrd";
170 :     int imgSizeX = 500;
171 :     int imgSizeY = 500;
172 :     double h = 0.005;
173 :     int stepNum = 10;
174 :    
175 :     me = argv[0];
176 :     airArray *mop;
177 :     mop = airMopNew();
178 :    
179 :     // Read in the data
180 :     char* err;
181 :     Nrrd *V = nrrdNew();
182 :     airMopAdd(mop, V, (airMopper)nrrdNuke, airMopAlways);
183 :     Nrrd *R = nrrdNew();
184 :     airMopAdd(mop, R, (airMopper)nrrdNuke, airMopAlways);
185 :     if (nrrdLoad(V, inFile, NULL)
186 :     || nrrdLoad(R, randFile, NULL)) {
187 :     airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
188 :     fprintf(stderr, "%s: Trouble reading:\n%s", me, err);
189 :     airMopError(mop);
190 :     return -1;
191 :     }
192 :    
193 :     Nrrd *nout = nrrdNew();
194 :     airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);
195 :     if (nrrdAlloc_va(nout, nrrdTypeDouble, 3,
196 :     AIR_CAST(size_t, 3),
197 :     AIR_CAST(size_t, 500),
198 :     AIR_CAST(size_t, 500))) {
199 :     airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
200 :     fprintf(stderr, "%s: Trouble allocating:\n%s", me, err);
201 :     airMopError(mop);
202 :     return -1;
203 :     }
204 :     double *out_data = AIR_CAST(double *, nout->data);
205 :    
206 :     double t0 = airTime(); // start timing
207 :     double pos0[2];
208 :     double forw[2];
209 :     double back[2];
210 :     double temp[2];
211 :     double sum;
212 :     double crl;
213 :     double s, l, c, m;
214 :     int step;
215 :     int num;
216 :     for (int yi = 0; yi < 500; yi++) {
217 :     pos0[1] = lerp(0.0, 8.0, 0, yi, 499);
218 :     for (int xi = 0; xi < 500; xi++) {
219 :     pos0[0] = lerp(0.0, 8.0, 0, xi, 499);
220 :     forw[0] = pos0[0];
221 :     forw[1] = pos0[1];
222 :     back[0] = pos0[0];
223 :     back[1] = pos0[1];
224 :     atR(R->data, pos0, &sum);
225 :     sum = -sum;
226 :     step = -1;
227 :     num = -1;
228 :     while (true) {
229 :     at(V->data, pos0, temp);
230 :     ELL_3V_SET(out_data + 3*(xi + 500*yi), temp[0], temp[1], temp[0]);
231 :     break;
232 :     }
233 :     }
234 :     }
235 :    
236 :     double totalTime = airTime() - t0;
237 :     printf("usr=%f\n", totalTime);
238 :    
239 :     if (nrrdSave(outFile, nout, NULL)) {
240 :     airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
241 :     fprintf(stderr, "%s: Trouble saving:\n%s", me, err);
242 :     airMopError(mop);
243 :     return -1;
244 :     }
245 :    
246 :     airMopOkay(mop);
247 :     return 0;
248 :     }

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0