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