src/rt/oocsort.c

#ifndef lint
static const char RCSid[] = "$Id$";
#endif


/* 
   =========================================================================
   N-way out-of-core merge sort for records with 3D keys.  Recursively
   subdivides input into N blocks until these are of sufficiently small size
   to be sorted in-core according to Z-order (Morton code), then N-way
   merging them out-of-core using a priority queue as the stack unwinds.
   
   Roland Schregle (roland.schregle@{hslu.ch, gmail.com}) 
   (c) Lucerne University of Applied Sciences and Arts, 
       supported by the Swiss National Science Foundation (SNSF, #147053)
   ==========================================================================
   
   $Id: oocsort.c,v 2.3 2016/05/17 17:39:47 rschregle Exp $
*/


#if !defined(_WIN32) && !defined(_WIN64) || defined(PMAP_OOC)
/* No Windoze support for now */

#include "oocsort.h"
#include "oocmorton.h"
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <fcntl.h>
#include <sys/wait.h>


/* Priority queue node */
typedef struct {
   OOC_MortonIdx     pri;        /* Record's priority (sort key) */
   unsigned          blk;        /* Block containing record */
} OOC_SortQueueNode;


/* Priority queue */
typedef struct {
   OOC_SortQueueNode *node;
   unsigned          len, tail;
} OOC_SortQueue;


/* Additional data for qsort() compare function. We resort to instancing
 * this as a global variable instead of passing it to the compare func via
 * qsort_r(), since the latter is a non-portable GNU extension. */
typedef struct {
   RREAL    *(*key)(const void*);   /* Callback to access 3D key in record */
   FVECT    bbOrg;                  /* Origin of bbox containing keys */
   RREAL    mortonScale;            /* Scale for generating Morton code */
} OOC_KeyData;

static OOC_KeyData keyData;


static int OOC_KeyCompare (const void *rec1, const void *rec2)
/* Comparison function for in-core quicksort */
{
   OOC_MortonIdx  pri1, pri2;
   
   if (!rec1 || !rec2)
      return 0;
      
   pri1 = OOC_Key2Morton(keyData.key(rec1), keyData.bbOrg,
                         keyData.mortonScale);
   pri2 = OOC_Key2Morton(keyData.key(rec2), keyData.bbOrg, 
                         keyData.mortonScale);
   
   if (pri1 < pri2)
      return -1;
   else if (pri1 > pri2)
      return 1;
   else 
      return 0;
}


static int OOC_SortRead (FILE *file, unsigned recSize, char *rec)
/* Read next record from file; return 0 and record in rec on success, 
 * else -1 */
{
   if (!rec || feof(file) || !fread(rec, recSize, 1, file))
      return -1;

   return 0;
}


static int OOC_SortPeek (FILE *file, unsigned recSize, char *rec)
/* Return next record from file WITHOUT advancing file position;
 * return 0 and record in rec on success, else -1 */
{
   const unsigned long filePos = ftell(file);

   if (OOC_SortRead(file, recSize, rec))
      return -1;
      
   /* Workaround; fseek(file, -recSize, SEEK_CUR) causes subsequent
    * fread()'s to fail until rewind() */
   rewind(file);    
   if (fseek(file, filePos, SEEK_SET) < 0)
      return -1;

   return 0;
}


static int OOC_SortWrite (FILE *file, unsigned recSize, const char *rec)
/* Output record to file; return 0 on success, else -1 */
{
   if (!rec || !fwrite(rec, recSize, 1, file))
      return -1;
   
   return 0;
}


#ifdef DEBUG_OOC_SORT
   static int OOC_SortQueueCheck (OOC_SortQueue *pq, unsigned root)
   /* Priority queue sanity check */
   {
      unsigned kid;
      
      if (root < pq -> tail) {
         if ((kid = (root << 1) + 1) < pq -> tail) {
            if (pq -> node [kid].pri < pq -> node [root].pri)
               return -1;
            else return OOC_SortQueueCheck(pq, kid);
         }
            
         if ((kid = kid + 1) < pq -> tail) {
            if (pq -> node [kid].pri < pq -> node [root].pri)
               return -1;
            else return OOC_SortQueueCheck(pq, kid);
         }
      }
      
      return 0;
   }
#endif


static int OOC_SortPush (OOC_SortQueue *pq, OOC_MortonIdx pri, unsigned blk)
/* Insert specified block index into priority queue; return block index 
 * or -1 if queue is full */
{
   OOC_SortQueueNode *pqn = pq -> node;
   unsigned          kid, root;
   
   if (pq -> tail >= pq -> len)
      /* Queue full */
      return -1;
   
   /* Append node at tail and re-sort */
   kid = pq -> tail++;
   
   while (kid) {
      root = (kid - 1) >> 1;
      
      /* Compare with parent node and swap if necessary, else terminate */
      if (pri < pqn [root].pri) {
         pqn [kid].pri = pqn [root].pri;
         pqn [kid].blk = pqn [root].blk;
         kid = root;
      }
      else break;
   }
      
   pqn [kid].pri = pri;
   pqn [kid].blk = blk;   

#ifdef DEBUG_OOC_SORT
   if (OOC_SortQueueCheck(pq, 0) < 0) {
      fprintf(stderr, "OOC_Sort: priority queue inconsistency\n");
      return -1;
   }
#endif
   
   return blk;
}


