ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/rholo3.c
Revision: 3.43
Committed: Thu Sep 30 14:10:19 2010 UTC (13 years, 6 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R2P2, rad5R0, rad5R1, rad4R2, rad4R1, rad4R2P1
Changes since 3.42: +7 -2 lines
Log Message:
Increased maximum "ambient" holodeck size from 2 GBytes to 10 GBytes

File Contents

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