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, 9 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

# User Rev Content
1 gwlarson 3.23 /* Copyright (c) 1998 Silicon Graphics, Inc. */
2    
3     #ifndef lint
4     static char SCCSid[] = "$SunId$ SGI";
5     #endif
6    
7 gregl 3.1 /*
8     * Routines for tracking beam compuatations
9     */
10    
11     #include "rholo.h"
12 gwlarson 3.28 #include <sys/types.h>
13 gregl 3.1
14 gwlarson 3.27 #ifndef NFRAG2CHUNK
15     #define NFRAG2CHUNK 4096 /* number of fragments to start chunking */
16     #endif
17    
18     #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.28 extern time_t time();
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 gwlarson 3.27 static int chunky=0; /* clump beams together on disk */
34 gregl 3.1
35    
36     int
37 gwlarson 3.23 beamcmp(b0, b1) /* comparison for compute order */
38 gregl 3.2 register PACKHEAD *b0, *b1;
39     {
40 gwlarson 3.27 BEAMI *bip0, *bip1;
41     register long c;
42     /* first check desired quantities */
43     if (chunky)
44 gwlarson 3.28 c = rchunk(b1->nr)*(rchunk(b0->nc)+1L) -
45     rchunk(b0->nr)*(rchunk(b1->nc)+1L);
46 gwlarson 3.27 else
47 gwlarson 3.28 c = b1->nr*(b0->nc+1L) - b0->nr*(b1->nc+1L);
48     if (c > 0) return(1);
49     if (c < 0) return(-1);
50 gwlarson 3.27 /* 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 gregl 3.2 }
67    
68    
69 gregl 3.14 int
70 gwlarson 3.23 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 gregl 3.14 register BEAM *b;
83 gwlarson 3.23 register HDBEAMI *hb;
84 gregl 3.14 {
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 gwlarson 3.23 for (p->hd = 0; hdlist[p->hd] != hb->h; p->hd++)
101 gregl 3.14 if (hdlist[p->hd] == NULL)
102     error(CONSISTENCY, "unregistered holodeck in dispbeam");
103 gwlarson 3.23 p->bi = hb->b;
104 gregl 3.14 disp_packet(p); /* display it */
105     }
106    
107    
108 gregl 3.2 bundle_set(op, clist, nents) /* bundle set operation */
109     int op;
110 gwlarson 3.23 PACKHEAD *clist;
111 gregl 3.2 int nents;
112     {
113 gwlarson 3.23 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 gwlarson 3.25 qsort((char *)clist, nents, sizeof(PACKHEAD), beamidcmp);
121 gwlarson 3.23 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 gregl 3.20 }
147 gwlarson 3.24 /* record computed rays for uncommon beams */
148 gwlarson 3.23 for (csm = clist+nents; csm-- > clist; )
149     if (csm->nc < 0)
150     csm->nc = bnrays(hdlist[csm->hd], csm->bi);
151 gregl 3.20 /* complete list operations */
152 gregl 3.2 switch (op) {
153     case BS_NEW: /* new computation set */
154 gwlarson 3.21 listpos = 0; lastin = -1;
155 gregl 3.20 if (complen) /* free old list */
156 gregl 3.2 free((char *)complist);
157 gregl 3.20 complist = NULL;
158     if (!(complen = nents))
159 gregl 3.2 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 gregl 3.11 case BS_ADJ: /* adjust set quantities */
167 gregl 3.2 if (nents <= 0)
168     return;
169 gregl 3.20 sortcomplist(); /* sort updated list & new entries */
170 gregl 3.2 qsort((char *)clist, nents, sizeof(PACKHEAD), beamcmp);
171     /* what can't we satisfy? */
172 gwlarson 3.23 for (i = nents, csm = clist; i-- && csm->nr > csm->nc; csm++)
173 gregl 3.2 ;
174 gwlarson 3.23 n = csm - clist;
175 gwlarson 3.22 if (op == BS_ADJ) { /* don't regenerate adjusted beams */
176 gwlarson 3.24 for (++i; i-- && csm->nr > 0; csm++)
177 gwlarson 3.22 ;
178 gwlarson 3.24 nents = csm - clist;
179 gwlarson 3.22 }
180 gregl 3.2 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 gregl 3.20 return; /* already done */
198 gregl 3.2 default:
199     error(CONSISTENCY, "bundle_set called with unknown operation");
200     }
201 gwlarson 3.27 if (outdev == NULL || !nents) /* nothing to display? */
202 gwlarson 3.22 return;
203     /* load and display beams we have */
204 gwlarson 3.23 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 gregl 3.11 }
209 gwlarson 3.23 hdloadbeams(hbarr, nents, dispbeam);
210     free((char *)hbarr);
211 gwlarson 3.27 if (hdfragflags&FF_READ) {
212     listpos = 0;
213     lastin = -1; /* need to re-sort list */
214     }
215 gregl 3.2 return;
216     memerr:
217     error(SYSTEM, "out of memory in bundle_set");
218     }
219    
220    
221 gregl 3.13 double
222     beamvolume(hp, bi) /* compute approximate volume of a beam */
223 gregl 3.1 HOLO *hp;
224 gregl 3.13 int bi;
225 gregl 3.1 {
226 gregl 3.13 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 gregl 3.1 }
244 gregl 3.13 /* compute difference vector */
245     VSUM(diffv, cent[1], cent[0], -1.0);
246     for (i = 0; i < 2; i++) { /* compute volume contributions */
247 gregl 3.15 vol[i] = 0.5*DOT(crossp[i], diffv);
248 gregl 3.13 if (vol[i] < 0.) vol[i] = -vol[i];
249     }
250     return(vol[0] + vol[1]); /* return total volume */
251 gregl 3.1 }
252    
253    
254     init_global() /* initialize global ray computation */
255     {
256     long wtotal = 0;
257     double frac;
258 gregl 3.13 int i;
259     register int j, k;
260 gregl 3.18 /* free old list and empty queue */
261     if (complen > 0) {
262 gregl 3.3 free((char *)complist);
263 gregl 3.18 done_packets(flush_queue());
264     }
265 gwlarson 3.28 /* reseed random number generator */
266     srandom(time(NULL));
267 gregl 3.1 /* 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 gregl 3.13 for (j = 0; hdlist[j] != NULL; j++) {
277 gregl 3.19 frac = 512. * VLEN(hdlist[j]->wg[0]) *
278     VLEN(hdlist[j]->wg[1]) *
279     VLEN(hdlist[j]->wg[2]);
280 gregl 3.1 for (i = nbeams(hdlist[j]); i > 0; i--) {
281     complist[k].hd = j;
282     complist[k].bi = i;
283 gregl 3.13 complist[k].nr = frac*beamvolume(hdlist[j], i) + 0.5;
284 gregl 3.18 complist[k].nc = bnrays(hdlist[j], i);
285 gregl 3.1 wtotal += complist[k++].nr;
286     }
287 gregl 3.13 }
288 gregl 3.1 /* adjust weights */
289 gregl 3.12 if (vdef(DISKSPACE))
290 gregl 3.1 frac = 1024.*1024.*vflt(DISKSPACE) / (wtotal*sizeof(RAYVAL));
291 gregl 3.12 else
292     frac = 1024.*1024.*16384. / (wtotal*sizeof(RAYVAL));
293     while (k--)
294 gwlarson 3.24 complist[k].nr = frac*complist[k].nr + 0.5;
295 gwlarson 3.21 listpos = 0; lastin = -1; /* perform initial sort */
296     sortcomplist();
297 gwlarson 3.25 /* no view vicinity */
298     myeye.rng = 0;
299 gregl 3.1 }
300    
301    
302     mergeclists(cdest, cl1, n1, cl2, n2) /* merge two sorted lists */
303 gregl 3.16 register PACKHEAD *cdest;
304     register PACKHEAD *cl1, *cl2;
305 gregl 3.1 int n1, n2;
306     {
307 gregl 3.16 register int cmp;
308 gregl 3.1
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 gwlarson 3.27 int listlen;
329 gregl 3.2 register int i;
330    
331 gregl 3.3 if (complen <= 0) /* check to see if there is even a list */
332 gregl 3.2 return;
333 gwlarson 3.27 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 gwlarson 3.28 #ifdef DEBUG
346     else
347     fprintf(stderr, "sortcomplist: %d fragments\n",
348     listlen);
349     #endif
350 gregl 3.6 if (lastin < 0 || listpos*4 >= complen*3)
351 gregl 3.1 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 gregl 3.2 /* drop satisfied requests */
363 gregl 3.11 for (i = complen; i-- && complist[i].nr <= complist[i].nc; )
364 gregl 3.2 ;
365 gregl 3.4 if (i < 0) {
366     free((char *)complist);
367     complist = NULL;
368     complen = 0;
369     } else if (i < complen-1) {
370 gregl 3.2 list2 = (PACKHEAD *)realloc((char *)complist,
371     (i+1)*sizeof(PACKHEAD));
372 gwlarson 3.26 if (list2 != NULL)
373 gregl 3.2 complist = list2;
374 gwlarson 3.26 complen = i+1;
375 gregl 3.2 }
376     listpos = 0; lastin = i;
377 gregl 3.1 }
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 gwlarson 3.27 * up, and if we get further in our compute list than this, we re-sort the
386 gregl 3.11 * list and start again from the beginning. Since
387     * a merge sort is used, the sorting costs are minimal.
388 gregl 3.1 */
389 gregl 3.17 next_packet(p, n) /* prepare packet for computation */
390 gregl 3.1 register PACKET *p;
391 gregl 3.17 int n;
392 gregl 3.1 {
393     register int i;
394    
395 gregl 3.10 if (listpos > lastin) /* time to sort the list */
396     sortcomplist();
397 gregl 3.1 if (complen <= 0)
398     return(0);
399     p->hd = complist[listpos].hd;
400     p->bi = complist[listpos].bi;
401 gregl 3.11 p->nc = complist[listpos].nc;
402     p->nr = complist[listpos].nr - p->nc;
403 gregl 3.1 if (p->nr <= 0)
404     return(0);
405 gwlarson 3.28 DCHECK(n < 1 | n > RPACKSIZ,
406     CONSISTENCY, "next_packet called with bad n value");
407 gregl 3.17 if (p->nr > n)
408     p->nr = n;
409 gregl 3.11 complist[listpos].nc += p->nr; /* find where this one would go */
410 gwlarson 3.27 if (hdgetbeam(hdlist[p->hd], p->bi) != NULL)
411     hdfreefrag(hdlist[p->hd], p->bi);
412 gregl 3.11 while (lastin > listpos &&
413     beamcmp(complist+lastin, complist+listpos) > 0)
414 gregl 3.1 lastin--;
415     listpos++;
416     return(1);
417     }