static int OOC_SortPop (OOC_SortQueue *pq)
/* Remove head of priority queue and return its block index 
   or -1 if queue empty */
{
   OOC_SortQueueNode *pqn = pq -> node;
   OOC_MortonIdx     pri;
   unsigned          kid, kid2, root = 0, blk, res;
   
   if (!pq -> tail)
      /* Queue empty */
      return -1;
   
   res = pqn -> blk;
   pri = pqn [--pq -> tail].pri;
   blk = pqn [pq -> tail].blk;

   /* Replace head with tail node and re-sort */
   while ((kid = (root << 1) + 1) < pq -> tail) {
      /* Compare with smaller kid and swap if necessary, else terminate */
      if ((kid2 = (kid + 1)) < pq -> tail && pqn [kid2].pri < pqn [kid].pri)
         kid = kid2;
            
      if (pri > pqn [kid].pri) {
         pqn [root].pri = pqn [kid].pri;
         pqn [root].blk = pqn [kid].blk;
      }
      else break;
         
      root = kid;
   }
   
   pqn [root].pri = pri;
   pqn [root].blk = blk;

#ifdef DEBUG_OOC_SORT
   if (OOC_SortQueueCheck(pq, 0) < 0) {
      fprintf(stderr, "OOC_Sort: priority queue inconsistency\n");
      return -1;
   }
#endif
   
   return res;
}


/* Active subprocess counter updated by parent process; must persist when
 * recursing into OOC_SortRecurse(), hence global */
static unsigned procCnt = 0;

static int OOC_SortRecurse (FILE *in, unsigned long blkLo, 
                            unsigned long blkHi, FILE *out, 
                            unsigned numBlk, unsigned long maxBlkSize, 
                            unsigned numProc, unsigned recSize, 
                            char *sortBuf, OOC_SortQueue *pqueue, 
                            const OOC_KeyData *keyData)
