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