ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rholo3.c
Revision: 3.28
Committed: Sat Jan 9 09:17:10 1999 UTC (25 years, 3 months ago) by gwlarson
Content type: text/plain
Branch: MAIN
Changes since 3.27: +17 -8 lines
Log Message:
added random number seed to init_global()

File Contents

# Content
1 /* Copyright (c) 1998 Silicon Graphics, Inc. */
2
3 #ifndef lint
4 static char SCCSid[] = "$SunId$ SGI";
5 #endif
6
7 /*
8 * Routines for tracking beam compuatations
9 */
10
11 #include "rholo.h"
12 #include <sys/types.h>
13
14 #ifndef NFRAG2CHUNK
15 #define NFRAG2CHUNK 4096 /* number of fragments to start chunking */
16 #endif
17
18 #ifndef abs
19 #define abs(x) ((x) > 0 ? (x) : -(x))
20 #endif
21 #ifndef sgn
22 #define sgn(x) ((x) > 0 ? 1 : (x) < 0 ? -1 : 0)
23 #endif
24
25 #define rchunk(n) (((n)+(RPACKSIZ/2))/RPACKSIZ)
26
27 extern time_t time();
28
29 static PACKHEAD *complist=NULL; /* list of beams to compute */
30 static int complen=0; /* length of complist */
31 static int listpos=0; /* current list position for next_packet */
32 static int lastin= -1; /* last ordered position in list */
33 static int chunky=0; /* clump beams together on disk */
34
35
36 int
37 beamcmp(b0, b1) /* comparison for compute order */
38 register PACKHEAD *b0, *b1;
39 {
40 BEAMI *bip0, *bip1;
41 register long c;
42 /* first check desired quantities */
43 if (chunky)
44 c = rchunk(b1->nr)*(rchunk(b0->nc)+1L) -
45 rchunk(b0->nr)*(rchunk(b1->nc)+1L);
46 else
47 c = b1->nr*(b0->nc+1L) - b0->nr*(b1->nc+1L);
48 if (c > 0) return(1);
49 if (c < 0) return(-1);
50 /* only one file, so skip the following: */
51 #if 0
52 /* next, check file descriptors */
53 c = hdlist[b0->hd]->fd - hdlist[b1->hd]->fd;
54 if (c) return(c);
55 #endif
56 /* finally, check file positions */
57 bip0 = &hdlist[b0->hd]->bi[b0->bi];
58 bip1 = &hdlist[b1->hd]->bi[b1->bi];
59 /* put diskless beams last */
60 if (!bip0->nrd)
61 return(bip1->nrd > 0);
62 if (!bip1->nrd)
63 return(-1);
64 c = bip0->fo - bip1->fo;
65 return(c < 0 ? -1 : c > 0);
66 }
67
68
69 int
70 beamidcmp(b0, b1) /* comparison for beam searching */
71 register PACKHEAD *b0, *b1;
72 {
73 register int c = b0->hd - b1->hd;
74
75 if (c) return(c);
76 return(b0->bi - b1->bi);
77 }
78
79
80 int
81 dispbeam(b, hb) /* display a holodeck beam */
82 register BEAM *b;
83 register HDBEAMI *hb;
84 {
85 static int n = 0;
86 static PACKHEAD *p = NULL;
87
88 if (b == NULL)
89 return;
90 if (b->nrm > n) { /* (re)allocate packet holder */
91 n = b->nrm;
92 if (p == NULL) p = (PACKHEAD *)malloc(packsiz(n));
93 else p = (PACKHEAD *)realloc((char *)p, packsiz(n));
94 if (p == NULL)
95 error(SYSTEM, "out of memory in dispbeam");
96 }
97 /* assign packet fields */
98 bcopy((char *)hdbray(b), (char *)packra(p), b->nrm*sizeof(RAYVAL));
99 p->nr = p->nc = b->nrm;
100 for (p->hd = 0; hdlist[p->hd] != hb->h; p->hd++)
101 if (hdlist[p->hd] == NULL)
102 error(CONSISTENCY, "unregistered holodeck in dispbeam");
103 p->bi = hb->b;
104 disp_packet(p); /* display it */
105 }
106
107
108 bundle_set(op, clist, nents) /* bundle set operation */
109 int op;
110 PACKHEAD *clist;
111 int nents;
112 {
113 int oldnr, n;
114 HDBEAMI *hbarr;
115 register PACKHEAD *csm;
116 register int i;
117 /* search for common members */
118 for (csm = clist+nents; csm-- > clist; )
119 csm->nc = -1;
120 qsort((char *)clist, nents, sizeof(PACKHEAD), beamidcmp);
121 for (i = 0; i < complen; i++) {
122 csm = (PACKHEAD *)bsearch((char *)(complist+i), (char *)clist,
123 nents, sizeof(PACKHEAD), beamidcmp);
124 if (csm == NULL)
125 continue;
126 oldnr = complist[i].nr;
127 csm->nc = complist[i].nc;
128 switch (op) {
129 case BS_ADD: /* add to count */
130 complist[i].nr += csm->nr;
131 csm->nr = 0;
132 break;
133 case BS_ADJ: /* reset count */
134 complist[i].nr = csm->nr;
135 csm->nr = 0;
136 break;
137 case BS_DEL: /* delete count */
138 if (csm->nr == 0 || csm->nr >= complist[i].nr)
139 complist[i].nr = 0;
140 else
141 complist[i].nr -= csm->nr;
142 break;
143 }
144 if (complist[i].nr != oldnr)
145 lastin = -1; /* flag sort */
146 }
147 /* record computed rays for uncommon beams */
148 for (csm = clist+nents; csm-- > clist; )
149 if (csm->nc < 0)
150 csm->nc = bnrays(hdlist[csm->hd], csm->bi);
151 /* complete list operations */
152 switch (op) {
153 case BS_NEW: /* new computation set */
154 listpos = 0; lastin = -1;
155 if (complen) /* free old list */
156 free((char *)complist);
157 complist = NULL;
158 if (!(complen = nents))
159 return;
160 complist = (PACKHEAD *)malloc(nents*sizeof(PACKHEAD));
161 if (complist == NULL)
162 goto memerr;
163 bcopy((char *)clist, (char *)complist, nents*sizeof(PACKHEAD));
164 break;
165 case BS_ADD: /* add to computation set */
166 case BS_ADJ: /* adjust set quantities */
167 if (nents <= 0)
168 return;
169 sortcomplist(); /* sort updated list & new entries */
170 qsort((char *)clist, nents, sizeof(PACKHEAD), beamcmp);
171 /* what can't we satisfy? */
172 for (i = nents, csm = clist; i-- && csm->nr > csm->nc; csm++)
173 ;
174 n = csm - clist;
175 if (op == BS_ADJ) { /* don't regenerate adjusted beams */
176 for (++i; i-- && csm->nr > 0; csm++)
177 ;
178 nents = csm - clist;
179 }
180 if (n) { /* allocate space for merged list */
181 PACKHEAD *newlist;
182 newlist = (PACKHEAD *)malloc(
183 (complen+n)*sizeof(PACKHEAD) );
184 if (newlist == NULL)
185 goto memerr;
186 /* merge lists */
187 mergeclists(newlist, clist, n, complist, complen);
188 if (complen)
189 free((char *)complist);
190 complist = newlist;
191 complen += n;
192 }
193 listpos = 0;
194 lastin = complen-1; /* list is now sorted */
195 break;
196 case BS_DEL: /* delete from computation set */
197 return; /* already done */
198 default:
199 error(CONSISTENCY, "bundle_set called with unknown operation");
200 }
201 if (outdev == NULL || !nents) /* nothing to display? */
202 return;
203 /* load and display beams we have */
204 hbarr = (HDBEAMI *)malloc(nents*sizeof(HDBEAMI));
205 for (i = nents; i--; ) {
206 hbarr[i].h = hdlist[clist[i].hd];
207 hbarr[i].b = clist[i].bi;
208 }
209 hdloadbeams(hbarr, nents, dispbeam);
210 free((char *)hbarr);
211 if (hdfragflags&FF_READ) {
212 listpos = 0;
213 lastin = -1; /* need to re-sort list */
214 }
215 return;
216 memerr:
217 error(SYSTEM, "out of memory in bundle_set");
218 }
219
220
221 double
222 beamvolume(hp, bi) /* compute approximate volume of a beam */
223 HOLO *hp;
224 int bi;
225 {
226 GCOORD gc[2];
227 FVECT cp[4], edgeA, edgeB, cent[2];
228 FVECT v, crossp[2], diffv;
229 double vol[2];
230 register int i;
231 /* get grid coordinates */
232 if (!hdbcoord(gc, hp, bi))
233 error(CONSISTENCY, "bad beam index in beamvolume");
234 for (i = 0; i < 2; i++) { /* compute cell area vectors */
235 hdcell(cp, hp, gc+i);
236 VSUM(edgeA, cp[1], cp[0], -1.0);
237 VSUM(edgeB, cp[3], cp[1], -1.0);
238 fcross(crossp[i], edgeA, edgeB);
239 /* compute center */
240 cent[i][0] = 0.5*(cp[0][0] + cp[2][0]);
241 cent[i][1] = 0.5*(cp[0][1] + cp[2][1]);
242 cent[i][2] = 0.5*(cp[0][2] + cp[2][2]);
243 }
244 /* compute difference vector */
245 VSUM(diffv, cent[1], cent[0], -1.0);
246 for (i = 0; i < 2; i++) { /* compute volume contributions */
247 vol[i] = 0.5*DOT(crossp[i], diffv);
248 if (vol[i] < 0.) vol[i] = -vol[i];
249 }
250 return(vol[0] + vol[1]); /* return total volume */
251 }
252
253
254 init_global() /* initialize global ray computation */
255 {
256 long wtotal = 0;
257 double frac;
258 int i;
259 register int j, k;
260 /* free old list and empty queue */
261 if (complen > 0) {
262 free((char *)complist);
263 done_packets(flush_queue());
264 }
265 /* reseed random number generator */
266 srandom(time(NULL));
267 /* allocate beam list */
268 complen = 0;
269 for (j = 0; hdlist[j] != NULL; j++)
270 complen += nbeams(hdlist[j]);
271 complist = (PACKHEAD *)malloc(complen*sizeof(PACKHEAD));
272 if (complist == NULL)
273 error(SYSTEM, "out of memory in init_global");
274 /* compute beam weights */
275 k = 0;
276 for (j = 0; hdlist[j] != NULL; j++) {
277 frac = 512. * VLEN(hdlist[j]->wg[0]) *
278 VLEN(hdlist[j]->wg[1]) *
279 VLEN(hdlist[j]->wg[2]);
280 for (i = nbeams(hdlist[j]); i > 0; i--) {
281 complist[k].hd = j;
282 complist[k].bi = i;
283 complist[k].nr = frac*beamvolume(hdlist[j], i) + 0.5;
284 complist[k].nc = bnrays(hdlist[j], i);
285 wtotal += complist[k++].nr;
286 }
287 }
288 /* adjust weights */
289 if (vdef(DISKSPACE))
290 frac = 1024.*1024.*vflt(DISKSPACE) / (wtotal*sizeof(RAYVAL));
291 else
292 frac = 1024.*1024.*16384. / (wtotal*sizeof(RAYVAL));
293 while (k--)
294 complist[k].nr = frac*complist[k].nr + 0.5;
295 listpos = 0; lastin = -1; /* perform initial sort */
296 sortcomplist();
297 /* no view vicinity */
298 myeye.rng = 0;
299 }
300
301
302 mergeclists(cdest, cl1, n1, cl2, n2) /* merge two sorted lists */
303 register PACKHEAD *cdest;
304 register PACKHEAD *cl1, *cl2;
305 int n1, n2;
306 {
307 register int cmp;
308
309 while (n1 | n2) {
310 if (!n1) cmp = 1;
311 else if (!n2) cmp = -1;
312 else cmp = beamcmp(cl1, cl2);
313 if (cmp > 0) {
314 copystruct(cdest, cl2);
315 cl2++; n2--;
316 } else {
317 copystruct(cdest, cl1);
318 cl1++; n1--;
319 }
320 cdest++;
321 }
322 }
323
324
325 sortcomplist() /* fix our list order */
326 {
327 PACKHEAD *list2;
328 int listlen;
329 register int i;
330
331 if (complen <= 0) /* check to see if there is even a list */
332 return;
333 if (!chunky) /* check to see if fragment list is full */
334 if (!hdfragOK(hdlist[0]->fd, &listlen, NULL)
335 #if NFRAG2CHUNK
336 || listlen >= NFRAG2CHUNK
337 #endif
338 ) {
339 #ifdef DEBUG
340 error(WARNING, "using chunky comparison mode");
341 #endif
342 chunky++; /* use "chunky" comparison */
343 lastin = -1; /* need to re-sort list */
344 }
345 #ifdef DEBUG
346 else
347 fprintf(stderr, "sortcomplist: %d fragments\n",
348 listlen);
349 #endif
350 if (lastin < 0 || listpos*4 >= complen*3)
351 qsort((char *)complist, complen, sizeof(PACKHEAD), beamcmp);
352 else if (listpos) { /* else sort and merge sublist */
353 list2 = (PACKHEAD *)malloc(listpos*sizeof(PACKHEAD));
354 if (list2 == NULL)
355 error(SYSTEM, "out of memory in sortcomplist");
356 bcopy((char *)complist,(char *)list2,listpos*sizeof(PACKHEAD));
357 qsort((char *)list2, listpos, sizeof(PACKHEAD), beamcmp);
358 mergeclists(complist, list2, listpos,
359 complist+listpos, complen-listpos);
360 free((char *)list2);
361 }
362 /* drop satisfied requests */
363 for (i = complen; i-- && complist[i].nr <= complist[i].nc; )
364 ;
365 if (i < 0) {
366 free((char *)complist);
367 complist = NULL;
368 complen = 0;
369 } else if (i < complen-1) {
370 list2 = (PACKHEAD *)realloc((char *)complist,
371 (i+1)*sizeof(PACKHEAD));
372 if (list2 != NULL)
373 complist = list2;
374 complen = i+1;
375 }
376 listpos = 0; lastin = i;
377 }
378
379
380 /*
381 * The following routine works on the assumption that the bundle weights are
382 * more or less evenly distributed, such that computing a packet causes
383 * a given bundle to move way down in the computation order. We keep
384 * track of where the computed bundle with the highest priority would end
385 * up, and if we get further in our compute list than this, we re-sort the
386 * list and start again from the beginning. Since
387 * a merge sort is used, the sorting costs are minimal.
388 */
389 next_packet(p, n) /* prepare packet for computation */
390 register PACKET *p;
391 int n;
392 {
393 register int i;
394
395 if (listpos > lastin) /* time to sort the list */
396 sortcomplist();
397 if (complen <= 0)
398 return(0);
399 p->hd = complist[listpos].hd;
400 p->bi = complist[listpos].bi;
401 p->nc = complist[listpos].nc;
402 p->nr = complist[listpos].nr - p->nc;
403 if (p->nr <= 0)
404 return(0);
405 DCHECK(n < 1 | n > RPACKSIZ,
406 CONSISTENCY, "next_packet called with bad n value");
407 if (p->nr > n)
408 p->nr = n;
409 complist[listpos].nc += p->nr; /* find where this one would go */
410 if (hdgetbeam(hdlist[p->hd], p->bi) != NULL)
411 hdfreefrag(hdlist[p->hd], p->bi);
412 while (lastin > listpos &&
413 beamcmp(complist+lastin, complist+listpos) > 0)
414 lastin--;
415 listpos++;
416 return(1);
417 }