/* Recursive part of OOC_Sort(). Reads block of records from input file
 * within the interval [blkLo, blkHi] and divides it into numBlk blocks
 * until the size (in bytes) does not exceed maxBlkSize, then quicksorts
 * each block into a temporary file.  These files are then mergesorted via a
 * priority queue to the output file as the stack unwinds.  NOTE: Critical
 * sections are prepended with '!!!'
 *
 * Parameters are as follows:
 * in          Input file containing unsorted records (assumed to be open)
 * blkLo       Start of current block in input file, in number of records
 * blkHi       End of current block in input file, in number of records
 * out         Output file containing sorted records (assumed to be open)
 * numBlk      Number of blocks to divide into / merge from
 * maxBlkSize  Max block size and size of in-core sort buffer, in bytes
 * numProc     Number of parallel processes for in-core sort
 * recSize     Size of input records in bytes
 * sortBuf     Preallocated in-core sort buffer of size maxBlkSize
 * pqueue      Preallocated priority queue of length numBlk for block merge
 * keyData     Aggregate data for Morton code generation and comparison
 */
{
   const unsigned long  blkLen = blkHi - blkLo + 1, 
                        blkSize = blkLen * recSize;
   int                  stat;
   
   if (!blkLen || blkHi < blkLo)
      return 0;
      
   if (blkSize > maxBlkSize) {
      unsigned long     blkLo2 = blkLo, blkHi2 = blkLo, blkSize2;
      const double      blkLen2 = (double)blkLen / numBlk;
      FILE              *blkFile [numBlk];         /* Violates C89!  */
      char              rec [recSize];             /* Ditto          */
      OOC_MortonIdx     pri;
      int               b, pid;
#ifdef DEBUG_OOC_SORT      
      unsigned long     pqCnt = 0;
#endif      
      
      /* ======================================================
       * Block too large for in-core sort -> divide into numBlk
       * subblocks and recurse
       * ====================================================== */

#ifdef DEBUG_OOC_SORT
      fprintf(stderr, "OOC_Sort: splitting block [%lu - %lu]\n", 
              blkLo, blkHi);
#endif

      for (b = 0; b < numBlk; b++) {
         /* Open temporary file as output for subblock */
         if (!(blkFile [b] = tmpfile())) {
            perror("OOC_Sort: failed opening temporary block file");
            return -1;
         }
         
         /* Subblock interval [blkLo2, blkHi2] of size blkSize2 */
         blkHi2 = blkLo + (b + 1) * blkLen2 - 1;
         blkSize2 = (blkHi2 - blkLo2 + 1) * recSize;

         if (blkSize2 <= maxBlkSize) {
            /* !!! Will be in-core sorted on recursion -> fork kid process,
             * !!! but don't fork more than numProc kids; wait if necessary */
            while (procCnt >= numProc && wait(&stat) >= 0) {
               if (!WIFEXITED(stat) || WEXITSTATUS(stat))
                  return -1;

               procCnt--;
            }

            /* Now fork kid process */
            if (!(pid = fork())) {
               /* Recurse on subblocks with new input filehandle; */
               if (OOC_SortRecurse(in, blkLo2, blkHi2, blkFile [b], numBlk,
                                   maxBlkSize, numProc, recSize, sortBuf,
                                   pqueue, keyData))
                  exit(-1);
                                                                         
               /* !!! Apparently the parent's tmpfile isn't deleted when the
                * !!! child exits (which is what we want), though some
                * !!! sources suggest using _Exit() instead; is this
                * !!! implementation specific?  */
               exit(0);
            }
            else if (pid < 0) {
               fprintf(stderr, "OOC_Sort: failed to fork subprocess\n");
               return -1;
            }

#ifdef DEBUG_OOC_FORK
            fprintf(stderr, "OOC_Sort: forking kid %d for block %d\n",
                    procCnt, b);
#endif            
      
            /* Parent continues here */
            procCnt++;
         }
         else {
            /* Recurse on subblock; without forking */
            if (OOC_SortRecurse(in, blkLo2, blkHi2, blkFile [b], numBlk,
                                maxBlkSize, numProc, recSize, sortBuf,
                                pqueue, keyData))
               return -1;
         }

         /* Prepare for next block */ 
         blkLo2 = blkHi2 + 1;
      }

      /* !!! Wait for any forked processes to terminate */
      while (procCnt && wait(&stat) >= 0) {
         if (!WIFEXITED(stat) || WEXITSTATUS(stat))
            return -1;

         procCnt--;
      }

      /* ===============================================================
       * Subblocks are now sorted; prepare priority queue by peeking and
       * enqueueing first record from corresponding temporary file
       * =============================================================== */

#ifdef DEBUG_OOC_SORT
      fprintf(stderr, "OOC_Sort: merging block [%lu - %lu]\n", blkLo, blkHi);
#endif

      for (b = 0; b < numBlk; b++) {
#ifdef DEBUG_OOC_SORT
         fseek(blkFile [b], 0, SEEK_END);
         if (ftell(blkFile [b]) % recSize) {
            fprintf(stderr, "OOC_Sort: truncated record in tmp block "
                    "file %d\n", b);
            return -1;
         }
          
         fprintf(stderr, "OOC_Sort: tmp block file %d contains %ld rec\n",
                 b, ftell(blkFile [b]) / recSize);
#endif

         rewind(blkFile [b]);

         if (OOC_SortPeek(blkFile [b], recSize, rec)) {
            fprintf(stderr, "OOC_Sort: error reading from block file\n");
            return -1;
         }      
         
         /* Enqueue record along with its Morton code as priority */
         pri = OOC_Key2Morton(keyData -> key(rec), keyData -> bbOrg, 
                              keyData -> mortonScale);
                              
         if (OOC_SortPush(pqueue, pri, b) < 0) {
            fprintf(stderr, "OOC_Sort: failed priority queue insertion\n");
            return -1;
         }
      }
      
      /* ==========================================================
       * Subblocks now sorted and priority queue filled; merge from
       * temporary files 
       * ========================================================== */

      do {
         /* Get next node in priority queue, read next record in corresponding
          * block, and send to output */
         b = OOC_SortPop(pqueue);
               
         if (b >= 0) {
            /* Priority queue non-empty */
            if (OOC_SortRead(blkFile [b], recSize, rec)) {
               /* Corresponding record should still be in the input block */
               fprintf(stderr, "OOC_Sort: unexpected end reading block file\n");
               return -1;
            }

            if (OOC_SortWrite(out, recSize, rec)) {
               fprintf(stderr, "OOC_Sort; error writing output file\n");
               return -1;
            }

#ifdef DEBUG_OOC_SORT
            pqCnt++;
#endif            
            
            /* Peek next record from same block and insert in priority queue */
            if (!OOC_SortPeek(blkFile [b], recSize, rec)) {
               /* Block not exhausted */
               pri = OOC_Key2Morton(keyData -> key(rec), keyData -> bbOrg, 
                                    keyData -> mortonScale);
               
               if (OOC_SortPush(pqueue, pri, b) < 0) {
                  fprintf(stderr, "OOC_Sort: failed priority queue insert\n");
                  return -1;
               }
            }
         }
      } while (b >= 0);
      
#ifdef DEBUG_OOC_SORT
      fprintf(stderr, "OOC_Sort: dequeued %lu rec\n", pqCnt);
      fprintf(stderr, "OOC_Sort: merged file contains %lu rec\n",
              ftell(out) / recSize);
#endif
      
      /* Priority queue now empty -> done; close temporary subblock files,
       * causing them to be automatically deleted. */
      for (b = 0; b < numBlk; b++) {
         fclose(blkFile [b]);
      }
         
      /* !!! Commit output file to disk before caller reads it; omitting
       * !!! this apparently leads to corrupted files (race condition?) when
       * !!! the caller tries to read them! */
      fflush(out);
      fsync(fileno(out));
   }

   else {   
      /* ======================================
       * Block is small enough for in-core sort
       * ====================================== */   
      int   ifd = fileno(in), ofd = fileno(out);

#ifdef DEBUG_OOC_SORT
      fprintf(stderr, "OOC_Sort: Proc %d (%d/%d) sorting block [%lu - %lu]\n", 
              getpid(), procCnt, numProc - 1, blkLo, blkHi);

#endif

      /* Atomically seek and read block into in-core sort buffer */
      /* !!! Unbuffered I/O via pread() avoids potential race conditions
       * !!! and buffer corruption which can occur with lseek()/fseek()
       * !!! followed by read()/fread(). */
      if (pread(ifd, sortBuf, blkSize, blkLo * recSize) != blkSize) {
         perror("OOC_Sort: error reading from input file");
         return -1;
      }            

      /* Quicksort block in-core and write to output file */
      qsort(sortBuf, blkLen, recSize, OOC_KeyCompare);

      if (write(ofd, sortBuf, blkSize) != blkSize) {
         perror("OOC_Sort: error writing to block file");
         return -1;
      }
      
      fsync(ofd);

#ifdef DEBUG_OOC_SORT
      fprintf(stderr, "OOC_Sort: proc %d wrote %ld records\n",
              getpid(), lseek(ofd, 0, SEEK_END) / recSize);
#endif      
   }      

   return 0;
}  

 
int OOC_Sort (FILE *in, FILE *out, unsigned numBlk, 
              unsigned long blkSize, unsigned numProc, unsigned recSize, 
              FVECT bbOrg, RREAL bbSize, RREAL *(*key)(const void*))
