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