1 #define CODE_lp_utils
2 
3 #include <string.h>
4 #include "commonlib.h"
5 #include "lp_lib.h"
6 #include "lp_utils.h"
7 #include <time.h>
8 #include <sys/timeb.h>
9 
10 #ifdef FORTIFY
11 # include "lp_fortify.h"
12 #endif
13 
14 #include "R.h"
15 
16 /*
17     Miscellaneous utilities as implemented for lp_solve v5.0+
18    ----------------------------------------------------------------------------------
19     Author:        Kjell Eikland
20     Contact:       kjell.eikland@broadpark.no
21     License terms: GLPL.
22 
23     Requires:      lp_utils.h, lp_lib.h
24 
25     Release notes:
26     v1.0.0  1 January 2003      Memory allocation, sorting, searching, time and
27                                 doubly linked list functions.
28     v1.1.0  15 May 2004         Added vector packing functionality
29     v1.2.0  10 January 2005     Added vector pushing/popping functionality
30                                 Modified return values and fixed problem in
31                                 linked list functions.
32 
33    ----------------------------------------------------------------------------------
34 */
35 
allocCHAR(lprec * lp,char ** ptr,int size,MYBOOL clear)36 STATIC MYBOOL allocCHAR(lprec *lp, char **ptr, int size, MYBOOL clear)
37 {
38   if(clear == TRUE)
39     *ptr = (char *) calloc(size, sizeof(**ptr));
40   else if(clear & AUTOMATIC) {
41     *ptr = (char *) realloc(*ptr, size * sizeof(**ptr));
42     if(clear & TRUE)
43       MEMCLEAR(*ptr, size);
44   }
45   else
46     *ptr = (char *) malloc(size * sizeof(**ptr));
47   if(((*ptr) == NULL) && (size > 0)) {
48     lp->report(lp, CRITICAL, "alloc of %d 'char' failed\n", size);
49     lp->spx_status = NOMEMORY;
50     return( FALSE );
51   }
52   else
53     return( TRUE );
54 }
allocMYBOOL(lprec * lp,MYBOOL ** ptr,int size,MYBOOL clear)55 STATIC MYBOOL allocMYBOOL(lprec *lp, MYBOOL **ptr, int size, MYBOOL clear)
56 {
57   if(clear == TRUE)
58     *ptr = (MYBOOL *) calloc(size, sizeof(**ptr));
59   else if(clear & AUTOMATIC) {
60     *ptr = (MYBOOL *) realloc(*ptr, size * sizeof(**ptr));
61     if(clear & TRUE)
62       MEMCLEAR(*ptr, size);
63   }
64   else
65     *ptr = (MYBOOL *) malloc(size * sizeof(**ptr));
66   if(((*ptr) == NULL) && (size > 0)) {
67     lp->report(lp, CRITICAL, "alloc of %d 'MYBOOL' failed\n", size);
68     lp->spx_status = NOMEMORY;
69     return( FALSE );
70   }
71   else
72     return( TRUE );
73 }
allocINT(lprec * lp,int ** ptr,int size,MYBOOL clear)74 STATIC MYBOOL allocINT(lprec *lp, int **ptr, int size, MYBOOL clear)
75 {
76   if(clear == TRUE)
77     *ptr = (int *) calloc(size, sizeof(**ptr));
78   else if(clear & AUTOMATIC) {
79     *ptr = (int *) realloc(*ptr, size * sizeof(**ptr));
80     if(clear & TRUE)
81       MEMCLEAR(*ptr, size);
82   }
83   else
84     *ptr = (int *) malloc(size * sizeof(**ptr));
85   if(((*ptr) == NULL) && (size > 0)) {
86     lp->report(lp, CRITICAL, "alloc of %d 'INT' failed\n", size);
87     lp->spx_status = NOMEMORY;
88     return( FALSE );
89   }
90   else
91     return( TRUE );
92 }
allocREAL(lprec * lp,REAL ** ptr,int size,MYBOOL clear)93 STATIC MYBOOL allocREAL(lprec *lp, REAL **ptr, int size, MYBOOL clear)
94 {
95   if(clear == TRUE)
96     *ptr = (REAL *) calloc(size, sizeof(**ptr));
97   else if(clear & AUTOMATIC) {
98     *ptr = (REAL *) realloc(*ptr, size * sizeof(**ptr));
99     if(clear & TRUE)
100       MEMCLEAR(*ptr, size);
101   }
102   else
103     *ptr = (REAL *) malloc(size * sizeof(**ptr));
104   if(((*ptr) == NULL) && (size > 0)) {
105     lp->report(lp, CRITICAL, "alloc of %d 'REAL' failed\n", size);
106     lp->spx_status = NOMEMORY;
107     return( FALSE );
108   }
109   else
110     return( TRUE );
111 }
allocLREAL(lprec * lp,LREAL ** ptr,int size,MYBOOL clear)112 STATIC MYBOOL allocLREAL(lprec *lp, LREAL **ptr, int size, MYBOOL clear)
113 {
114   if(clear == TRUE)
115     *ptr = (LREAL *) calloc(size, sizeof(**ptr));
116   else if(clear & AUTOMATIC) {
117     *ptr = (LREAL *) realloc(*ptr, size * sizeof(**ptr));
118     if(clear & TRUE)
119       MEMCLEAR(*ptr, size);
120   }
121   else
122     *ptr = (LREAL *) malloc(size * sizeof(**ptr));
123   if(((*ptr) == NULL) && (size > 0)) {
124     lp->report(lp, CRITICAL, "alloc of %d 'LREAL' failed\n", size);
125     lp->spx_status = NOMEMORY;
126     return( FALSE );
127   }
128   else
129     return( TRUE );
130 }
131 
allocFREE(lprec * lp,void ** ptr)132 STATIC MYBOOL allocFREE(lprec *lp, void **ptr)
133 {
134   MYBOOL status = TRUE;
135 
136   if(*ptr != NULL) {
137     free(*ptr);
138     *ptr = NULL;
139   }
140   else {
141     status = FALSE;
142     lp->report(lp, CRITICAL, "free() failed on line %d of file %s\n",
143                              __LINE__, __FILE__);
144   }
145   return(status);
146 }
147 
148 /* Do hoops to provide debugging info with FORTIFY */
149 #undef CODE_lp_utils
150 #include "lp_utils.h"
151 /* alloc-routines should always be before this line! */
152 
153 #if 0
154 !if !defined INLINE
155 void set_biton(MYBOOL *bitarray, int item)
156 {
157   bitarray[item / 8] |= (1 << (item % 8));
158 }
set_bitoff(MYBOOL * bitarray,int item)159 void set_bitoff(MYBOOL *bitarray, int item)
160 {
161   bitarray[item / 8] &= ~(1 << (item % 8));
162 }
is_biton(MYBOOL * bitarray,int item)163 MYBOOL is_biton(MYBOOL *bitarray, int item)
164 {
165   return( (MYBOOL) ((bitarray[item / 8] & (1 << (item % 8))) != 0) );
166 }
167 !endif
168 #endif
169 int comp_bits(MYBOOL *bitarray1, MYBOOL *bitarray2, int items)
170 {
171   int            i, items4, left = 0, right = 0;
172   MYBOOL         comp1;
173   unsigned long comp4;
174 
175   /* Convert items count to 8-bit representation, if necessary */
176   if(items > 0) {
177     i = items % 8;
178     items /= 8;
179     if(i)
180       items++;
181   }
182   else
183     items = -items;
184 
185   /* Do the wide unsigned integer part for speed */
186   items4 = items / sizeof(unsigned long);
187   i = 0;
188   while(i < items4) {
189     comp4 = ((unsigned long *) bitarray1)[i] &  ~((unsigned long *) bitarray2)[i];
190     if(comp4)
191       left++;
192     comp4 = ((unsigned long *) bitarray2)[i] &  ~((unsigned long *) bitarray1)[i];
193     if(comp4)
194       right++;
195     i++;
196   }
197 
198   /* Do the trailing slow narrow unsigned integer part */
199   i *= sizeof(unsigned long);
200   i++;
201   while(i < items) {
202     comp1 = bitarray1[i] & ~bitarray2[i];
203     if(comp1)
204       left++;
205     comp1 = bitarray2[i] & ~bitarray1[i];
206     if(comp1)
207       right++;
208     i++;
209   }
210 
211   /* Determine set comparison outcomes */
212   if((left > 0) && (right == 0))         /* array1 is a superset of array2 */
213     i = 1;
214   else if((left == 0) && (right > 0))   /* array2 is a superset of array1 */
215     i = -1;
216   else if((left == 0) && (right == 0))  /* array1 and array2 are identical */
217     i = 0;
218   else
219     i = -2;                              /* indicate all other outcomes */
220   return( i );
221 }
222 
223 
mempool_create(lprec * lp)224 STATIC workarraysrec *mempool_create(lprec *lp)
225 {
226   workarraysrec *temp;
227   temp = (workarraysrec *) calloc(1, sizeof(workarraysrec));
228   temp->lp = lp;
229   return( temp );
230 }
mempool_obtainVector(workarraysrec * mempool,int count,int unitsize)231 STATIC char *mempool_obtainVector(workarraysrec *mempool, int count, int unitsize)
232 {
233   char   *newmem = NULL;
234   MYBOOL *bnewmem = NULL;
235   int    *inewmem = NULL, size, i, ib, ie, memMargin = 0;
236   REAL   *rnewmem = NULL;
237 
238   /* First find the iso-sized window (binary search) */
239   size = count*unitsize;
240   memMargin += size;
241   ib = 0;
242   ie = mempool->count-1;
243   while(ie >= ib) {
244     i = (ib+ie) / 2;
245     if(abs(mempool->vectorsize[i]) > memMargin)
246       ie = i-1;
247     else if(abs(mempool->vectorsize[i]) < size)
248       ib = i+1;
249     else {
250       /* Find the beginning of the exact-sized array group */
251       do {
252         ib = i;
253         i--;
254       } while((i >= 0) && (abs(mempool->vectorsize[i]) >= size));
255       break;
256     }
257   }
258 
259   /* Check if we have a preallocated unused array of sufficient size */
260   ie = mempool->count-1;
261   for(i = ib; i <= ie; i++)
262     if(mempool->vectorsize[i] < 0)
263       break;
264 
265   /* Obtain and activate existing, unused vector if we are permitted */
266   if(i <= ie) {
267 #ifdef Paranoia
268     if((mempool->vectorsize[i] > 0) || (abs(mempool->vectorsize[i]) < size)) {
269       lprec *lp = mempool->lp;
270       lp->report(lp, SEVERE, "mempool_obtainVector: Invalid %s existing vector selected\n",
271                              (ie < 0 ? "too small" : "occupied"));
272       lp->spx_status = NOMEMORY;
273       lp->bb_break = TRUE;
274       return( newmem );
275     }
276 #endif
277     newmem = mempool->vectorarray[i];
278     mempool->vectorsize[i] *= -1;
279   }
280 
281   /* Otherwise allocate a new vector */
282   else if(unitsize == sizeof(MYBOOL)) {
283     allocMYBOOL(mempool->lp, &bnewmem, count, TRUE);
284     newmem = (char *) bnewmem;
285   }
286   else if(unitsize == sizeof(int)) {
287     allocINT(mempool->lp, &inewmem, count, TRUE);
288     newmem = (char *) inewmem;
289   }
290   else if(unitsize == sizeof(REAL)) {
291     allocREAL(mempool->lp, &rnewmem, count, TRUE);
292     newmem = (char *) rnewmem;
293   }
294 
295   /* Insert into master array if necessary (maintain sort by ascending size) */
296   if((i > ie) && (newmem != NULL)) {
297     mempool->count++;
298     if(mempool->count >= mempool->size) {
299       mempool->size += 10;
300       mempool->vectorarray = (char **) realloc(mempool->vectorarray,
301                                      sizeof(*(mempool->vectorarray))*mempool->size);
302       mempool->vectorsize  = (int *) realloc(mempool->vectorsize,
303                                      sizeof(*(mempool->vectorsize))*mempool->size);
304     }
305     ie++;
306     i = ie + 1;
307     if(i < mempool->count) {
308       MEMMOVE(mempool->vectorarray+i, mempool->vectorarray+ie, 1);
309       MEMMOVE(mempool->vectorsize+i,  mempool->vectorsize+ie,  1);
310     }
311     mempool->vectorarray[ie] = newmem;
312     mempool->vectorsize[ie]  = size;
313   }
314 
315   return( newmem );
316 }
mempool_releaseVector(workarraysrec * mempool,char * memvector,MYBOOL forcefree)317 STATIC MYBOOL mempool_releaseVector(workarraysrec *mempool, char *memvector, MYBOOL forcefree)
318 {
319   int i;
320 
321 #if 0
322   forcefree = TRUE;
323 #endif
324 
325   for(i = mempool->count-1; i >= 0; i--)
326     if(mempool->vectorarray[i] == memvector)
327       break;
328 
329   if((i < 0) || (mempool->vectorsize[i] < 0))
330     return( FALSE );
331 
332   if(forcefree) {
333     FREE(mempool->vectorarray[i]);
334     mempool->count--;
335     for(; i < mempool->count; i++)
336       mempool->vectorarray[i] = mempool->vectorarray[i+1];
337   }
338   else
339     mempool->vectorsize[i] *= -1;
340 
341   return( TRUE );
342 }
mempool_free(workarraysrec ** mempool)343 STATIC MYBOOL mempool_free(workarraysrec **mempool)
344 {
345   int i = (*mempool)->count;
346 
347   while(i > 0) {
348     i--;
349     if((*mempool)->vectorsize[i] < 0)  /* Handle unused vectors */
350       (*mempool)->vectorsize[i] *= -1;
351     mempool_releaseVector(*mempool, (*mempool)->vectorarray[i], TRUE);
352   }
353   FREE((*mempool)->vectorarray);
354   FREE((*mempool)->vectorsize);
355   FREE(*mempool);
356   return( TRUE );
357 }
358 
cloneREAL(lprec * lp,REAL * origlist,int size)359 REAL *cloneREAL(lprec *lp, REAL *origlist, int size)
360 {
361   REAL *newlist;
362 
363   size += 1;
364   if(allocREAL(lp, &newlist, size, FALSE))
365     MEMCOPY(newlist, origlist, size);
366   return(newlist);
367 }
cloneMYBOOL(lprec * lp,MYBOOL * origlist,int size)368 MYBOOL *cloneMYBOOL(lprec *lp, MYBOOL *origlist, int size)
369 {
370   MYBOOL *newlist;
371 
372   size += 1;
373   if(allocMYBOOL(lp, &newlist, size, FALSE))
374     MEMCOPY(newlist, origlist, size);
375   return(newlist);
376 }
cloneINT(lprec * lp,int * origlist,int size)377 int *cloneINT(lprec *lp, int *origlist, int size)
378 {
379   int *newlist;
380 
381   size += 1;
382   if(allocINT(lp, &newlist, size, FALSE))
383     MEMCOPY(newlist, origlist, size);
384   return(newlist);
385 }
386 
roundVector(LREAL * myvector,int endpos,LREAL roundzero)387 STATIC void roundVector(LREAL *myvector, int endpos, LREAL roundzero)
388 {
389   if(roundzero > 0)
390     for(; endpos >= 0; myvector++, endpos--)
391       if(fabs(*myvector) < roundzero)
392         *myvector = 0;
393 }
394 
normalizeVector(REAL * myvector,int endpos)395 STATIC REAL normalizeVector(REAL *myvector, int endpos)
396 /* Scale the ingoing vector so that its norm is unit, and return the original length */
397 {
398   int  i;
399   REAL SSQ;
400 
401   /* Cumulate squares */
402   SSQ = 0;
403   for(i = 0; i <= endpos; myvector++, i++)
404     SSQ += (*myvector) * (*myvector);
405 
406   /* Normalize */
407   SSQ = sqrt(SSQ);
408   if(SSQ > 0)
409     for(myvector--; i > 0; myvector--, i--)
410       (*myvector) /= SSQ;
411 
412   return( SSQ );
413 }
414 
415 /* ---------------------------------------------------------------------------------- */
416 /* Other general utilities                                                            */
417 /* ---------------------------------------------------------------------------------- */
418 
swapINT(int * item1,int * item2)419 STATIC void swapINT(int *item1, int *item2)
420 {
421   int hold = *item1;
422   *item1 = *item2;
423   *item2 = hold;
424 }
425 
swapREAL(REAL * item1,REAL * item2)426 STATIC void swapREAL(REAL *item1, REAL *item2)
427 {
428   REAL hold = *item1;
429   *item1 = *item2;
430   *item2 = hold;
431 }
432 
swapPTR(void ** item1,void ** item2)433 STATIC void swapPTR(void **item1, void **item2)
434 {
435   void *hold;
436   hold = *item1;
437   *item1 = *item2;
438   *item2 = hold;
439 }
440 
441 
restoreINT(REAL valREAL,REAL epsilon)442 STATIC REAL restoreINT(REAL valREAL, REAL epsilon)
443 {
444   REAL valINT, fracREAL, fracABS;
445 
446   fracREAL = modf(valREAL, &valINT);
447   fracABS = fabs(fracREAL);
448   if(fracABS < epsilon)
449     return(valINT);
450   else if(fracABS > 1-epsilon) {
451     if(fracREAL < 0)
452       return(valINT-1);
453     else
454       return(valINT+1);
455   }
456   return(valREAL);
457 }
458 
roundToPrecision(REAL value,REAL precision)459 STATIC REAL roundToPrecision(REAL value, REAL precision)
460 {
461 #if 1
462   REAL  vmod;
463   int   vexp2, vexp10;
464   LLONG sign;
465 
466   if(precision == 0)
467     return(value);
468 
469   sign  = my_sign(value);
470   value = fabs(value);
471 
472   /* Round to integer if possible */
473   if(value < precision)
474     return( 0 );
475   else if(value == floor(value))
476     return( value*sign );
477   else if((value < (REAL) MAXINT64) &&
478      (modf((REAL) (value+precision), &vmod) < precision)) {
479     sign *= (LLONG) (value+precision);
480     return( (REAL) sign );
481   }
482 
483   /* Optionally round with base 2 representation for additional precision */
484 #define roundPrecisionBase2
485 #ifdef roundPrecisionBase2
486   value = frexp(value, &vexp2);
487 #else
488   vexp2 = 0;
489 #endif
490 
491   /* Convert to desired precision */
492   vexp10 = (int) log10(value);
493   precision *= pow(10.0, vexp10);
494   modf(value/precision+0.5, &value);
495   value *= sign*precision;
496 
497   /* Restore base 10 representation if base 2 was active */
498   if(vexp2 != 0)
499     value = ldexp(value, vexp2);
500 #endif
501 
502   return( value );
503 }
504 
505 
506 /* ---------------------------------------------------------------------------------- */
507 /* Searching function specialized for lp_solve                                        */
508 /* ---------------------------------------------------------------------------------- */
searchFor(int target,int * attributes,int size,int offset,MYBOOL absolute)509 STATIC int searchFor(int target, int *attributes, int size, int offset, MYBOOL absolute)
510 {
511   int beginPos, endPos;
512   int newPos, match;
513 
514  /* Set starting and ending index offsets */
515   beginPos = offset;
516   endPos = beginPos + size - 1;
517 
518  /* Do binary search logic based on a sorted attribute vector */
519   newPos = (beginPos + endPos) / 2;
520   match = attributes[newPos];
521   if(absolute)
522     match = abs(match);
523   while(endPos - beginPos > LINEARSEARCH) {
524     if(match < target) {
525       beginPos = newPos + 1;
526       newPos = (beginPos + endPos) / 2;
527       match = attributes[newPos];
528       if(absolute)
529         match = abs(match);
530     }
531     else if(match > target) {
532       endPos = newPos - 1;
533       newPos = (beginPos + endPos) / 2;
534       match = attributes[newPos];
535       if(absolute)
536         match = abs(match);
537     }
538     else {
539       beginPos = newPos;
540       endPos = newPos;
541     }
542   }
543 
544  /* Do linear (unsorted) search logic */
545   if(endPos - beginPos <= LINEARSEARCH) {
546     match = attributes[beginPos];
547     if(absolute)
548       match = abs(match);
549       while((beginPos < endPos) && (match != target)) {
550         beginPos++;
551         match = attributes[beginPos];
552         if(absolute)
553           match = abs(match);
554       }
555       if(match == target)
556         endPos = beginPos;
557   }
558 
559  /* Return the index if a match was found, or signal failure with a -1 */
560   if((beginPos == endPos) && (match == target))
561     return(beginPos);
562   else
563     return(-1);
564 
565 }
566 
567 
568 /* ---------------------------------------------------------------------------------- */
569 /* Other supporting math routines                                                     */
570 /* ---------------------------------------------------------------------------------- */
571 
isINT(lprec * lp,REAL value)572 STATIC MYBOOL isINT(lprec *lp, REAL value)
573 {
574 #if 0
575   return( (MYBOOL) (modf(fabs(value)+lp->epsint, &value) < 2*lp->epsint) );
576 #elif 1
577   value = fabs(value)+lp->epsint;
578   return( (MYBOOL) (my_reldiff(value, floor(value)) < 2*lp->epsint) );
579 #elif 0
580   static REAL hold;
581   value = fabs(value);
582   hold = pow(10, MIN(-2, log10(value+1)+log10(lp->epsint)));
583   return( (MYBOOL) (modf(value+lp->epsint, &value) < 2*hold) );
584 #elif 0
585   value -= (REAL)floor(value);
586   return( (MYBOOL) ((value < lp->epsint) || (value > (1 - lp->epsint)) );
587 #else
588   value += lp->epsint;
589   return( (MYBOOL) (fabs(value-floor(value)) < 2*lp->epsint) );
590 #endif
591 }
592 
593 STATIC MYBOOL isOrigFixed(lprec *lp, int varno)
594 {
595   return( (MYBOOL) (lp->orig_upbo[varno] - lp->orig_lowbo[varno] <= lp->epsmachine) );
596 }
597 
598 STATIC void chsign_bounds(REAL *lobound, REAL *upbound)
599 {
600   REAL temp;
601   temp = *upbound;
602   if(fabs(*lobound) > 0)
603     *upbound = -(*lobound);
604   else
605     *upbound = 0;
606   if(fabs(temp) > 0)
607     *lobound = -temp;
608   else
609     *lobound = 0;
610 }
611 
612 
613 /* ---------------------------------------------------------------------------------- */
614 /* Define randomization routine                                                       */
615 /* ---------------------------------------------------------------------------------- */
616 STATIC REAL rand_uniform(lprec *lp, REAL range)
617 {
618   static MYBOOL randomized = FALSE;
619 
620   if(!randomized) {
621     GetRNGstate();
622 /* Original code:   srand((unsigned) time( NULL )); */
623     randomized = TRUE;
624   }
625 /* We need to call Put when we're done. So...every time? */
626 
627   range *= (REAL) unif_rand();
628   PutRNGstate();
629 
630 /* Original code: range *= (REAL) rand() / (REAL) RAND_MAX; */
631   return( range );
632 }
633 
634 
635 /* ---------------------------------------------------------------------------------- */
636 /* Define routines for doubly linked lists of integers                                */
637 /* ---------------------------------------------------------------------------------- */
638 
639 STATIC int createLink(int size, LLrec **linkmap, MYBOOL *usedpos)
640 {
641   int i, j;
642   MYBOOL reverse;
643 
644   *linkmap = (LLrec *) calloc(1, sizeof(**linkmap));
645   if(*linkmap == NULL)
646     return( -1 );
647 
648   reverse = (MYBOOL) (size < 0);
649   if(reverse)
650     size = -size;
651   (*linkmap)->map = (int *) calloc(2*(size + 1), sizeof(int));
652   if((*linkmap)->map == NULL)
653     return( -1 );
654 
655   (*linkmap)->size = size;
656   j = 0;
657   if(usedpos == NULL)
658     (*linkmap)->map[0] = 0;
659   else {
660     for(i = 1; i <= size; i++)
661       if(!usedpos[i] ^ reverse) {
662         /* Set the forward link */
663         (*linkmap)->map[j] = i;
664         /* Set the backward link */
665         (*linkmap)->map[size+i] = j;
666         j = i;
667         if((*linkmap)->count == 0)
668           (*linkmap)->firstitem = i;
669         (*linkmap)->lastitem = i;
670         (*linkmap)->count++;
671       }
672   }
673   (*linkmap)->map[2*size+1] = j;
674 
675   return( (*linkmap)->count );
676 }
677 
678 STATIC MYBOOL freeLink(LLrec **linkmap)
679 {
680   MYBOOL status = TRUE;
681 
682   if((linkmap == NULL) || (*linkmap == NULL))
683     status = FALSE;
684   else {
685     if((*linkmap)->map != NULL)
686       free((*linkmap)->map);
687     free(*linkmap);
688     *linkmap = NULL;
689   }
690   return( status );
691 }
692 
693 STATIC int sizeLink(LLrec *linkmap)
694 {
695   return(linkmap->size);
696 }
697 
698 STATIC MYBOOL isActiveLink(LLrec *linkmap, int itemnr)
699 {
700   if((linkmap->map[itemnr] != 0) ||
701      (linkmap->map[linkmap->size+itemnr] != 0) ||
702      (linkmap->map[0] == itemnr))
703     return( TRUE );
704   else
705     return( FALSE );
706 }
707 
708 STATIC int countActiveLink(LLrec *linkmap)
709 {
710   return(linkmap->count);
711 }
712 
713 STATIC int countInactiveLink(LLrec *linkmap)
714 {
715   return(linkmap->size-linkmap->count);
716 }
717 
718 STATIC int firstActiveLink(LLrec *linkmap)
719 {
720   return(linkmap->map[0]);
721 }
722 
723 STATIC int lastActiveLink(LLrec *linkmap)
724 {
725   return(linkmap->map[2*linkmap->size+1]);
726 }
727 
728 STATIC MYBOOL appendLink(LLrec *linkmap, int newitem)
729 {
730   int k, size;
731   size = linkmap->size;
732 
733   if(linkmap->map[newitem] != 0)
734     return( FALSE );
735 
736   /* Link forward */
737   k = linkmap->map[2*size+1];
738   linkmap->map[k] = newitem;
739 
740   /* Link backward */
741   linkmap->map[size+newitem] = k;
742   linkmap->map[2*size+1] = newitem;
743 
744   /* Update count and return */
745   if(linkmap->count == 0)
746     linkmap->firstitem = newitem;
747   linkmap->lastitem = newitem;
748   linkmap->count++;
749 
750   return( TRUE );
751 }
752 
753 STATIC MYBOOL insertLink(LLrec *linkmap, int afteritem, int newitem)
754 {
755   int k, size;
756 
757   size = linkmap->size;
758 
759   if(linkmap->map[newitem] != 0)
760     return( FALSE );
761 
762   if(afteritem == linkmap->map[2*size+1])
763     appendLink(linkmap, newitem);
764   else {
765     /* Link forward */
766     k = linkmap->map[afteritem];
767     linkmap->map[afteritem] = newitem;
768     linkmap->map[newitem] = k;
769 
770     /* Link backward */
771     linkmap->map[size+k] = newitem;
772     linkmap->map[size+newitem] = afteritem;
773 
774     /* Update count */
775     SETMIN(linkmap->firstitem, newitem);
776     SETMAX(linkmap->lastitem, newitem);
777     linkmap->count++;
778   }
779 
780   return( TRUE );
781 }
782 
783 STATIC MYBOOL setLink(LLrec *linkmap, int newitem)
784 {
785   if(isActiveLink(linkmap, newitem))
786     return( FALSE );
787   else
788     return( insertLink(linkmap, prevActiveLink(linkmap, newitem), newitem) );
789 }
790 
791 STATIC MYBOOL fillLink(LLrec *linkmap)
792 {
793   int k, size;
794   size = linkmap->size;
795 
796   k = firstActiveLink(linkmap);
797   if(k != 0)
798     return( FALSE );
799   for(k = 1; k <= size; k++)
800     appendLink(linkmap, k);
801   return( TRUE );
802 }
803 
804 STATIC int nextActiveLink(LLrec *linkmap, int backitemnr)
805 {
806   if((backitemnr < 0) || (backitemnr > linkmap->size))
807     return( -1 );
808   else {
809     if(backitemnr < linkmap->lastitem)
810     while((backitemnr > linkmap->firstitem) && (linkmap->map[backitemnr] == 0))
811       backitemnr--;
812     return(linkmap->map[backitemnr]);
813   }
814 }
815 
816 STATIC int prevActiveLink(LLrec *linkmap, int forwitemnr)
817 {
818   if((forwitemnr <= 0) || (forwitemnr > linkmap->size+1))
819     return( -1 );
820   else {
821     if(forwitemnr > linkmap->lastitem)
822       return( linkmap->lastitem);
823     if(forwitemnr > linkmap->firstitem) {
824       forwitemnr += linkmap->size;
825       while((forwitemnr < linkmap->size + linkmap->lastitem) && (linkmap->map[forwitemnr] == 0))
826         forwitemnr++;
827     }
828     else
829       forwitemnr += linkmap->size;
830     return(linkmap->map[forwitemnr]);
831   }
832 }
833 
834 STATIC int firstInactiveLink(LLrec *linkmap)
835 {
836   int i, n;
837 
838   if(countInactiveLink(linkmap) == 0)
839     return( 0 );
840   n = 1;
841   i = firstActiveLink(linkmap);
842   while(i == n) {
843     n++;
844     i = nextActiveLink(linkmap, i);
845   }
846   return( n );
847 }
848 
849 STATIC int lastInactiveLink(LLrec *linkmap)
850 {
851   int i, n;
852 
853   if(countInactiveLink(linkmap) == 0)
854     return( 0 );
855   n = linkmap->size;
856   i = lastActiveLink(linkmap);
857   while(i == n) {
858     n--;
859     i = prevActiveLink(linkmap, i);
860   }
861   return( n );
862 }
863 
864 STATIC int nextInactiveLink(LLrec *linkmap, int backitemnr)
865 {
866   do {
867     backitemnr++;
868   } while((backitemnr <= linkmap->size) && isActiveLink(linkmap, backitemnr));
869   if(backitemnr <= linkmap->size)
870     return( backitemnr );
871   else
872     return( 0 );
873 }
874 
875 STATIC int prevInactiveLink(LLrec *linkmap, int forwitemnr)
876 {
877   return( 0 );
878 }
879 
880 STATIC int removeLink(LLrec *linkmap, int itemnr)
881 {
882   int size, prevnr, nextnr = -1;
883 
884   size = linkmap->size;
885   if((itemnr <= 0) || (itemnr > size))
886     return( nextnr );
887 #ifdef Paranoia
888   if(!isActiveLink(linkmap, itemnr))
889     return( nextnr );
890 #endif
891 
892   /* Get link data at the specified position */
893   nextnr = linkmap->map[itemnr];
894   prevnr = linkmap->map[size+itemnr];
895   if(itemnr == linkmap->firstitem)
896     linkmap->firstitem = nextnr;
897   if(itemnr == linkmap->lastitem)
898     linkmap->lastitem = prevnr;
899 
900   /* Update forward link */
901   linkmap->map[prevnr] = linkmap->map[itemnr];
902   linkmap->map[itemnr] = 0;
903 
904   /* Update backward link */
905   if(nextnr == 0)
906     linkmap->map[2*size+1] = prevnr;
907   else
908     linkmap->map[size+nextnr] = linkmap->map[size+itemnr];
909   linkmap->map[size+itemnr] = 0;
910 
911   /* Decrement the count */
912   linkmap->count--;
913 
914   /* Return the next active item */
915   return( nextnr );
916 }
917 
918 STATIC LLrec *cloneLink(LLrec *sourcemap, int newsize, MYBOOL freesource)
919 {
920   LLrec *testmap = NULL;
921 
922   if((newsize == sourcemap->size) || (newsize <= 0)) {
923     createLink(sourcemap->size, &testmap, NULL);
924     MEMCOPY(testmap->map, sourcemap->map, 2*(sourcemap->size+1));
925     testmap->firstitem = sourcemap->firstitem;
926     testmap->lastitem = sourcemap->lastitem;
927     testmap->size = sourcemap->size;
928     testmap->count = sourcemap->count;
929   }
930   else {
931     int j;
932 
933     createLink(newsize, &testmap, NULL);
934     for(j = firstActiveLink(sourcemap); (j != 0) && (j <= newsize); j = nextActiveLink(sourcemap, j))
935       appendLink(testmap, j);
936   }
937   if(freesource)
938     freeLink(&sourcemap);
939 
940   return(testmap);
941 }
942 
943 STATIC int compareLink(LLrec *linkmap1, LLrec *linkmap2)
944 {
945   int test;
946 
947   test = memcmp(&linkmap1->size, &linkmap2->size, sizeof(int));
948   if(test == 0)
949     test = memcmp(&linkmap1->count, &linkmap2->count, sizeof(int));
950     if(test == 0)
951       test = memcmp(linkmap1->map, linkmap2->map, sizeof(int)*(2*linkmap1->size+1));
952 
953   return( test );
954 }
955 
956 STATIC MYBOOL verifyLink(LLrec *linkmap, int itemnr, MYBOOL doappend)
957 {
958   LLrec *testmap;
959 
960   testmap = cloneLink(linkmap, -1, FALSE);
961   if(doappend) {
962     appendLink(testmap, itemnr);
963     removeLink(testmap, itemnr);
964   }
965   else {
966     int previtem = prevActiveLink(testmap, itemnr);
967     removeLink(testmap, itemnr);
968     insertLink(testmap, previtem, itemnr);
969   }
970   itemnr = compareLink(linkmap, testmap);
971   freeLink(&testmap);
972   return((MYBOOL) (itemnr == 0));
973 }
974 
975 /* Packed vector routines */
976 STATIC PVrec *createPackedVector(int size, REAL *values, int *workvector)
977 {
978   int      i, k;
979   REGISTER REAL  ref;
980   PVrec    *newPV = NULL;
981   MYBOOL   localWV = (MYBOOL) (workvector == NULL);
982 
983   if(localWV)
984     workvector = (int *) malloc((size+1)*sizeof(*workvector));
985 
986   /* Tally equal-valued vector entries - also check if it is worth compressing */
987   k = 0;
988   workvector[k] = 1;
989   ref = values[1];
990   for(i = 2; i <= size; i++) {
991     if(fabs(ref - values[i]) > DEF_EPSMACHINE) {
992       k++;
993       workvector[k] = i;
994       ref = values[i];
995     }
996   }
997   if(k > size / 2) {
998     if(localWV)
999       FREE(workvector);
1000     return( newPV );
1001   }
1002 
1003   /* Create the packing object, adjust the position vector and allocate value vector */
1004   newPV = (PVrec *) malloc(sizeof(*newPV));
1005   k++;                            /* Adjust from index to to count */
1006   newPV->count = k;
1007   if(localWV)
1008     newPV->startpos = (int *) realloc(workvector, (k + 1)*sizeof(*(newPV->startpos)));
1009   else {
1010     newPV->startpos = (int *) malloc((k + 1)*sizeof(*(newPV->startpos)));
1011     MEMCOPY(newPV->startpos, workvector, k);
1012   }
1013   newPV->startpos[k] = size + 1;  /* Store terminal index + 1 for searching purposes */
1014   newPV->value = (REAL *) malloc(k*sizeof(*(newPV->value)));
1015 
1016   /* Fill the values vector before returning */
1017   for(i = 0; i < k; i++)
1018     newPV->value[i] = values[newPV->startpos[i]];
1019 
1020   return( newPV );
1021 }
1022 
1023 STATIC MYBOOL unpackPackedVector(PVrec *PV, REAL **target)
1024 {
1025   int      i, ii, k;
1026   REGISTER REAL ref;
1027 
1028   /* Test for validity of the target and create it if necessary */
1029   if(target == NULL)
1030     return( FALSE );
1031   if(*target == NULL)
1032     allocREAL(NULL, target, PV->startpos[PV->count], FALSE);
1033 
1034   /* Expand the packed vector into the target */
1035   i = PV->startpos[0];
1036   for(k = 0; k < PV->count; k++) {
1037     ii = PV->startpos[k+1];
1038     ref = PV->value[k];
1039     while (i < ii) {
1040       (*target)[i] = ref;
1041       i++;
1042     }
1043   }
1044   return( TRUE );
1045 }
1046 
1047 STATIC REAL getvaluePackedVector(PVrec *PV, int index)
1048 {
1049   index = searchFor(index, PV->startpos, PV->count, 0, FALSE);
1050   index = abs(index)-1;
1051   if(index >= 0)
1052     return( PV->value[index] );
1053   else
1054     return( 0 );
1055 }
1056 
1057 STATIC MYBOOL freePackedVector(PVrec **PV)
1058 {
1059   if((PV == NULL) || (*PV == NULL))
1060     return( FALSE );
1061 
1062   FREE((*PV)->value);
1063   FREE((*PV)->startpos);
1064   FREE(*PV);
1065   return( TRUE );
1066 }
1067 
1068 STATIC void pushPackedVector(PVrec *PV, PVrec *parent)
1069 {
1070   PV->parent = parent;
1071 }
1072 
1073 STATIC PVrec *popPackedVector(PVrec *PV)
1074 {
1075   PVrec *parent = PV->parent;
1076   freePackedVector(&PV);
1077   return( parent );
1078 }
1079 
1080