/* Sort records in inFile and append to outFile by subdividing inFile into
 * small blocks, sorting these in-core, and merging them out-of-core via a
 * priority queue.
 *
 * This is implemented as a recursive (numBlk)-way sort; the input is
 * successively split into numBlk smaller blocks until these are of size <=
 * blkSize, i.e. small enough for in-core sorting, then merging the sorted
 * blocks as the stack unwinds. The in-core sort is parallelised over
 * numProx processes. 
 *
 * Parameters are as follows:
 * inFile      Opened input file containing unsorted records
 * outFile     Opened output file containing sorted records
 * numBlk      Number of blocks to divide into / merge from 
 * blkSize     Max block size and size of in-core sort buffer, in bytes 
 * numProc     Number of parallel processes for in-core sort 
 * recSize     Size of input records in bytes
 * bbOrg       Origin of bounding box containing record keys for Morton code
 * bbSize      Extent of bounding box containing record keys for Morton code
 * key         Callback to access 3D coords from records for Morton code
 */
{
   unsigned long     numRec;
   OOC_SortQueue     pqueue;
   char              *sortBuf = NULL;
   int               stat;
   
   if (numBlk < 1)
      numBlk = 1;

   /* Open input file & get size in number of records */
   if (fseek(in, 0, SEEK_END) < 0) {
      fputs("OOC_Sort: failed seek in input file\n", stderr);
      return -1;
   }
   
   if (!(numRec = ftell(in) / recSize)) {
      fputs("OOC_Sort: empty input file\n", stderr);
      return -1;
   }

   /* Allocate & init priority queue */
   if (!(pqueue.node = calloc(numBlk, sizeof(OOC_SortQueueNode)))) {
      fputs("OOC_Sort: failure allocating priority queue\n", stderr);
      return -1;
   }
   pqueue.tail = 0;
   pqueue.len = numBlk;

   /* Allocate in-core sort buffa */
   if (!(sortBuf = malloc(blkSize))) {
      fprintf(stderr, "OOC_Sort: failure allocating in-core sort buffer");
      return -1;
   }
   
   /* Set up key data to pass to qsort()'s comparison func */
   keyData.key = key;
   keyData.mortonScale = OOC_MORTON_MAX / bbSize;
   VCOPY(keyData.bbOrg, bbOrg);
     
   stat = OOC_SortRecurse(in, 0, numRec - 1, out, numBlk, blkSize, numProc,
                          recSize, sortBuf, &pqueue, &keyData);

   /* Cleanup */
   free(pqueue.node);
   free(sortBuf);
   
   return stat;           
}

#endif /* NIX / PMAP_OOC */
Revision:	2.4
Committed:	Mon Aug 14 21:12:10 2017 UTC (6 years, 9 months ago) by rschregle
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad5R2, rad5R3, rad5R1
Changes since 2.3:	+10 -2 lines
Log Message:	Updated photon map code for Windows; no multproc or ooC for now
#	User	Rev	Content
1	rschregle	2.4	#ifndef lint
2			static const char RCSid[] = "$Id$";
3			#endif
4
5
6	greg	2.1	/*
7	rschregle	2.3	=========================================================================
8			N-way out-of-core merge sort for records with 3D keys. Recursively
9			subdivides input into N blocks until these are of sufficiently small size
10			to be sorted in-core according to Z-order (Morton code), then N-way
11			merging them out-of-core using a priority queue as the stack unwinds.
12
13			Roland Schregle (roland.schregle@{hslu.ch, gmail.com})
14			(c) Lucerne University of Applied Sciences and Arts,
15			supported by the Swiss National Science Foundation (SNSF, #147053)
16			==========================================================================
17	greg	2.1
18	rschregle	2.4	$Id: oocsort.c,v 2.3 2016/05/17 17:39:47 rschregle Exp $
19	greg	2.1	*/
20
21
22	rschregle	2.4	#if !defined(_WIN32) && !defined(_WIN64) \|\| defined(PMAP_OOC)
23			/* No Windoze support for now */
24	rschregle	2.3
25	greg	2.1	#include "oocsort.h"
26	rschregle	2.3	#include "oocmorton.h"
27	greg	2.1	#include <stdio.h>
28			#include <stdlib.h>
29	rschregle	2.3	#include <unistd.h>
30			#include <fcntl.h>
31			#include <sys/wait.h>
32
33	greg	2.1
34			/* Priority queue node */
35			typedef struct {
36	rschregle	2.3	OOC_MortonIdx pri; /* Record's priority (sort key) */
37			unsigned blk; /* Block containing record */
38			} OOC_SortQueueNode;
39
40	greg	2.1
41	rschregle	2.3	/* Priority queue */
42	greg	2.1	typedef struct {
43	rschregle	2.3	OOC_SortQueueNode *node;
44			unsigned len, tail;
45			} OOC_SortQueue;
46	greg	2.1
47	rschregle	2.3
48			/* Additional data for qsort() compare function. We resort to instancing
49			* this as a global variable instead of passing it to the compare func via
50			* qsort_r(), since the latter is a non-portable GNU extension. */
51	greg	2.1	typedef struct {
52	rschregle	2.3	RREAL (key)(const void); / Callback to access 3D key in record */
53			FVECT bbOrg; /* Origin of bbox containing keys */
54			RREAL mortonScale; /* Scale for generating Morton code */
55			} OOC_KeyData;
56	greg	2.1
57	rschregle	2.3	static OOC_KeyData keyData;
58	greg	2.1
59
60
61	rschregle	2.3	static int OOC_KeyCompare (const void rec1, const void rec2)
62			/* Comparison function for in-core quicksort */
63	greg	2.1	{
64	rschregle	2.3	OOC_MortonIdx pri1, pri2;
65	greg	2.1
66			if (!rec1 \|\| !rec2)
67			return 0;
68	rschregle	2.3
69			pri1 = OOC_Key2Morton(keyData.key(rec1), keyData.bbOrg,
70			keyData.mortonScale);
71			pri2 = OOC_Key2Morton(keyData.key(rec2), keyData.bbOrg,
72			keyData.mortonScale);
73	greg	2.1
74			if (pri1 < pri2)
75			return -1;
76			else if (pri1 > pri2)
77			return 1;
78			else
79			return 0;
80			}
81
82
83	rschregle	2.3
84			static int OOC_SortRead (FILE file, unsigned recSize, char rec)
85			/* Read next record from file; return 0 and record in rec on success,
86			* else -1 */
87	greg	2.1	{
88	rschregle	2.3	if (!rec \|\| feof(file) \|\| !fread(rec, recSize, 1, file))
89			return -1;
90
91			return 0;
92	greg	2.1	}
93
94
95	rschregle	2.3
96			static int OOC_SortPeek (FILE file, unsigned recSize, char rec)
97			/* Return next record from file WITHOUT advancing file position;
98			* return 0 and record in rec on success, else -1 */
99	greg	2.1	{
100	rschregle	2.3	const unsigned long filePos = ftell(file);
101
102			if (OOC_SortRead(file, recSize, rec))
103			return -1;
104
105			/* Workaround; fseek(file, -recSize, SEEK_CUR) causes subsequent
106			* fread()'s to fail until rewind() */
107			rewind(file);
108			if (fseek(file, filePos, SEEK_SET) < 0)
109			return -1;
110
111			return 0;
112	greg	2.1	}
113
114
115	rschregle	2.3
116			static int OOC_SortWrite (FILE file, unsigned recSize, const char rec)
117			/* Output record to file; return 0 on success, else -1 */
118	greg	2.1	{
119	rschregle	2.3	if (!rec \|\| !fwrite(rec, recSize, 1, file))
120			return -1;
121	greg	2.1
122	rschregle	2.3	return 0;
123	greg	2.1	}
124
125
126	rschregle	2.3
127			#ifdef DEBUG_OOC_SORT
128			static int OOC_SortQueueCheck (OOC_SortQueue *pq, unsigned root)
129	greg	2.1	/* Priority queue sanity check */
130			{
131			unsigned kid;
132
133			if (root < pq -> tail) {
134	rschregle	2.3	if ((kid = (root << 1) + 1) < pq -> tail) {
135	greg	2.1	if (pq -> node [kid].pri < pq -> node [root].pri)
136			return -1;
137	rschregle	2.3	else return OOC_SortQueueCheck(pq, kid);
138			}
139	greg	2.1
140	rschregle	2.3	if ((kid = kid + 1) < pq -> tail) {
141	greg	2.1	if (pq -> node [kid].pri < pq -> node [root].pri)
142			return -1;
143	rschregle	2.3	else return OOC_SortQueueCheck(pq, kid);
144			}
145	greg	2.1	}
146
147			return 0;
148			}
149			#endif
150
151
152	rschregle	2.3
153			static int OOC_SortPush (OOC_SortQueue *pq, OOC_MortonIdx pri, unsigned blk)
154	greg	2.1	/* Insert specified block index into priority queue; return block index
155			* or -1 if queue is full */
156			{
157	rschregle	2.3	OOC_SortQueueNode *pqn = pq -> node;
158			unsigned kid, root;
159	greg	2.1
160			if (pq -> tail >= pq -> len)
161			/* Queue full */
162			return -1;
163
164			/* Append node at tail and re-sort */
165			kid = pq -> tail++;
166
167			while (kid) {
168	rschregle	2.3	root = (kid - 1) >> 1;
169	greg	2.1
170			/* Compare with parent node and swap if necessary, else terminate */
171			if (pri < pqn [root].pri) {
172			pqn [kid].pri = pqn [root].pri;
173			pqn [kid].blk = pqn [root].blk;
174			kid = root;
175			}
176			else break;
177			}
178
179			pqn [kid].pri = pri;
180			pqn [kid].blk = blk;
181
182	rschregle	2.3	#ifdef DEBUG_OOC_SORT
183			if (OOC_SortQueueCheck(pq, 0) < 0) {
184	greg	2.1	fprintf(stderr, "OOC_Sort: priority queue inconsistency\n");
185			return -1;
186			}
187			#endif
188
189			return blk;
190			}
191
192
193	rschregle	2.3
194			static int OOC_SortPop (OOC_SortQueue *pq)
195	greg	2.1	/* Remove head of priority queue and return its block index
196	rschregle	2.3	or -1 if queue empty */
197	greg	2.1	{
198	rschregle	2.3	OOC_SortQueueNode *pqn = pq -> node;
199			OOC_MortonIdx pri;
200			unsigned kid, kid2, root = 0, blk, res;
201	greg	2.1
202			if (!pq -> tail)
203			/* Queue empty */
204			return -1;
205
206			res = pqn -> blk;
207			pri = pqn [--pq -> tail].pri;
208			blk = pqn [pq -> tail].blk;
209
210			/* Replace head with tail node and re-sort */
211			while ((kid = (root << 1) + 1) < pq -> tail) {
212	rschregle	2.3	/* Compare with smaller kid and swap if necessary, else terminate */
213	greg	2.1	if ((kid2 = (kid + 1)) < pq -> tail && pqn [kid2].pri < pqn [kid].pri)
214			kid = kid2;
215
216			if (pri > pqn [kid].pri) {
217			pqn [root].pri = pqn [kid].pri;
218			pqn [root].blk = pqn [kid].blk;
219			}
220			else break;
221
222			root = kid;
223			}
224
225			pqn [root].pri = pri;
226			pqn [root].blk = blk;
227
228	rschregle	2.3	#ifdef DEBUG_OOC_SORT
229			if (OOC_SortQueueCheck(pq, 0) < 0) {
230	greg	2.1	fprintf(stderr, "OOC_Sort: priority queue inconsistency\n");
231			return -1;
232			}
233			#endif
234
235			return res;
236			}
237
238	rschregle	2.3
239
240			/* Active subprocess counter updated by parent process; must persist when
241			* recursing into OOC_SortRecurse(), hence global */
242			static unsigned procCnt = 0;
243
244			static int OOC_SortRecurse (FILE *in, unsigned long blkLo,
245			unsigned long blkHi, FILE *out,
246			unsigned numBlk, unsigned long maxBlkSize,
247			unsigned numProc, unsigned recSize,
248			char sortBuf, OOC_SortQueue pqueue,
249			const OOC_KeyData *keyData)
250			/* Recursive part of OOC_Sort(). Reads block of records from input file
251			* within the interval [blkLo, blkHi] and divides it into numBlk blocks
252			* until the size (in bytes) does not exceed maxBlkSize, then quicksorts
253			* each block into a temporary file. These files are then mergesorted via a
254			* priority queue to the output file as the stack unwinds. NOTE: Critical
255			* sections are prepended with '!!!'
256			*
257			* Parameters are as follows:
258			* in Input file containing unsorted records (assumed to be open)
259			* blkLo Start of current block in input file, in number of records
260			* blkHi End of current block in input file, in number of records
261			* out Output file containing sorted records (assumed to be open)
262			* numBlk Number of blocks to divide into / merge from
263			* maxBlkSize Max block size and size of in-core sort buffer, in bytes
264			* numProc Number of parallel processes for in-core sort
265			* recSize Size of input records in bytes
266			* sortBuf Preallocated in-core sort buffer of size maxBlkSize
267			* pqueue Preallocated priority queue of length numBlk for block merge
268			* keyData Aggregate data for Morton code generation and comparison
269			*/
270	greg	2.1	{
271	rschregle	2.3	const unsigned long blkLen = blkHi - blkLo + 1,
272			blkSize = blkLen * recSize;
273			int stat;
274	greg	2.1
275	rschregle	2.3	if (!blkLen \|\| blkHi < blkLo)
276			return 0;
277
278			if (blkSize > maxBlkSize) {
279			unsigned long blkLo2 = blkLo, blkHi2 = blkLo, blkSize2;
280			const double blkLen2 = (double)blkLen / numBlk;
281			FILE blkFile [numBlk]; / Violates C89! */
282			char rec [recSize]; /* Ditto */
283			OOC_MortonIdx pri;
284			int b, pid;
285			#ifdef DEBUG_OOC_SORT
286			unsigned long pqCnt = 0;
287			#endif
288
289			/* ======================================================
290			* Block too large for in-core sort -> divide into numBlk
291			* subblocks and recurse
292			* ====================================================== */
293
294			#ifdef DEBUG_OOC_SORT
295			fprintf(stderr, "OOC_Sort: splitting block [%lu - %lu]\n",
296			blkLo, blkHi);
297			#endif
298
299			for (b = 0; b < numBlk; b++) {
300			/* Open temporary file as output for subblock */
301			if (!(blkFile [b] = tmpfile())) {
302			perror("OOC_Sort: failed opening temporary block file");
303			return -1;
304			}
305
306			/* Subblock interval [blkLo2, blkHi2] of size blkSize2 */
307			blkHi2 = blkLo + (b + 1) * blkLen2 - 1;
308			blkSize2 = (blkHi2 - blkLo2 + 1) * recSize;
309
310			if (blkSize2 <= maxBlkSize) {
311			/* !!! Will be in-core sorted on recursion -> fork kid process,
312			* !!! but don't fork more than numProc kids; wait if necessary */
313			while (procCnt >= numProc && wait(&stat) >= 0) {
314			if (!WIFEXITED(stat) \|\| WEXITSTATUS(stat))
315			return -1;
316
317			procCnt--;
318			}
319
320			/* Now fork kid process */
321			if (!(pid = fork())) {
322			/* Recurse on subblocks with new input filehandle; */
323			if (OOC_SortRecurse(in, blkLo2, blkHi2, blkFile [b], numBlk,
324			maxBlkSize, numProc, recSize, sortBuf,
325			pqueue, keyData))
326			exit(-1);
327
328			/* !!! Apparently the parent's tmpfile isn't deleted when the
329			* !!! child exits (which is what we want), though some
330			* !!! sources suggest using _Exit() instead; is this
331			* !!! implementation specific? */
332			exit(0);
333			}
334			else if (pid < 0) {
335			fprintf(stderr, "OOC_Sort: failed to fork subprocess\n");
336			return -1;
337			}
338
339			#ifdef DEBUG_OOC_FORK
340			fprintf(stderr, "OOC_Sort: forking kid %d for block %d\n",
341			procCnt, b);
342			#endif
343
344			/* Parent continues here */
345			procCnt++;
346			}
347			else {
348			/* Recurse on subblock; without forking */
349			if (OOC_SortRecurse(in, blkLo2, blkHi2, blkFile [b], numBlk,
350			maxBlkSize, numProc, recSize, sortBuf,
351			pqueue, keyData))
352			return -1;
353			}
354
355			/* Prepare for next block */
356			blkLo2 = blkHi2 + 1;
357			}
358
359			/* !!! Wait for any forked processes to terminate */
360			while (procCnt && wait(&stat) >= 0) {
361			if (!WIFEXITED(stat) \|\| WEXITSTATUS(stat))
362			return -1;
363
364			procCnt--;
365			}
366
367			/* ===============================================================
368			* Subblocks are now sorted; prepare priority queue by peeking and
369			* enqueueing first record from corresponding temporary file
370			* =============================================================== */
371
372			#ifdef DEBUG_OOC_SORT
373			fprintf(stderr, "OOC_Sort: merging block [%lu - %lu]\n", blkLo, blkHi);
374			#endif
375	greg	2.1
376	rschregle	2.3	for (b = 0; b < numBlk; b++) {
377			#ifdef DEBUG_OOC_SORT
378			fseek(blkFile [b], 0, SEEK_END);
379			if (ftell(blkFile [b]) % recSize) {
380			fprintf(stderr, "OOC_Sort: truncated record in tmp block "
381			"file %d\n", b);
382			return -1;
383			}
384
385			fprintf(stderr, "OOC_Sort: tmp block file %d contains %ld rec\n",
386			b, ftell(blkFile [b]) / recSize);
387	greg	2.1	#endif
388
389	rschregle	2.3	rewind(blkFile [b]);
390
391			if (OOC_SortPeek(blkFile [b], recSize, rec)) {
392			fprintf(stderr, "OOC_Sort: error reading from block file\n");
393			return -1;
394			}
395
396			/* Enqueue record along with its Morton code as priority */
397			pri = OOC_Key2Morton(keyData -> key(rec), keyData -> bbOrg,
398			keyData -> mortonScale);
399
400			if (OOC_SortPush(pqueue, pri, b) < 0) {
401			fprintf(stderr, "OOC_Sort: failed priority queue insertion\n");
402			return -1;
403			}
404	greg	2.1	}
405
406	rschregle	2.3	/* ==========================================================
407			* Subblocks now sorted and priority queue filled; merge from
408			* temporary files
409			* ========================================================== */
410
411			do {
412			/* Get next node in priority queue, read next record in corresponding
413			* block, and send to output */
414			b = OOC_SortPop(pqueue);
415
416			if (b >= 0) {
417			/* Priority queue non-empty */
418			if (OOC_SortRead(blkFile [b], recSize, rec)) {
419			/* Corresponding record should still be in the input block */
420			fprintf(stderr, "OOC_Sort: unexpected end reading block file\n");
421			return -1;
422			}
423
424			if (OOC_SortWrite(out, recSize, rec)) {
425			fprintf(stderr, "OOC_Sort; error writing output file\n");
426			return -1;
427			}
428
429			#ifdef DEBUG_OOC_SORT
430			pqCnt++;
431			#endif
432
433			/* Peek next record from same block and insert in priority queue */
434			if (!OOC_SortPeek(blkFile [b], recSize, rec)) {
435			/* Block not exhausted */
436			pri = OOC_Key2Morton(keyData -> key(rec), keyData -> bbOrg,
437			keyData -> mortonScale);
438
439			if (OOC_SortPush(pqueue, pri, b) < 0) {
440			fprintf(stderr, "OOC_Sort: failed priority queue insert\n");
441			return -1;
442			}
443			}
444			}
445			} while (b >= 0);
446
447			#ifdef DEBUG_OOC_SORT
448			fprintf(stderr, "OOC_Sort: dequeued %lu rec\n", pqCnt);
449			fprintf(stderr, "OOC_Sort: merged file contains %lu rec\n",
450			ftell(out) / recSize);
451			#endif
452	greg	2.1
453	rschregle	2.3	/* Priority queue now empty -> done; close temporary subblock files,
454			* causing them to be automatically deleted. */
455			for (b = 0; b < numBlk; b++) {
456			fclose(blkFile [b]);
457	greg	2.1	}
458
459	rschregle	2.3	/* !!! Commit output file to disk before caller reads it; omitting
460			* !!! this apparently leads to corrupted files (race condition?) when
461			* !!! the caller tries to read them! */
462			fflush(out);
463			fsync(fileno(out));
464			}
465
466			else {
467			/* ======================================
468			* Block is small enough for in-core sort
469			* ====================================== */
470			int ifd = fileno(in), ofd = fileno(out);
471
472			#ifdef DEBUG_OOC_SORT
473			fprintf(stderr, "OOC_Sort: Proc %d (%d/%d) sorting block [%lu - %lu]\n",
474			getpid(), procCnt, numProc - 1, blkLo, blkHi);
475
476			#endif
477
478			/* Atomically seek and read block into in-core sort buffer */
479			/* !!! Unbuffered I/O via pread() avoids potential race conditions
480			* !!! and buffer corruption which can occur with lseek()/fseek()
481			* !!! followed by read()/fread(). */
482			if (pread(ifd, sortBuf, blkSize, blkLo * recSize) != blkSize) {
483			perror("OOC_Sort: error reading from input file");
484			return -1;
485			}
486
487			/* Quicksort block in-core and write to output file */
488			qsort(sortBuf, blkLen, recSize, OOC_KeyCompare);
489
490			if (write(ofd, sortBuf, blkSize) != blkSize) {
491			perror("OOC_Sort: error writing to block file");
492	greg	2.1	return -1;
493			}
494
495	rschregle	2.3	fsync(ofd);
496
497			#ifdef DEBUG_OOC_SORT
498			fprintf(stderr, "OOC_Sort: proc %d wrote %ld records\n",
499			getpid(), lseek(ofd, 0, SEEK_END) / recSize);
500			#endif
501			}
502
503			return 0;
504			}
505
506
507	greg	2.1
508	rschregle	2.3	int OOC_Sort (FILE in, FILE out, unsigned numBlk,
509			unsigned long blkSize, unsigned numProc, unsigned recSize,
510			FVECT bbOrg, RREAL bbSize, RREAL (key)(const void*))
511			/* Sort records in inFile and append to outFile by subdividing inFile into
512			* small blocks, sorting these in-core, and merging them out-of-core via a
513			* priority queue.
514			*
515			* This is implemented as a recursive (numBlk)-way sort; the input is
516			* successively split into numBlk smaller blocks until these are of size <=
517			* blkSize, i.e. small enough for in-core sorting, then merging the sorted
518			* blocks as the stack unwinds. The in-core sort is parallelised over
519			* numProx processes.
520			*
521			* Parameters are as follows:
522			* inFile Opened input file containing unsorted records
523			* outFile Opened output file containing sorted records
524			* numBlk Number of blocks to divide into / merge from
525			* blkSize Max block size and size of in-core sort buffer, in bytes
526			* numProc Number of parallel processes for in-core sort
527			* recSize Size of input records in bytes
528			* bbOrg Origin of bounding box containing record keys for Morton code
529			* bbSize Extent of bounding box containing record keys for Morton code
530			* key Callback to access 3D coords from records for Morton code
531			*/
532			{
533			unsigned long numRec;
534			OOC_SortQueue pqueue;
535			char *sortBuf = NULL;
536			int stat;
537
538			if (numBlk < 1)
539			numBlk = 1;
540	greg	2.1
541	rschregle	2.3	/* Open input file & get size in number of records */
542			if (fseek(in, 0, SEEK_END) < 0) {
543			fputs("OOC_Sort: failed seek in input file\n", stderr);
544	greg	2.1	return -1;
545			}
546
547	rschregle	2.3	if (!(numRec = ftell(in) / recSize)) {
548			fputs("OOC_Sort: empty input file\n", stderr);
549			return -1;
550	greg	2.1	}
551	rschregle	2.3
552			/* Allocate & init priority queue */
553			if (!(pqueue.node = calloc(numBlk, sizeof(OOC_SortQueueNode)))) {
554			fputs("OOC_Sort: failure allocating priority queue\n", stderr);
555	greg	2.1	return -1;
556			}
557	rschregle	2.3	pqueue.tail = 0;
558			pqueue.len = numBlk;
559
560			/* Allocate in-core sort buffa */
561			if (!(sortBuf = malloc(blkSize))) {
562			fprintf(stderr, "OOC_Sort: failure allocating in-core sort buffer");
563	greg	2.1	return -1;
564			}
565
566	rschregle	2.3	/* Set up key data to pass to qsort()'s comparison func */
567			keyData.key = key;
568			keyData.mortonScale = OOC_MORTON_MAX / bbSize;
569			VCOPY(keyData.bbOrg, bbOrg);
570
571			stat = OOC_SortRecurse(in, 0, numRec - 1, out, numBlk, blkSize, numProc,
572			recSize, sortBuf, &pqueue, &keyData);
573	greg	2.1
574	rschregle	2.3	/* Cleanup */
575	greg	2.1	free(pqueue.node);
576	rschregle	2.3	free(sortBuf);
577	greg	2.1
578	rschregle	2.3	return stat;
579			}
580	rschregle	2.4
581			#endif /* NIX / PMAP_OOC */