12345qiupeng
/
gtsam
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

remove lpsolve

release/4.3a0
thduynguyen 2014-12-16 11:28:20 -05:00
parent b39d14872a
commit 1694a56194
6 changed files with 0 additions and 2908 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

992
gtsam/3rdparty/lp_solve_5.5/shared/commonlib.c vendored
View File

@ -1,992 +0,0 @@
#include <sys/types.h>
#if defined INTEGERTIME || defined CLOCKTIME || defined PosixTime
# include <time.h>
#elif defined EnhTime
# include <windows.h>
#else
# include <sys/timeb.h>
#endif
#include <stdlib.h>
#include <stdio.h>
#ifdef WIN32
# include <io.h> /* Used in file search functions */
#endif
#include <ctype.h>
#include <string.h>
#include <float.h>
#include <math.h>
#include "commonlib.h"
#ifdef FORTIFY
# include "lp_fortify.h"
#endif
#if defined FPUexception
/* FPU exception masks */
unsigned int clearFPU()
{
return( _clearfp() );
}
unsigned int resetFPU(unsigned int mask)
{
_clearfp();
mask = _controlfp( mask, 0xfffff);
return( mask );
}
unsigned int catchFPU(unsigned int mask)
{
/* Always call _clearfp before enabling/unmasking a FPU exception */
unsigned int u = _clearfp();
/* Set the new mask by not-and'ing it with the previous settings */
u = _controlfp(0, 0);
mask = u & ~mask;
mask = _controlfp(mask, _MCW_EM);
/* Return the previous mask */
return( u );
}
#endif
/* Math operator equivalence function */
int intpow(int base, int exponent)
{
int result = 1;
while(exponent > 0) {
result *= base;
exponent--;
}
while(exponent < 0) {
result /= base;
exponent++;
}
return( result );
}
int mod(int n, int d)
{
return(n % d);
}
/* Some string functions */
void strtoup(char *s)
{
if(s != NULL)
while (*s) {
*s = toupper(*s);
s++;
}
}
void strtolo(char *s)
{
if(s != NULL)
while (*s) {
*s = tolower(*s);
s++;
}
}
void strcpyup(char *t, char *s)
{
if((s != NULL) && (t != NULL)) {
while (*s) {
*t = toupper(*s);
t++;
s++;
}
*t = '\0';
}
}
void strcpylo(char *t, char *s)
{
if((s != NULL) && (t != NULL)) {
while (*s) {
*t = tolower(*s);
t++;
s++;
}
*t = '\0';
}
}
/* Unix library naming utility function */
MYBOOL so_stdname(char *stdname, char *descname, int buflen)
{
char *ptr;
if((descname == NULL) || (stdname == NULL) || (((int) strlen(descname)) >= buflen - 6))
return( FALSE );
strcpy(stdname, descname);
if((ptr = strrchr(descname, '/')) == NULL)
ptr = descname;
else
ptr++;
stdname[(int) (ptr - descname)] = 0;
if(strncmp(ptr, "lib", 3))
strcat(stdname, "lib");
strcat(stdname, ptr);
if(strcmp(stdname + strlen(stdname) - 3, ".so"))
strcat(stdname, ".so");
return( TRUE );
}
/* Return the greatest common divisor of a and b, or -1 if it is
not defined. Return through the pointer arguments the integers
such that gcd(a,b) = c*a + b*d. */
int gcd(LLONG a, LLONG b, int *c, int *d)
{
LLONG q,r,t;
int cret,dret,C,D,rval, sgn_a = 1,sgn_b = 1, swap = 0;
if((a == 0) || (b == 0))
return( -1 );
/* Use local multiplier instances, if necessary */
if(c == NULL)
c = &cret;
if(d == NULL)
d = &dret;
/* Normalize so that 0 < a <= b */
if(a < 0){
a = -a;
sgn_a = -1;
}
if(b < 0){
b = -b;
sgn_b = -1;
}
if(b < a){
t = b;
b = a;
a = t;
swap = 1;
}
/* Now a <= b and both >= 1. */
q = b/a;
r = b - a*q;
if(r == 0) {
if(swap){
*d = 1;
*c = 0;
}
else {
*c = 1;
*d = 0;
}
*c = sgn_a*(*c);
*d = sgn_b*(*d);
return( (int) a );
}
rval = gcd(a,r,&C,&D);
if(swap){
*d = (int) (C-D*q);
*c = D;
}
else {
*d = D;
*c = (int) (C-D*q);
}
*c = sgn_a*(*c);
*d = sgn_b*(*d);
return( rval );
}
/* Array search functions */
int findIndex(int target, int *attributes, int count, int offset)
{
int focusPos, beginPos, endPos;
int focusAttrib, beginAttrib, endAttrib;
/* Set starting and ending index offsets */
beginPos = offset;
endPos = beginPos + count - 1;
if(endPos < beginPos)
return(-1);
/* Do binary search logic based on a sorted (decending) attribute vector */
focusPos = (beginPos + endPos) / 2;
beginAttrib = attributes[beginPos];
focusAttrib = attributes[focusPos];
endAttrib = attributes[endPos];
while(endPos - beginPos > LINEARSEARCH) {
if(beginAttrib == target) {
focusAttrib = beginAttrib;
endPos = beginPos;
}
else if(endAttrib == target) {
focusAttrib = endAttrib;
beginPos = endPos;
}
else if(focusAttrib < target) {
beginPos = focusPos + 1;
beginAttrib = attributes[beginPos];
focusPos = (beginPos + endPos) / 2;
focusAttrib = attributes[focusPos];
}
else if(focusAttrib > target) {
endPos = focusPos - 1;
endAttrib = attributes[endPos];
focusPos = (beginPos + endPos) / 2;
focusAttrib = attributes[focusPos];
}
else {
beginPos = focusPos;
endPos = focusPos;
}
}
/* Do linear (unsorted) search logic */
if(endPos - beginPos <= LINEARSEARCH) {
/* CPU intensive loop; provide alternative evaluation models */
#if defined DOFASTMATH
/* Do fast pointer arithmetic */
int *attptr = attributes + beginPos;
while((beginPos < endPos) && ((*attptr) < target)) {
beginPos++;
attptr++;
}
focusAttrib = (*attptr);
#else
/* Do traditional indexed access */
focusAttrib = attributes[beginPos];
while((beginPos < endPos) && (focusAttrib < target)) {
beginPos++;
focusAttrib = attributes[beginPos];
}
#endif
}
/* Return the index if a match was found, or signal failure with a -1 */
if(focusAttrib == target) /* Found; return retrieval index */
return(beginPos);
else if(focusAttrib > target) /* Not found; last item */
return(-beginPos);
else if(beginPos > offset+count-1)
return(-(endPos+1)); /* Not found; end of list */
else
return(-(beginPos+1)); /* Not found; intermediate point */
}
int findIndexEx(void *target, void *attributes, int count, int offset, int recsize, findCompare_func findCompare, MYBOOL ascending)
{
int focusPos, beginPos, endPos, compare, order;
void *focusAttrib, *beginAttrib, *endAttrib;
/* Set starting and ending index offsets */
beginPos = offset;
endPos = beginPos + count - 1;
if(endPos < beginPos)
return(-1);
order = (ascending ? -1 : 1);
/* Do binary search logic based on a sorted attribute vector */
focusPos = (beginPos + endPos) / 2;
beginAttrib = CMP_ATTRIBUTES(beginPos);
focusAttrib = CMP_ATTRIBUTES(focusPos);
endAttrib = CMP_ATTRIBUTES(endPos);
compare = 0;
while(endPos - beginPos > LINEARSEARCH) {
if(findCompare(target, beginAttrib) == 0) {
focusAttrib = beginAttrib;
endPos = beginPos;
}
else if(findCompare(target, endAttrib) == 0) {
focusAttrib = endAttrib;
beginPos = endPos;
}
else {
compare = findCompare(target, focusAttrib)*order;
if(compare < 0) {
beginPos = focusPos + 1;
beginAttrib = CMP_ATTRIBUTES(beginPos);
focusPos = (beginPos + endPos) / 2;
focusAttrib = CMP_ATTRIBUTES(focusPos);
}
else if(compare > 0) {
endPos = focusPos - 1;
endAttrib = CMP_ATTRIBUTES(endPos);
focusPos = (beginPos + endPos) / 2;
focusAttrib = CMP_ATTRIBUTES(focusPos);
}
else {
beginPos = focusPos;
endPos = focusPos;
}
}
}
/* Do linear (unsorted) search logic */
if(endPos - beginPos <= LINEARSEARCH) {
/* Do traditional indexed access */
focusAttrib = CMP_ATTRIBUTES(beginPos);
if(beginPos == endPos)
compare = findCompare(target, focusAttrib)*order;
else
while((beginPos < endPos) &&
((compare = findCompare(target, focusAttrib)*order) < 0)) {
beginPos++;
focusAttrib = CMP_ATTRIBUTES(beginPos);
}
}
/* Return the index if a match was found, or signal failure with a -1 */
if(compare == 0) /* Found; return retrieval index */
return(beginPos);
else if(compare > 0) /* Not found; last item */
return(-beginPos);
else if(beginPos > offset+count-1)
return(-(endPos+1)); /* Not found; end of list */
else
return(-(beginPos+1)); /* Not found; intermediate point */
}
/* Simple sorting and searching comparison "operators" */
int CMP_CALLMODEL compareCHAR(const void *current, const void *candidate)
{
return( CMP_COMPARE( *(char *) current, *(char *) candidate ) );
}
int CMP_CALLMODEL compareINT(const void *current, const void *candidate)
{
return( CMP_COMPARE( *(int *) current, *(int *) candidate ) );
}
int CMP_CALLMODEL compareREAL(const void *current, const void *candidate)
{
return( CMP_COMPARE( *(REAL *) current, *(REAL *) candidate ) );
}
/* Heap sort function (procedurally based on the Numerical Recipes version,
but expanded and generalized to hande any object with the use of
qsort-style comparison operator). An expanded version is also implemented,
where interchanges are reflected in a caller-initialized integer "tags" list. */
void hpsort(void *attributes, int count, int offset, int recsize, MYBOOL descending, findCompare_func findCompare)
{
register int i, j, k, ir, order;
register char *hold, *base;
char *save;
if(count < 2)
return;
offset -= 1;
attributes = CMP_ATTRIBUTES(offset);
base = CMP_ATTRIBUTES(1);
save = (char *) malloc(recsize);
if(descending)
order = -1;
else
order = 1;
k = (count >> 1) + 1;
ir = count;
for(;;) {
if(k > 1) {
MEMCOPY(save, CMP_ATTRIBUTES(--k), recsize);
}
else {
hold = CMP_ATTRIBUTES(ir);
MEMCOPY(save, hold, recsize);
MEMCOPY(hold, base, recsize);
if(--ir == 1) {
MEMCOPY(base, save, recsize);
break;
}
}
i = k;
j = k << 1;
while(j <= ir) {
hold = CMP_ATTRIBUTES(j);
if( (j < ir) && (findCompare(hold, CMP_ATTRIBUTES(j+1))*order < 0) ) {
hold += recsize;
j++;
}
if(findCompare(save, hold)*order < 0) {
MEMCOPY(CMP_ATTRIBUTES(i), hold, recsize);
i = j;
j <<= 1;
}
else
break;
}
MEMCOPY(CMP_ATTRIBUTES(i), save, recsize);
}
FREE(save);
}
void hpsortex(void *attributes, int count, int offset, int recsize, MYBOOL descending, findCompare_func findCompare, int *tags)
{
if(count < 2)
return;
if(tags == NULL) {
hpsort(attributes, count, offset, recsize, descending, findCompare);
return;
}
else {
register int i, j, k, ir, order;
register char *hold, *base;
char *save;
int savetag;
offset -= 1;
attributes = CMP_ATTRIBUTES(offset);
tags += offset;
base = CMP_ATTRIBUTES(1);
save = (char *) malloc(recsize);
if(descending)
order = -1;
else
order = 1;
k = (count >> 1) + 1;
ir = count;
for(;;) {
if(k > 1) {
MEMCOPY(save, CMP_ATTRIBUTES(--k), recsize);
savetag = tags[k];
}
else {
hold = CMP_ATTRIBUTES(ir);
MEMCOPY(save, hold, recsize);
MEMCOPY(hold, base, recsize);
savetag = tags[ir];
tags[ir] = tags[1];
if(--ir == 1) {
MEMCOPY(base, save, recsize);
tags[1] = savetag;
break;
}
}
i = k;
j = k << 1;
while(j <= ir) {
hold = CMP_ATTRIBUTES(j);
if( (j < ir) && (findCompare(hold, CMP_ATTRIBUTES(j+1))*order < 0) ) {
hold += recsize;
j++;
}
if(findCompare(save, hold)*order < 0) {
MEMCOPY(CMP_ATTRIBUTES(i), hold, recsize);
tags[i] = tags[j];
i = j;
j <<= 1;
}
else
break;
}
MEMCOPY(CMP_ATTRIBUTES(i), save, recsize);
tags[i] = savetag;
}
FREE(save);
}
}
/* This is a "specialized generic" version of C.A.R Hoare's Quick Sort algorithm.
It will handle arrays that are already sorted, and arrays with duplicate keys.
There are two versions here; one extended conventional with optional tag data
for each sortable value, and a version for the QSORTrec format. The QSORTrec
format i.a. includes the ability for to do linked list sorting. If the passed
comparison operator is NULL, the comparison is assumed to be for integers. */
#define QS_IS_switch LINEARSEARCH /* Threshold for switching to insertion sort */
void qsortex_swap(void *attributes, int l, int r, int recsize,
void *tags, int tagsize, char *save, char *savetag)
{
MEMCOPY(save, CMP_ATTRIBUTES(l), recsize);
MEMCOPY(CMP_ATTRIBUTES(l), CMP_ATTRIBUTES(r), recsize);
MEMCOPY(CMP_ATTRIBUTES(r), save, recsize);
if(tags != NULL) {
MEMCOPY(savetag, CMP_TAGS(l), tagsize);
MEMCOPY(CMP_TAGS(l), CMP_TAGS(r), tagsize);
MEMCOPY(CMP_TAGS(r), savetag, tagsize);
}
}
int qsortex_sort(void *attributes, int l, int r, int recsize, int sortorder, findCompare_func findCompare,
void *tags, int tagsize, char *save, char *savetag)
{
register int i, j, nmove = 0;
char *v;
/* Perform the a fast QuickSort */
if((r-l) > QS_IS_switch) {
i = (r+l)/2;
/* Tri-Median Method */
if(sortorder*findCompare(CMP_ATTRIBUTES(l), CMP_ATTRIBUTES(i)) > 0)
{ nmove++; qsortex_swap(attributes, l,i, recsize, tags, tagsize, save, savetag); }
if(sortorder*findCompare(CMP_ATTRIBUTES(l), CMP_ATTRIBUTES(r)) > 0)
{ nmove++; qsortex_swap(attributes, l,r, recsize, tags, tagsize, save, savetag); }
if(sortorder*findCompare(CMP_ATTRIBUTES(i), CMP_ATTRIBUTES(r)) > 0)
{ nmove++; qsortex_swap(attributes, i,r, recsize, tags, tagsize, save, savetag); }
j = r-1;
qsortex_swap(attributes, i,j, recsize, tags, tagsize, save, savetag);
i = l;
v = CMP_ATTRIBUTES(j);
for(;;) {
while(sortorder*findCompare(CMP_ATTRIBUTES(++i), v) < 0);
while(sortorder*findCompare(CMP_ATTRIBUTES(--j), v) > 0);
if(j < i) break;
nmove++; qsortex_swap(attributes, i,j, recsize, tags, tagsize, save, savetag);
}
nmove++; qsortex_swap(attributes, i,r-1, recsize, tags, tagsize, save, savetag);
nmove += qsortex_sort(attributes, l,j, recsize, sortorder, findCompare, tags, tagsize, save, savetag);
nmove += qsortex_sort(attributes, i+1,r, recsize, sortorder, findCompare, tags, tagsize, save, savetag);
}
return( nmove );
}
int qsortex_finish(void *attributes, int lo0, int hi0, int recsize, int sortorder, findCompare_func findCompare,
void *tags, int tagsize, char *save, char *savetag)
{
int i, j, nmove = 0;
/* This is actually InsertionSort, which is faster for local sorts */
for(i = lo0+1; i <= hi0; i++) {
/* Save bottom-most item */
MEMCOPY(save, CMP_ATTRIBUTES(i), recsize);
if(tags != NULL)
MEMCOPY(savetag, CMP_TAGS(i), tagsize);
/* Shift down! */
j = i;
while ((j > lo0) && (sortorder*findCompare(CMP_ATTRIBUTES(j-1), save) > 0)) {
MEMCOPY(CMP_ATTRIBUTES(j), CMP_ATTRIBUTES(j-1), recsize);
if(tags != NULL)
MEMCOPY(CMP_TAGS(j), CMP_TAGS(j-1), tagsize);
j--;
nmove++;
}
/* Store bottom-most item at the top */
MEMCOPY(CMP_ATTRIBUTES(j), save, recsize);
if(tags != NULL)
MEMCOPY(CMP_TAGS(j), savetag, tagsize);
}
return( nmove );
}
int qsortex(void *attributes, int count, int offset, int recsize, MYBOOL descending, findCompare_func findCompare, void *tags, int tagsize)
{
int iswaps = 0, sortorder = (descending ? -1 : 1);
char *save = NULL, *savetag = NULL;
/* Check and initialize to zero-based arrays */
if(count <= 1)
goto Finish;
attributes = (void *) CMP_ATTRIBUTES(offset);
save = (char *) malloc(recsize);
if((tagsize <= 0) && (tags != NULL))
tags = NULL;
else if(tags != NULL) {
tags = (void *) CMP_TAGS(offset);
savetag = (char *) malloc(tagsize);
}
count--;
/* Perform sort */
iswaps = qsortex_sort(attributes, 0, count, recsize, sortorder, findCompare, tags, tagsize, save, savetag);
#if QS_IS_switch > 0
iswaps += qsortex_finish(attributes, 0, count, recsize, sortorder, findCompare, tags, tagsize, save, savetag);
#endif
Finish:
FREE(save);
FREE(savetag);
return( iswaps );
}
#undef QS_IS_switch
/* This is a "specialized generic" version of C.A.R Hoare's Quick Sort algorithm.
It will handle arrays that are already sorted, and arrays with duplicate keys.
The implementation here requires the user to pass a comparison operator and
assumes that the array passed has the QSORTrec format, which i.a. includes
the ability for to do linked list sorting. If the passed comparison operator
is NULL, the comparison is assumed to be for integers. */
#define QS_IS_switch 4 /* Threshold for switching to insertion sort */
void QS_swap(UNIONTYPE QSORTrec a[], int i, int j)
{
UNIONTYPE QSORTrec T = a[i];
a[i] = a[j];
a[j] = T;
}
int QS_addfirst(UNIONTYPE QSORTrec a[], void *mydata)
{
a[0].pvoid2.ptr = mydata;
return( 0 );
}
int QS_append(UNIONTYPE QSORTrec a[], int ipos, void *mydata)
{
if(ipos <= 0)
ipos = QS_addfirst(a, mydata);
else
a[ipos].pvoid2.ptr = mydata;
return( ipos );
}
void QS_replace(UNIONTYPE QSORTrec a[], int ipos, void *mydata)
{
a[ipos].pvoid2.ptr = mydata;
}
void QS_insert(UNIONTYPE QSORTrec a[], int ipos, void *mydata, int epos)
{
for(; epos > ipos; epos--)
a[epos] = a[epos-1];
a[ipos].pvoid2.ptr = mydata;
}
void QS_delete(UNIONTYPE QSORTrec a[], int ipos, int epos)
{
for(; epos > ipos; epos--)
a[epos] = a[epos-1];
}
int QS_sort(UNIONTYPE QSORTrec a[], int l, int r, findCompare_func findCompare)
{
register int i, j, nmove = 0;
UNIONTYPE QSORTrec v;
/* Perform the a fast QuickSort */
if((r-l) > QS_IS_switch) {
i = (r+l)/2;
/* Tri-Median Method */
if(findCompare((char *) &a[l], (char *) &a[i]) > 0)
{ nmove++; QS_swap(a,l,i); }
if(findCompare((char *) &a[l], (char *) &a[r]) > 0)
{ nmove++; QS_swap(a,l,r); }
if(findCompare((char *) &a[i], (char *) &a[r]) > 0)
{ nmove++; QS_swap(a,i,r); }
j = r-1;
QS_swap(a,i,j);
i = l;
v = a[j];
for(;;) {
while(findCompare((char *) &a[++i], (char *) &v) < 0);
while(findCompare((char *) &a[--j], (char *) &v) > 0);
if(j < i) break;
nmove++; QS_swap (a,i,j);
}
nmove++; QS_swap(a,i,r-1);
nmove += QS_sort(a,l,j,findCompare);
nmove += QS_sort(a,i+1,r,findCompare);
}
return( nmove );
}
int QS_finish(UNIONTYPE QSORTrec a[], int lo0, int hi0, findCompare_func findCompare)
{
int i, j, nmove = 0;
UNIONTYPE QSORTrec v;
/* This is actually InsertionSort, which is faster for local sorts */
for(i = lo0+1; i <= hi0; i++) {
/* Save bottom-most item */
v = a[i];
/* Shift down! */
j = i;
while ((j > lo0) && (findCompare((char *) &a[j-1], (char *) &v) > 0)) {
a[j] = a[j-1];
j--;
nmove++;
}
/* Store bottom-most item at the top */
a[j] = v;
}
return( nmove );
}
MYBOOL QS_execute(UNIONTYPE QSORTrec a[], int count, findCompare_func findCompare, int *nswaps)
{
int iswaps = 0;
/* Check and initialize */
if(count <= 1)
goto Finish;
count--;
/* Perform sort */
iswaps = QS_sort(a, 0, count, findCompare);
#if QS_IS_switch > 0
iswaps += QS_finish(a, 0, count, findCompare);
#endif
Finish:
if(nswaps != NULL)
*nswaps = iswaps;
return( TRUE );
}
/* Simple specialized bubble/insertion sort functions */
int sortByREAL(int *item, REAL *weight, int size, int offset, MYBOOL unique)
{
int i, ii, saveI;
REAL saveW;
for(i = 1; i < size; i++) {
ii = i+offset-1;
while ((ii >= offset) && (weight[ii] >= weight[ii+1])) {
if(weight[ii] == weight[ii+1]) {
if(unique)
return(item[ii]);
}
else {
saveI = item[ii];
saveW = weight[ii];
item[ii] = item[ii+1];
weight[ii] = weight[ii+1];
item[ii+1] = saveI;
weight[ii+1] = saveW;
}
ii--;
}
}
return(0);
}
int sortByINT(int *item, int *weight, int size, int offset, MYBOOL unique)
{
int i, ii, saveI;
int saveW;
for(i = 1; i < size; i++) {
ii = i+offset-1;
while ((ii >= offset) && (weight[ii] >= weight[ii+1])) {
if(weight[ii] == weight[ii+1]) {
if(unique)
return(item[ii]);
}
else {
saveI = item[ii];
saveW = weight[ii];
item[ii] = item[ii+1];
weight[ii] = weight[ii+1];
item[ii+1] = saveI;
weight[ii+1] = saveW;
}
ii--;
}
}
return(0);
}
REAL sortREALByINT(REAL *item, int *weight, int size, int offset, MYBOOL unique)
{
int i, ii, saveW;
REAL saveI;
for(i = 1; i < size; i++) {
ii = i+offset-1;
while ((ii >= offset) && (weight[ii] >= weight[ii+1])) {
if(weight[ii] == weight[ii+1]) {
if(unique)
return(item[ii]);
}
else {
saveI = item[ii];
saveW = weight[ii];
item[ii] = item[ii+1];
weight[ii] = weight[ii+1];
item[ii+1] = saveI;
weight[ii+1] = saveW;
}
ii--;
}
}
return(0);
}
/* Time and message functions */
double timeNow(void)
{
#ifdef INTEGERTIME
return((double)time(NULL));
#elif defined CLOCKTIME
return((double)clock()/CLOCKS_PER_SEC /* CLK_TCK */);
#elif defined PosixTime
struct timespec t;
# if 0
clock_gettime(CLOCK_REALTIME, &t);
return( (double) t.tv_sec + (double) t.tv_nsec/1.0e9 );
# else
static double timeBase;
clock_gettime(CLOCK_MONOTONIC, &t);
if(timeBase == 0)
timeBase = clockNow() - ((double) t.tv_sec + (double) t.tv_nsec/1.0e9);
return( timeBase + (double) t.tv_sec + (double) t.tv_nsec/1.0e9 );
# endif
#elif defined EnhTime
static LARGE_INTEGER freq;
static double timeBase;
LARGE_INTEGER now;
QueryPerformanceCounter(&now);
if(timeBase == 0) {
QueryPerformanceFrequency(&freq);
timeBase = clockNow() - (double) now.QuadPart/(double) freq.QuadPart;
}
return( timeBase + (double) now.QuadPart/(double) freq.QuadPart );
#else
struct timeb buf;
ftime(&buf);
return((double)buf.time+((double) buf.millitm)/1000.0);
#endif
}
/* Miscellaneous reporting functions */
/* List a vector of INT values for the given index range */
void blockWriteINT(FILE *output, char *label, int *myvector, int first, int last)
{
int i, k = 0;
fprintf(output, label);
fprintf(output, "\n");
for(i = first; i <= last; i++) {
fprintf(output, " %5d", myvector[i]);
k++;
if(k % 12 == 0) {
fprintf(output, "\n");
k = 0;
}
}
if(k % 12 != 0)
fprintf(output, "\n");
}
/* List a vector of MYBOOL values for the given index range */
void blockWriteBOOL(FILE *output, char *label, MYBOOL *myvector, int first, int last, MYBOOL asRaw)
{
int i, k = 0;
fprintf(output, label);
fprintf(output, "\n");
for(i = first; i <= last; i++) {
if(asRaw)
fprintf(output, " %1d", myvector[i]);
else
fprintf(output, " %5s", my_boolstr(myvector[i]));
k++;
if(k % 36 == 0) {
fprintf(output, "\n");
k = 0;
}
}
if(k % 36 != 0)
fprintf(output, "\n");
}
/* List a vector of REAL values for the given index range */
void blockWriteREAL(FILE *output, char *label, REAL *myvector, int first, int last)
{
int i, k = 0;
fprintf(output, label);
fprintf(output, "\n");
for(i = first; i <= last; i++) {
fprintf(output, " %18g", myvector[i]);
k++;
if(k % 4 == 0) {
fprintf(output, "\n");
k = 0;
}
}
if(k % 4 != 0)
fprintf(output, "\n");
}
/* CONSOLE vector and matrix printing routines */
void printvec( int n, REAL *x, int modulo )
{
int i;
if (modulo <= 0) modulo = 5;
for (i = 1; i<=n; i++) {
if(mod(i, modulo) == 1)
printf("\n%2d:%12g", i, x[i]);
else
printf(" %2d:%12g", i, x[i]);
}
if(i % modulo != 0) printf("\n");
}
void printmatUT( int size, int n, REAL *U, int modulo)
{
int i, ll;
ll = 0;
for(i = 1; i<=n; i++) {
printvec(n-i+1, &U[ll], modulo);
ll += size-i+1;
}
}
void printmatSQ( int size, int n, REAL *X, int modulo)
{
int i, ll;
ll = 0;
for(i = 1; i<=n; i++) {
printvec(n, &X[ll], modulo);
ll += size;
}
}
/* Miscellaneous file functions */
#if defined _MSC_VER
/* Check MS versions before 7 */
#if _MSC_VER < 1300
# define intptr_t long
#endif
int fileCount( char *filemask )
{
struct _finddata_t c_file;
intptr_t hFile;
int count = 0;
/* Find first .c file in current directory */
if( (hFile = _findfirst( filemask, &c_file )) == -1L )
;
/* Iterate over all matching names */
else {
while( _findnext( hFile, &c_file ) == 0 )
count++;
_findclose( hFile );
}
return( count );
}
MYBOOL fileSearchPath( char *envvar, char *searchfile, char *foundpath )
{
char pathbuffer[_MAX_PATH];
_searchenv( searchfile, envvar, pathbuffer );
if(pathbuffer[0] == '\0')
return( FALSE );
else {
if(foundpath != NULL)
strcpy(foundpath, pathbuffer);
return( TRUE );
}
}
#endif

331
gtsam/3rdparty/lp_solve_5.5/shared/commonlib.h vendored
View File

@ -1,331 +0,0 @@
#ifndef HEADER_commonlib
#define HEADER_commonlib
#include <stdlib.h>
#include <stdio.h>
static char SpaceChars[3] = {" " "\7"};
static char NumChars[14] = {"0123456789-+."};
#define BIGNUMBER 1.0e+30
#define TINYNUMBER 1.0e-04
#define MACHINEPREC 2.22e-16
#define MATHPREC 1.0e-16
#define ERRLIMIT 1.0e-06
#ifndef LINEARSEARCH
#define LINEARSEARCH 5
#endif
#if 0
#define INTEGERTIME
#endif
/* ************************************************************************ */
/* Define loadable library function headers */
/* ************************************************************************ */
#if (defined WIN32) || (defined WIN64)
#define my_LoadLibrary(name) LoadLibrary(name)
#define my_GetProcAddress(handle, name) GetProcAddress(handle, name)
#define my_FreeLibrary(handle) FreeLibrary(handle); \
handle = NULL
#else
#define my_LoadLibrary(name) dlopen(name, RTLD_LAZY)
#define my_GetProcAddress(handle, name) dlsym(handle, name)
#define my_FreeLibrary(handle) dlclose(handle); \
handle = NULL
#endif
/* ************************************************************************ */
/* Define sizes of standard number types */
/* ************************************************************************ */
#ifndef LLONG
#if defined __BORLANDC__
#define LLONG __int64
#elif !defined _MSC_VER || _MSC_VER >= 1310
#define LLONG long long
#else
#define LLONG __int64
#endif
#endif
#ifndef MYBOOL
#if 0
#define MYBOOL unsigned int
#else
#define MYBOOL unsigned char
#endif
#endif
#ifndef REAL
#define REAL double
#endif
#ifndef BLAS_prec
#define BLAS_prec "d" /* The BLAS precision prefix must correspond to the REAL type */
#endif
#ifndef REALXP
#if 1
#define REALXP long double /* Set local accumulation variable as long double */
#else
#define REALXP REAL /* Set local accumulation as default precision */
#endif
#endif
#ifndef my_boolstr
#define my_boolstr(x) (!(x) ? "FALSE" : "TRUE")
#endif
#ifndef NULL
#define NULL 0
#endif
#ifndef FALSE
#define FALSE 0
#define TRUE 1
#endif
#ifndef DEF_STRBUFSIZE
#define DEF_STRBUFSIZE 512
#endif
#ifndef MAXINT32
#define MAXINT32 2147483647
#endif
#ifndef MAXUINT32
#define MAXUINT32 4294967295
#endif
#ifndef MAXINT64
#if defined _LONGLONG || defined __LONG_LONG_MAX__ || defined LLONG_MAX
#define MAXINT64 9223372036854775807ll
#else
#define MAXINT64 9223372036854775807l
#endif
#endif
#ifndef MAXUINT64
#if defined _LONGLONG || defined __LONG_LONG_MAX__ || defined LLONG_MAX
#define MAXUINT64 18446744073709551616ll
#else
#define MAXUINT64 18446744073709551616l
#endif
#endif
#ifndef DOFASTMATH
#define DOFASTMATH
#endif
#ifndef CALLOC
#define CALLOC(ptr, nr)\
if(!((void *) ptr = calloc((size_t)(nr), sizeof(*ptr))) && nr) {\
printf("calloc of %d bytes failed on line %d of file %s\n",\
(size_t) nr * sizeof(*ptr), __LINE__, __FILE__);\
}
#endif
#ifndef MALLOC
#define MALLOC(ptr, nr)\
if(!((void *) ptr = malloc((size_t)((size_t) (nr) * sizeof(*ptr)))) && nr) {\
printf("malloc of %d bytes failed on line %d of file %s\n",\
(size_t) nr * sizeof(*ptr), __LINE__, __FILE__);\
}
#endif
#ifndef REALLOC
#define REALLOC(ptr, nr)\
if(!((void *) ptr = realloc(ptr, (size_t)((size_t) (nr) * sizeof(*ptr)))) && nr) {\
printf("realloc of %d bytes failed on line %d of file %s\n",\
(size_t) nr * sizeof(*ptr), __LINE__, __FILE__);\
}
#endif
#ifndef FREE
#define FREE(ptr)\
if((void *) ptr != NULL) {\
free(ptr);\
ptr = NULL; \
}
#endif
#ifndef MEMCOPY
#define MEMCOPY(nptr, optr, nr)\
memcpy((nptr), (optr), (size_t)((size_t)(nr) * sizeof(*(optr))))
#endif
#ifndef MEMMOVE
#define MEMMOVE(nptr, optr, nr)\
memmove((nptr), (optr), (size_t)((size_t)(nr) * sizeof(*(optr))))
#endif
#ifndef MEMALLOCCOPY
#define MEMALLOCCOPY(nptr, optr, nr)\
{MALLOC(nptr, (size_t)(nr));\
MEMCOPY(nptr, optr, (size_t)(nr));}
#endif
#ifndef STRALLOCCOPY
#define STRALLOCCOPY(nstr, ostr)\
{nstr = (char *) malloc((size_t) (strlen(ostr) + 1));\
strcpy(nstr, ostr);}
#endif
#ifndef MEMCLEAR
/*#define useMMX*/
#ifdef useMMX
#define MEMCLEAR(ptr, nr)\
mem_set((ptr), '\0', (size_t)((size_t)(nr) * sizeof(*(ptr))))
#else
#define MEMCLEAR(ptr, nr)\
memset((ptr), '\0', (size_t)((size_t)(nr) * sizeof(*(ptr))))
#endif
#endif
#define MIN(x, y) ((x) < (y) ? (x) : (y))
#define MAX(x, y) ((x) > (y) ? (x) : (y))
#define SETMIN(x, y) if((x) > (y)) x = y
#define SETMAX(x, y) if((x) < (y)) x = y
#define LIMIT(lo, x, hi) ((x < (lo) ? lo : ((x) > hi ? hi : x)))
#define BETWEEN(x, a, b) (MYBOOL) (((x)-(a)) * ((x)-(b)) <= 0)
#define IF(t, x, y) ((t) ? (x) : (y))
#define SIGN(x) ((x) < 0 ? -1 : 1)
#define DELTA_SIZE(newSize, oldSize) ((int) ((newSize) * MIN(1.33, pow(1.5, fabs((double)newSize)/((oldSize+newSize)+1)))))
#ifndef CMP_CALLMODEL
#if (defined WIN32) || (defined WIN64)
#define CMP_CALLMODEL _cdecl
#else
#define CMP_CALLMODEL
#endif
#endif
typedef int (CMP_CALLMODEL findCompare_func)(const void *current, const void *candidate);
#define CMP_COMPARE(current, candidate) ( current < candidate ? -1 : (current > candidate ? 1 : 0) )
#define CMP_ATTRIBUTES(item) (((char *) attributes)+(item)*recsize)
#define CMP_TAGS(item) (((char *) tags)+(item)*tagsize)
#ifndef UNIONTYPE
#ifdef __cplusplus
#define UNIONTYPE
#else
#define UNIONTYPE union
#endif
#endif
/* This defines a 16 byte sort record (in both 32 and 64 bit OS-es) */
typedef struct _QSORTrec1
{
void *ptr;
void *ptr2;
} QSORTrec1;
typedef struct _QSORTrec2
{
void *ptr;
double realval;
} QSORTrec2;
typedef struct _QSORTrec3
{
void *ptr;
int intval;
int intpar1;
} QSORTrec3;
typedef struct _QSORTrec4
{
REAL realval;
int intval;
int intpar1;
} QSORTrec4;
typedef struct _QSORTrec5
{
double realval;
long int longval;
} QSORTrec5;
typedef struct _QSORTrec6
{
double realval;
double realpar1;
} QSORTrec6;
typedef struct _QSORTrec7
{
int intval;
int intpar1;
int intpar2;
int intpar3;
} QSORTrec7;
union QSORTrec
{
QSORTrec1 pvoid2;
QSORTrec2 pvoidreal;
QSORTrec3 pvoidint2;
QSORTrec4 realint2;
QSORTrec5 reallong;
QSORTrec6 real2;
QSORTrec7 int4;
};
#ifdef __cplusplus
extern "C" {
#endif
int intpow(int base, int exponent);
int mod(int n, int d);
void strtoup(char *s);
void strtolo(char *s);
void strcpyup(char *t, char *s);
void strcpylo(char *t, char *s);
MYBOOL so_stdname(char *stdname, char *descname, int buflen);
int gcd(LLONG a, LLONG b, int *c, int *d);
int findIndex(int target, int *attributes, int count, int offset);
int findIndexEx(void *target, void *attributes, int count, int offset, int recsize, findCompare_func findCompare, MYBOOL ascending);
void qsortex_swap(void *attributes, int l, int r, int recsize,
void *tags, int tagsize, char *save, char *savetag);
int qsortex(void *attributes, int count, int offset, int recsize, MYBOOL descending, findCompare_func findCompare, void *tags, int tagsize);
int CMP_CALLMODEL compareCHAR(const void *current, const void *candidate);
int CMP_CALLMODEL compareINT(const void *current, const void *candidate);
int CMP_CALLMODEL compareREAL(const void *current, const void *candidate);
void hpsort(void *attributes, int count, int offset, int recsize, MYBOOL descending, findCompare_func findCompare);
void hpsortex(void *attributes, int count, int offset, int recsize, MYBOOL descending, findCompare_func findCompare, int *tags);
void QS_swap(UNIONTYPE QSORTrec a[], int i, int j);
int QS_addfirst(UNIONTYPE QSORTrec a[], void *mydata);
int QS_append(UNIONTYPE QSORTrec a[], int ipos, void *mydata);
void QS_replace(UNIONTYPE QSORTrec a[], int ipos, void *mydata);
void QS_insert(UNIONTYPE QSORTrec a[], int ipos, void *mydata, int epos);
void QS_delete(UNIONTYPE QSORTrec a[], int ipos, int epos);
MYBOOL QS_execute(UNIONTYPE QSORTrec a[], int count, findCompare_func findCompare, int *nswaps);
int sortByREAL(int *item, REAL *weight, int size, int offset, MYBOOL unique);
int sortByINT(int *item, int *weight, int size, int offset, MYBOOL unique);
REAL sortREALByINT(REAL *item, int *weight, int size, int offset, MYBOOL unique);
double timeNow(void);
void blockWriteBOOL(FILE *output, char *label, MYBOOL *myvector, int first, int last, MYBOOL asRaw);
void blockWriteINT(FILE *output, char *label, int *myvector, int first, int last);
void blockWriteREAL(FILE *output, char *label, REAL *myvector, int first, int last);
void printvec( int n, REAL *x, int modulo );
void printmatSQ( int size, int n, REAL *X, int modulo );
void printmatUT( int size, int n, REAL *U, int modulo );
unsigned int catchFPU(unsigned int mask);
#if defined _MSC_VER
int fileCount( char *filemask );
MYBOOL fileSearchPath( char *envvar, char *searchfile, char *foundpath );
#endif
#ifdef __cplusplus
}
#endif
#endif /* HEADER_commonlib */

495
gtsam/3rdparty/lp_solve_5.5/shared/mmio.c vendored
View File

@ -1,495 +0,0 @@
/*
* Matrix Market I/O library for ANSI C
*
* See http://math.nist.gov/MatrixMarket for details.
*
* (Version 1.01, 5/2003)
*/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include "mmio.h"
int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_,
double **val_, int **I_, int **J_)
{
FILE *f;
MM_typecode matcode;
int M, N, nz;
int i;
double *val;
int *I, *J;
if ((f = fopen(fname, "r")) == NULL)
return -1;
if (mm_read_banner(f, &matcode) != 0)
{
printf("mm_read_unsymetric: Could not process Matrix Market banner ");
printf(" in file [%s]\n", fname);
return -1;
}
if ( !(mm_is_real(matcode) && mm_is_matrix(matcode) &&
mm_is_sparse(matcode)))
{
fprintf(stderr, "Sorry, this application does not support ");
fprintf(stderr, "Market Market type: [%s]\n",
mm_typecode_to_str(matcode));
return -1;
}
/* find out size of sparse matrix: M, N, nz .... */
if (mm_read_mtx_crd_size(f, &M, &N, &nz) !=0)
{
fprintf(stderr, "read_unsymmetric_sparse(): could not parse matrix size.\n");
return -1;
}
*M_ = M;
*N_ = N;
*nz_ = nz;
/* reseve memory for matrices */
I = (int *) malloc(nz * sizeof(int));
J = (int *) malloc(nz * sizeof(int));
val = (double *) malloc(nz * sizeof(double));
*val_ = val;
*I_ = I;
*J_ = J;
/* NOTE: when reading in doubles, ANSI C requires the use of the "l" */
/* specifier as in "%lg", "%lf", "%le", otherwise errors will occur */
/* (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15) */
for (i=0; i<nz; i++)
{
fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i]);
I[i]--; /* adjust from 1-based to 0-based */
J[i]--;
}
fclose(f);
return 0;
}
int mm_is_valid(MM_typecode matcode)
{
if (!mm_is_matrix(matcode)) return 0;
if (mm_is_dense(matcode) && mm_is_pattern(matcode)) return 0;
if (mm_is_real(matcode) && mm_is_hermitian(matcode)) return 0;
if (mm_is_pattern(matcode) && (mm_is_hermitian(matcode) ||
mm_is_skew(matcode))) return 0;
return 1;
}
int mm_read_banner(FILE *f, MM_typecode *matcode)
{
char line[MM_MAX_LINE_LENGTH];
char banner[MM_MAX_TOKEN_LENGTH];
char mtx[MM_MAX_TOKEN_LENGTH];
char crd[MM_MAX_TOKEN_LENGTH];
char data_type[MM_MAX_TOKEN_LENGTH];
char storage_scheme[MM_MAX_TOKEN_LENGTH];
char *p;
mm_clear_typecode(matcode);
if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL)
return MM_PREMATURE_EOF;
if (sscanf(line, "%s %s %s %s %s", banner, mtx, crd, data_type,
storage_scheme) != 5)
return MM_PREMATURE_EOF;
/* convert to lower case */
for (p=mtx; *p!='\0'; *p= (char) tolower(*p),p++);
for (p=crd; *p!='\0'; *p= (char) tolower(*p),p++);
for (p=data_type; *p!='\0'; *p= (char) tolower(*p),p++);
for (p=storage_scheme; *p!='\0'; *p= (char) tolower(*p),p++);
/* check for banner */
if (strncmp(banner, MatrixMarketBanner, strlen(MatrixMarketBanner)) != 0)
return MM_NO_HEADER;
/* first field should be "mtx" */
if (strcmp(mtx, MM_MTX_STR) != 0)
return MM_UNSUPPORTED_TYPE;
mm_set_matrix(matcode);
/* second field describes whether this is a sparse matrix (in coordinate
storgae) or a dense array */
if (strcmp(crd, MM_SPARSE_STR) == 0)
mm_set_sparse(matcode);
else
if (strcmp(crd, MM_DENSE_STR) == 0)
mm_set_dense(matcode);
else
return MM_UNSUPPORTED_TYPE;
/* third field */
if (strcmp(data_type, MM_REAL_STR) == 0)
mm_set_real(matcode);
else
if (strcmp(data_type, MM_COMPLEX_STR) == 0)
mm_set_complex(matcode);
else
if (strcmp(data_type, MM_PATTERN_STR) == 0)
mm_set_pattern(matcode);
else
if (strcmp(data_type, MM_INT_STR) == 0)
mm_set_integer(matcode);
else
return MM_UNSUPPORTED_TYPE;
/* fourth field */
if (strcmp(storage_scheme, MM_GENERAL_STR) == 0)
mm_set_general(matcode);
else
if (strcmp(storage_scheme, MM_SYMM_STR) == 0)
mm_set_symmetric(matcode);
else
if (strcmp(storage_scheme, MM_HERM_STR) == 0)
mm_set_hermitian(matcode);
else
if (strcmp(storage_scheme, MM_SKEW_STR) == 0)
mm_set_skew(matcode);
else
return MM_UNSUPPORTED_TYPE;
return 0;
}
int mm_write_mtx_crd_size(FILE *f, int M, int N, int nz)
{
if (fprintf(f, "%d %d %d\n", M, N, nz) < 0)
return MM_COULD_NOT_WRITE_FILE;
else
return 0;
}
int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz )
{
char line[MM_MAX_LINE_LENGTH];
int num_items_read;
/* set return null parameter values, in case we exit with errors */
*M = *N = *nz = 0;
/* now continue scanning until you reach the end-of-comments */
do {
if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL)
return MM_PREMATURE_EOF;
} while (line[0] == '%');
/* line[] is either blank, har M,N or M,N,nz */
if (sscanf(line, "%d %d %d", M, N, nz) >= 2)
return 0;
else
do {
num_items_read = fscanf(f, "%d %d %d", M, N, nz);
if (num_items_read == EOF) return MM_PREMATURE_EOF;
} while (num_items_read < 2);
return 0;
}
int mm_read_mtx_array_size(FILE *f, int *M, int *N)
{
char line[MM_MAX_LINE_LENGTH];
int num_items_read;
/* set return null parameter values, in case we exit with errors */
*M = *N = 0;
/* now continue scanning until you reach the end-of-comments */
do
{
if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL)
return MM_PREMATURE_EOF;
}while (line[0] == '%');
/* line[] is either blank or has M,N, nz */
if (sscanf(line, "%d %d", M, N) == 2)
return 0;
else /* we have a blank line */
do
{
num_items_read = fscanf(f, "%d %d", M, N);
if (num_items_read == EOF) return MM_PREMATURE_EOF;
}
while (num_items_read != 2);
return 0;
}
int mm_write_mtx_array_size(FILE *f, int M, int N)
{
if (fprintf(f, "%d %d\n", M, N) < 0)
return MM_COULD_NOT_WRITE_FILE;
else
return 0;
}
/*-------------------------------------------------------------------------*/
/******************************************************************/
/* use when I[], J[], and val[]J, and val[] are already allocated */
/******************************************************************/
int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int I[], int J[],
double val[], MM_typecode matcode)
{
int i;
if (mm_is_complex(matcode))
{
for (i=0; i<nz; i++)
if (fscanf(f, "%d %d %lg %lg", &I[i], &J[i], &val[2*i], &val[2*i+1])
!= 4) return MM_PREMATURE_EOF;
}
else if (mm_is_real(matcode))
{
for (i=0; i<nz; i++)
{
if (fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i])
!= 3) return MM_PREMATURE_EOF;
}
}
else if (mm_is_pattern(matcode))
{
for (i=0; i<nz; i++)
if (fscanf(f, "%d %d", &I[i], &J[i])
!= 2) return MM_PREMATURE_EOF;
}
else
return MM_UNSUPPORTED_TYPE;
return 0;
}
int mm_read_mtx_crd_entry(FILE *f, int *I, int *J,
double *real, double *imag, MM_typecode matcode)
{
if (mm_is_complex(matcode))
{
if (fscanf(f, "%d %d %lg %lg", I, J, real, imag)
!= 4) return MM_PREMATURE_EOF;
}
else if (mm_is_real(matcode))
{
if (fscanf(f, "%d %d %lg\n", I, J, real)
!= 3) return MM_PREMATURE_EOF;
}
else if (mm_is_pattern(matcode))
{
if (fscanf(f, "%d %d", I, J) != 2) return MM_PREMATURE_EOF;
}
else
return MM_UNSUPPORTED_TYPE;
return 0;
}
/************************************************************************
mm_read_mtx_crd() fills M, N, nz, array of values, and return
type code, e.g. 'MCRS'
if matrix is complex, values[] is of size 2*nz,
(nz pairs of real/imaginary values)
************************************************************************/
int mm_read_mtx_crd(char *fname, int *M, int *N, int *nz, int **I, int **J,
double **val, MM_typecode *matcode)
{
int ret_code;
FILE *f;
if (strcmp(fname, "stdin") == 0) f=stdin;
else
if ((f = fopen(fname, "r")) == NULL)
return MM_COULD_NOT_READ_FILE;
if ((ret_code = mm_read_banner(f, matcode)) != 0)
return ret_code;
if (!(mm_is_valid(*matcode) && mm_is_sparse(*matcode) &&
mm_is_matrix(*matcode)))
return MM_UNSUPPORTED_TYPE;
if ((ret_code = mm_read_mtx_crd_size(f, M, N, nz)) != 0)
return ret_code;
*I = (int *) malloc(*nz * sizeof(int));
*J = (int *) malloc(*nz * sizeof(int));
*val = NULL;
if (mm_is_complex(*matcode))
{
*val = (double *) malloc(*nz * 2 * sizeof(double));
ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val,
*matcode);
if (ret_code != 0) return ret_code;
}
else if (mm_is_real(*matcode))
{
*val = (double *) malloc(*nz * sizeof(double));
ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val,
*matcode);
if (ret_code != 0) return ret_code;
}
else if (mm_is_pattern(*matcode))
{
ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val,
*matcode);
if (ret_code != 0) return ret_code;
}
if (f != stdin) fclose(f);
return 0;
}
int mm_write_banner(FILE *f, MM_typecode matcode)
{
char *str = mm_typecode_to_str(matcode);
int ret_code;
ret_code = fprintf(f, "%s %s\n", MatrixMarketBanner, str);
/* free(str); This is a bug from the official distribution - KE fixed */
if (ret_code < 0 )
return MM_COULD_NOT_WRITE_FILE;
else
return 0;
}
int mm_write_mtx_crd(char fname[], int M, int N, int nz, int I[], int J[],
double val[], MM_typecode matcode)
{
FILE *f;
int i;
if (strcmp(fname, "stdout") == 0)
f = stdout;
else
if ((f = fopen(fname, "w")) == NULL)
return MM_COULD_NOT_WRITE_FILE;
/* print banner followed by typecode */
fprintf(f, "%s ", MatrixMarketBanner);
fprintf(f, "%s\n", mm_typecode_to_str(matcode));
/* print matrix sizes and nonzeros */
fprintf(f, "%d %d %d\n", M, N, nz);
/* print values */
if (mm_is_pattern(matcode))
for (i=0; i<nz; i++)
fprintf(f, "%d %d\n", I[i], J[i]);
else
if (mm_is_real(matcode))
for (i=0; i<nz; i++)
fprintf(f, "%d %d %20.16g\n", I[i], J[i], val[i]);
else
if (mm_is_complex(matcode))
for (i=0; i<nz; i++)
fprintf(f, "%d %d %20.16g %20.16g\n", I[i], J[i], val[2*i],
val[2*i+1]);
else
{
if (f != stdout) fclose(f);
return MM_UNSUPPORTED_TYPE;
}
if (f !=stdout) fclose(f);
return 0;
}
char *mm_typecode_to_str(MM_typecode matcode)
{
static char buffer[MM_MAX_LINE_LENGTH];
char *types[4];
/* check for MTX type */
if (mm_is_matrix(matcode))
types[0] = MM_MTX_STR;
else
return NULL;
/* check for CRD or ARR matrix */
if (mm_is_sparse(matcode))
types[1] = MM_SPARSE_STR;
else
if (mm_is_dense(matcode))
types[1] = MM_DENSE_STR;
else
return NULL;
/* check for element data type */
if (mm_is_real(matcode))
types[2] = MM_REAL_STR;
else
if (mm_is_complex(matcode))
types[2] = MM_COMPLEX_STR;
else
if (mm_is_pattern(matcode))
types[2] = MM_PATTERN_STR;
else
if (mm_is_integer(matcode))
types[2] = MM_INT_STR;
else
return NULL;
/* check for symmetry type */
if (mm_is_general(matcode))
types[3] = MM_GENERAL_STR;
else
if (mm_is_symmetric(matcode))
types[3] = MM_SYMM_STR;
else
if (mm_is_hermitian(matcode))
types[3] = MM_HERM_STR;
else
if (mm_is_skew(matcode))
types[3] = MM_SKEW_STR;
else
return NULL;
sprintf(buffer,"%s %s %s %s", types[0], types[1], types[2], types[3]);
return & buffer[0];
}

133
gtsam/3rdparty/lp_solve_5.5/shared/mmio.h vendored
View File

@ -1,133 +0,0 @@
/*
* Matrix Market I/O library for ANSI C
*
* See http://math.nist.gov/MatrixMarket for details.
*
*
*/
#ifndef MM_IO_H
#define MM_IO_H
#define MM_MAX_LINE_LENGTH 1025
#define MatrixMarketBanner "%%MatrixMarket"
#define MM_MAX_TOKEN_LENGTH 64
typedef char MM_typecode[4];
char *mm_typecode_to_str(MM_typecode matcode);
int mm_read_banner(FILE *f, MM_typecode *matcode);
int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz);
int mm_read_mtx_array_size(FILE *f, int *M, int *N);
int mm_write_banner(FILE *f, MM_typecode matcode);
int mm_write_mtx_crd_size(FILE *f, int M, int N, int nz);
int mm_write_mtx_array_size(FILE *f, int M, int N);
/********************* MM_typecode query fucntions ***************************/
#define mm_is_matrix(typecode) ((typecode)[0]=='M')
#define mm_is_sparse(typecode) ((typecode)[1]=='C')
#define mm_is_coordinate(typecode)((typecode)[1]=='C')
#define mm_is_dense(typecode) ((typecode)[1]=='A')
#define mm_is_array(typecode) ((typecode)[1]=='A')
#define mm_is_complex(typecode) ((typecode)[2]=='C')
#define mm_is_real(typecode) ((typecode)[2]=='R')
#define mm_is_pattern(typecode) ((typecode)[2]=='P')
#define mm_is_integer(typecode) ((typecode)[2]=='I')
#define mm_is_symmetric(typecode)((typecode)[3]=='S')
#define mm_is_general(typecode) ((typecode)[3]=='G')
#define mm_is_skew(typecode) ((typecode)[3]=='K')
#define mm_is_hermitian(typecode)((typecode)[3]=='H')
int mm_is_valid(MM_typecode matcode); /* too complex for a macro */
/********************* MM_typecode modify fucntions ***************************/
#define mm_set_matrix(typecode) ((*typecode)[0]='M')
#define mm_set_coordinate(typecode) ((*typecode)[1]='C')
#define mm_set_array(typecode) ((*typecode)[1]='A')
#define mm_set_dense(typecode) mm_set_array(typecode)
#define mm_set_sparse(typecode) mm_set_coordinate(typecode)
#define mm_set_complex(typecode)((*typecode)[2]='C')
#define mm_set_real(typecode) ((*typecode)[2]='R')
#define mm_set_pattern(typecode)((*typecode)[2]='P')
#define mm_set_integer(typecode)((*typecode)[2]='I')
#define mm_set_symmetric(typecode)((*typecode)[3]='S')
#define mm_set_general(typecode)((*typecode)[3]='G')
#define mm_set_skew(typecode) ((*typecode)[3]='K')
#define mm_set_hermitian(typecode)((*typecode)[3]='H')
#define mm_clear_typecode(typecode) ((*typecode)[0]=(*typecode)[1]= \
(*typecode)[2]=' ',(*typecode)[3]='G')
#define mm_initialize_typecode(typecode) mm_clear_typecode(typecode)
/********************* Matrix Market error codes ***************************/
#define MM_COULD_NOT_READ_FILE 11
#define MM_PREMATURE_EOF 12
#define MM_NOT_MTX 13
#define MM_NO_HEADER 14
#define MM_UNSUPPORTED_TYPE 15
#define MM_LINE_TOO_LONG 16
#define MM_COULD_NOT_WRITE_FILE 17
/******************** Matrix Market internal definitions ********************
MM_matrix_typecode: 4-character sequence
ojbect sparse/ data storage
dense type scheme
string position: [0] [1] [2] [3]
Matrix typecode: M(atrix) C(oord) R(eal) G(eneral)
A(array) C(omplex) H(ermitian)
P(attern) S(ymmetric)
I(nteger) K(kew)
***********************************************************************/
#define MM_MTX_STR "matrix"
#define MM_ARRAY_STR "array"
#define MM_DENSE_STR "array"
#define MM_COORDINATE_STR "coordinate"
#define MM_SPARSE_STR "coordinate"
#define MM_COMPLEX_STR "complex"
#define MM_REAL_STR "real"
#define MM_INT_STR "integer"
#define MM_GENERAL_STR "general"
#define MM_SYMM_STR "symmetric"
#define MM_HERM_STR "hermitian"
#define MM_SKEW_STR "skew-symmetric"
#define MM_PATTERN_STR "pattern"
/* high level routines */
int mm_write_mtx_crd(char fname[], int M, int N, int nz, int I[], int J[],
double val[], MM_typecode matcode);
int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int I[], int J[],
double val[], MM_typecode matcode);
int mm_read_mtx_crd_entry(FILE *f, int *I, int *J, double *real, double *img,
MM_typecode matcode);
int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_,
double **val_, int **I_, int **J_);
#endif

825
gtsam/3rdparty/lp_solve_5.5/shared/myblas.c vendored
View File

@ -1,825 +0,0 @@
#include <stdlib.h>
#include <stdio.h>
/*#include <memory.h>*/
#include <string.h>
#include <math.h>
#include "myblas.h"
#ifdef FORTIFY
# include "lp_fortify.h"
#endif
/* ************************************************************************ */
/* Initialize BLAS interfacing routines */
/* ************************************************************************ */
MYBOOL mustinitBLAS = TRUE;
#ifdef WIN32
HINSTANCE hBLAS = NULL;
#else
void *hBLAS = NULL;
#endif
/* ************************************************************************ */
/* Function pointers for external BLAS library (C base 0) */
/* ************************************************************************ */
BLAS_dscal_func *BLAS_dscal;
BLAS_dcopy_func *BLAS_dcopy;
BLAS_daxpy_func *BLAS_daxpy;
BLAS_dswap_func *BLAS_dswap;
BLAS_ddot_func *BLAS_ddot;
BLAS_idamax_func *BLAS_idamax;
BLAS_dload_func *BLAS_dload;
BLAS_dnormi_func *BLAS_dnormi;
/* ************************************************************************ */
/* Define the BLAS interfacing routines */
/* ************************************************************************ */
void init_BLAS(void)
{
if(mustinitBLAS) {
load_BLAS(NULL);
mustinitBLAS = FALSE;
}
}
MYBOOL is_nativeBLAS(void)
{
#ifdef LoadableBlasLib
return( (MYBOOL) (hBLAS == NULL) );
#else
return( TRUE );
#endif
}
MYBOOL load_BLAS(char *libname)
{
MYBOOL result = TRUE;
#ifdef LoadableBlasLib
if(hBLAS != NULL) {
#ifdef WIN32
FreeLibrary(hBLAS);
#else
dlclose(hBLAS);
#endif
hBLAS = NULL;
}
#endif
if(libname == NULL) {
if(!mustinitBLAS && is_nativeBLAS())
return( FALSE );
BLAS_dscal = my_dscal;
BLAS_dcopy = my_dcopy;
BLAS_daxpy = my_daxpy;
BLAS_dswap = my_dswap;
BLAS_ddot = my_ddot;
BLAS_idamax = my_idamax;
BLAS_dload = my_dload;
BLAS_dnormi = my_dnormi;
if(mustinitBLAS)
mustinitBLAS = FALSE;
}
else {
#ifdef LoadableBlasLib
#ifdef WIN32
/* Get a handle to the Windows DLL module. */
hBLAS = LoadLibrary(libname);
/* If the handle is valid, try to get the function addresses. */
result = (MYBOOL) (hBLAS != NULL);
if(result) {
BLAS_dscal = (BLAS_dscal_func *) GetProcAddress(hBLAS, BLAS_prec "scal");
BLAS_dcopy = (BLAS_dcopy_func *) GetProcAddress(hBLAS, BLAS_prec "copy");
BLAS_daxpy = (BLAS_daxpy_func *) GetProcAddress(hBLAS, BLAS_prec "axpy");
BLAS_dswap = (BLAS_dswap_func *) GetProcAddress(hBLAS, BLAS_prec "swap");
BLAS_ddot = (BLAS_ddot_func *) GetProcAddress(hBLAS, BLAS_prec "dot");
BLAS_idamax = (BLAS_idamax_func *) GetProcAddress(hBLAS, "i" BLAS_prec "amax");
#if 0
BLAS_dload = (BLAS_dload_func *) GetProcAddress(hBLAS, BLAS_prec "load");
BLAS_dnormi = (BLAS_dnormi_func *) GetProcAddress(hBLAS, BLAS_prec "normi");
#endif
}
#else
/* First standardize UNIX .SO library name format. */
char blasname[260], *ptr;
strcpy(blasname, libname);
if((ptr = strrchr(libname, '/')) == NULL)
ptr = libname;
else
ptr++;
blasname[(int) (ptr - libname)] = 0;
if(strncmp(ptr, "lib", 3))
strcat(blasname, "lib");
strcat(blasname, ptr);
if(strcmp(blasname + strlen(blasname) - 3, ".so"))
strcat(blasname, ".so");
/* Get a handle to the module. */
hBLAS = dlopen(blasname, RTLD_LAZY);
/* If the handle is valid, try to get the function addresses. */
result = (MYBOOL) (hBLAS != NULL);
if(result) {
BLAS_dscal = (BLAS_dscal_func *) dlsym(hBLAS, BLAS_prec "scal");
BLAS_dcopy = (BLAS_dcopy_func *) dlsym(hBLAS, BLAS_prec "copy");
BLAS_daxpy = (BLAS_daxpy_func *) dlsym(hBLAS, BLAS_prec "axpy");
BLAS_dswap = (BLAS_dswap_func *) dlsym(hBLAS, BLAS_prec "swap");
BLAS_ddot = (BLAS_ddot_func *) dlsym(hBLAS, BLAS_prec "dot");
BLAS_idamax = (BLAS_idamax_func *) dlsym(hBLAS, "i" BLAS_prec "amax");
#if 0
BLAS_dload = (BLAS_dload_func *) dlsym(hBLAS, BLAS_prec "load");
BLAS_dnormi = (BLAS_dnormi_func *) dlsym(hBLAS, BLAS_prec "normi");
#endif
}
#endif
#endif
/* Do validation */
if(!result ||
((BLAS_dscal == NULL) ||
(BLAS_dcopy == NULL) ||
(BLAS_daxpy == NULL) ||
(BLAS_dswap == NULL) ||
(BLAS_ddot == NULL) ||
(BLAS_idamax == NULL) ||
(BLAS_dload == NULL) ||
(BLAS_dnormi == NULL))
) {
load_BLAS(NULL);
result = FALSE;
}
}
return( result );
}
MYBOOL unload_BLAS(void)
{
return( load_BLAS(NULL) );
}
/* ************************************************************************ */
/* Now define the unoptimized local BLAS functions */
/* ************************************************************************ */
void daxpy( int n, REAL da, REAL *dx, int incx, REAL *dy, int incy)
{
dx++;
dy++;
BLAS_daxpy( &n, &da, dx, &incx, dy, &incy);
}
void BLAS_CALLMODEL my_daxpy( int *_n, REAL *_da, REAL *dx, int *_incx, REAL *dy, int *_incy)
{
/* constant times a vector plus a vector.
uses unrolled loops for increments equal to one.
jack dongarra, linpack, 3/11/78.
modified 12/3/93, array[1] declarations changed to array[*] */
int i, ix, iy, m, mp1;
register REAL rda;
REAL da = *_da;
int n = *_n, incx = *_incx, incy = *_incy;
if (n <= 0) return;
if (da == 0.0) return;
dx--;
dy--;
ix = 1;
iy = 1;
if (incx < 0)
ix = (-n+1)*incx + 1;
if (incy < 0)
iy = (-n+1)*incy + 1;
rda = da;
/* CPU intensive loop; option to do pointer arithmetic */
#if defined DOFASTMATH
{
REAL *xptr, *yptr;
for (i = 1, xptr = dx + ix, yptr = dy + iy;
i <= n; i++, xptr += incx, yptr += incy)
(*yptr) += rda*(*xptr);
return;
}
#else
if (incx==1 && incy==1) goto x20;
/* code for unequal increments or equal increments not equal to 1 */
for (i = 1; i<=n; i++) {
dy[iy]+= rda*dx[ix];
ix+= incx;
iy+= incy;
}
return;
/* code for both increments equal to 1 */
/* clean-up loop */
x20:
m = n % 4;
if (m == 0) goto x40;
for (i = 1; i<=m; i++)
dy[i]+= rda*dx[i];
if(n < 4) return;
x40:
mp1 = m + 1;
for (i = mp1; i<=n; i=i+4) {
dy[i]+= rda*dx[i];
dy[i + 1]+= rda*dx[i + 1];
dy[i + 2]+= rda*dx[i + 2];
dy[i + 3]+= rda*dx[i + 3];
}
#endif
}
/* ************************************************************************ */
void dcopy( int n, REAL *dx, int incx, REAL *dy, int incy)
{
dx++;
dy++;
BLAS_dcopy( &n, dx, &incx, dy, &incy);
}
void BLAS_CALLMODEL my_dcopy (int *_n, REAL *dx, int *_incx, REAL *dy, int *_incy)
{
/* copies a vector, x, to a vector, y.
uses unrolled loops for increments equal to one.
jack dongarra, linpack, 3/11/78.
modified 12/3/93, array[1] declarations changed to array[*] */
int i, ix, iy, m, mp1;
int n = *_n, incx = *_incx, incy = *_incy;
if (n<=0) return;
dx--;
dy--;
ix = 1;
iy = 1;
if (incx<0)
ix = (-n+1)*incx + 1;
if (incy<0)
iy = (-n+1)*incy + 1;
/* CPU intensive loop; option to do pointer arithmetic */
#if defined DOFASTMATH
{
REAL *xptr, *yptr;
for (i = 1, xptr = dx + ix, yptr = dy + iy;
i <= n; i++, xptr += incx, yptr += incy)
(*yptr) = (*xptr);
return;
}
#else
if (incx==1 && incy==1) goto x20;
/* code for unequal increments or equal increments not equal to 1 */
for (i = 1; i<=n; i++) {
dy[iy] = dx[ix];
ix+= incx;
iy+= incy;
}
return;
/* code for both increments equal to 1 */
/* version with fast machine copy logic (requires memory.h or string.h) */
x20:
#if defined DOFASTMATH
MEMCOPY(&dy[1], &dx[1], n);
return;
#else
m = n % 7;
if (m == 0) goto x40;
for (i = 1; i<=m; i++)
dy[i] = dx[i];
if (n < 7) return;
x40:
mp1 = m + 1;
for (i = mp1; i<=n; i=i+7) {
dy[i] = dx[i];
dy[i + 1] = dx[i + 1];
dy[i + 2] = dx[i + 2];
dy[i + 3] = dx[i + 3];
dy[i + 4] = dx[i + 4];
dy[i + 5] = dx[i + 5];
dy[i + 6] = dx[i + 6];
}
#endif
#endif
}
/* ************************************************************************ */
void dscal (int n, REAL da, REAL *dx, int incx)
{
dx++;
BLAS_dscal (&n, &da, dx, &incx);
}
void BLAS_CALLMODEL my_dscal (int *_n, REAL *_da, REAL *dx, int *_incx)
{
/* Multiply a vector by a constant.
--Input--
N number of elements in input vector(s)
DA double precision scale factor
DX double precision vector with N elements
INCX storage spacing between elements of DX
--Output--
DX double precision result (unchanged if N.LE.0)
Replace double precision DX by double precision DA*DX.
For I = 0 to N-1, replace DX(IX+I*INCX) with DA * DX(IX+I*INCX),
where IX = 1 if INCX .GE. 0, else IX = 1+(1-N)*INCX. */
int i, ix, m, mp1;
register REAL rda;
REAL da = *_da;
int n = *_n, incx = *_incx;
if (n <= 0)
return;
rda = da;
dx--;
/* Optionally do fast pointer arithmetic */
#if defined DOFASTMATH
{
REAL *xptr;
for (i = 1, xptr = dx + 1; i <= n; i++, xptr += incx)
(*xptr) *= rda;
return;
}
#else
if (incx == 1)
goto x20;
/* Code for increment not equal to 1 */
ix = 1;
if (incx < 0)
ix = (-n+1)*incx + 1;
for(i = 1; i <= n; i++, ix += incx)
dx[ix] *= rda;
return;
/* Code for increment equal to 1. */
/* Clean-up loop so remaining vector length is a multiple of 5. */
x20:
m = n % 5;
if (m == 0) goto x40;
for( i = 1; i <= m; i++)
dx[i] *= rda;
if (n < 5)
return;
x40:
mp1 = m + 1;
for(i = mp1; i <= n; i += 5) {
dx[i] *= rda;
dx[i+1] *= rda;
dx[i+2] *= rda;
dx[i+3] *= rda;
dx[i+4] *= rda;
}
#endif
}
/* ************************************************************************ */
REAL ddot(int n, REAL *dx, int incx, REAL *dy, int incy)
{
dx++;
dy++;
return( BLAS_ddot (&n, dx, &incx, dy, &incy) );
}
REAL BLAS_CALLMODEL my_ddot(int *_n, REAL *dx, int *_incx, REAL *dy, int *_incy)
{
/* forms the dot product of two vectors.
uses unrolled loops for increments equal to one.
jack dongarra, linpack, 3/11/78.
modified 12/3/93, array[1] declarations changed to array[*] */
register REAL dtemp;
int i, ix, iy, m, mp1;
int n = *_n, incx = *_incx, incy = *_incy;
dtemp = 0.0;
if (n<=0)
return( (REAL) dtemp);
dx--;
dy--;
ix = 1;
iy = 1;
if (incx<0)
ix = (-n+1)*incx + 1;
if (incy<0)
iy = (-n+1)*incy + 1;
/* CPU intensive loop; option to do pointer arithmetic */
#if defined DOFASTMATH
{
REAL *xptr, *yptr;
for (i = 1, xptr = dx + ix, yptr = dy + iy;
i <= n; i++, xptr += incx, yptr += incy)
dtemp+= (*yptr)*(*xptr);
return(dtemp);
}
#else
if (incx==1 && incy==1) goto x20;
/* code for unequal increments or equal increments not equal to 1 */
for (i = 1; i<=n; i++) {
dtemp+= dx[ix]*dy[iy];
ix+= incx;
iy+= incy;
}
return(dtemp);
/* code for both increments equal to 1 */
/* clean-up loop */
x20:
m = n % 5;
if (m == 0) goto x40;
for (i = 1; i<=m; i++)
dtemp+= dx[i]*dy[i];
if (n < 5) goto x60;
x40:
mp1 = m + 1;
for (i = mp1; i<=n; i=i+5)
dtemp+= dx[i]*dy[i] + dx[i + 1]*dy[i + 1] +
dx[i + 2]*dy[i + 2] + dx[i + 3]*dy[i + 3] + dx[i + 4]*dy[i + 4];
x60:
return(dtemp);
#endif
}
/* ************************************************************************ */
void dswap( int n, REAL *dx, int incx, REAL *dy, int incy )
{
dx++;
dy++;
BLAS_dswap( &n, dx, &incx, dy, &incy );
}
void BLAS_CALLMODEL my_dswap( int *_n, REAL *dx, int *_incx, REAL *dy, int *_incy )
{
int i, ix, iy, m, mp1, ns;
REAL dtemp1, dtemp2, dtemp3;
int n = *_n, incx = *_incx, incy = *_incy;
if (n <= 0) return;
dx--;
dy--;
ix = 1;
iy = 1;
if (incx < 0)
ix = (-n+1)*incx + 1;
if (incy < 0)
iy = (-n+1)*incy + 1;
/* CPU intensive loop; option to do pointer arithmetic */
#if defined DOFASTMATH
{
REAL *xptr, *yptr;
for (i = 1, xptr = dx + ix, yptr = dy + iy;
i <= n; i++, xptr += incx, yptr += incy) {
dtemp1 = (*xptr);
(*xptr) = (*yptr);
(*yptr) = dtemp1;
}
return;
}
#else
if (incx == incy) {
if (incx <= 0) goto x5;
if (incx == 1) goto x20;
goto x60;
}
/* code for unequal or nonpositive increments. */
x5:
for (i = 1; i<=n; i++) {
dtemp1 = dx[ix];
dx[ix] = dy[iy];
dy[iy] = dtemp1;
ix+= incx;
iy+= incy;
}
return;
/* code for both increments equal to 1.
clean-up loop so remaining vector length is a multiple of 3. */
x20:
m = n % 3;
if (m == 0) goto x40;
for (i = 1; i<=m; i++) {
dtemp1 = dx[i];
dx[i] = dy[i];
dy[i] = dtemp1;
}
if (n < 3) return;
x40:
mp1 = m + 1;
for (i = mp1; i<=n; i=i+3) {
dtemp1 = dx[i];
dtemp2 = dx[i+1];
dtemp3 = dx[i+2];
dx[i] = dy[i];
dx[i+1] = dy[i+1];
dx[i+2] = dy[i+2];
dy[i] = dtemp1;
dy[i+1] = dtemp2;
dy[i+2] = dtemp3;
}
return;
/* code for equal, positive, non-unit increments. */
x60:
ns = n*incx;
for (i = 1; i<=ns; i=i+incx) {
dtemp1 = dx[i];
dx[i] = dy[i];
dy[i] = dtemp1;
}
#endif
}
/* ************************************************************************ */
void dload(int n, REAL da, REAL *dx, int incx)
{
dx++;
BLAS_dload (&n, &da, dx, &incx);
}
void BLAS_CALLMODEL my_dload (int *_n, REAL *_da, REAL *dx, int *_incx)
{
/* copies a scalar, a, to a vector, x.
uses unrolled loops when incx equals one.
To change the precision of this program, run the change
program on dload.f
Alternatively, to make a single precision version append a
comment character to the start of all lines between sequential
precision > double
and
end precision > double
comments and delete the comment character at the start of all
lines between sequential
precision > single
and
end precision > single
comments. To make a double precision version interchange
the append and delete operations in these instructions. */
int i, ix, m, mp1;
REAL da = *_da;
int n = *_n, incx = *_incx;
if (n<=0) return;
dx--;
if (incx==1) goto x20;
/* code for incx not equal to 1 */
ix = 1;
if (incx<0)
ix = (-n+1)*incx + 1;
for (i = 1; i<=n; i++) {
dx[ix] = da;
ix+= incx;
}
return;
/* code for incx equal to 1 and clean-up loop */
x20:
m = n % 7;
if (m == 0) goto x40;
for (i = 1; i<=m; i++)
dx[i] = da;
if (n < 7) return;
x40:
mp1 = m + 1;
for (i = mp1; i<=n; i=i+7) {
dx[i] = da;
dx[i + 1] = da;
dx[i + 2] = da;
dx[i + 3] = da;
dx[i + 4] = da;
dx[i + 5] = da;
dx[i + 6] = da;
}
}
/* ************************************************************************ */
int idamax( int n, REAL *x, int is )
{
x++;
return ( BLAS_idamax( &n, x, &is ) );
}
int BLAS_CALLMODEL my_idamax( int *_n, REAL *x, int *_is )
{
register REAL xmax, xtest;
int i, imax = 0;
#if !defined DOFASTMATH
int ii;
#endif
int n = *_n, is = *_is;
if((n < 1) || (is <= 0))
return(imax);
imax = 1;
if(n == 1)
return(imax);
#if defined DOFASTMATH
xmax = fabs(*x);
for (i = 2, x += is; i <= n; i++, x += is) {
xtest = fabs(*x);
if(xtest > xmax) {
xmax = xtest;
imax = i;
}
}
#else
x--;
ii = 1;
xmax = fabs(x[ii]);
for(i = 2, ii+ = is; i <= n; i++, ii+ = is) {
xtest = fabs(x[ii]);
if(xtest > xmax) {
xmax = xtest;
imax = i;
}
}
#endif
return(imax);
}
/* ************************************************************************ */
REAL dnormi( int n, REAL *x )
{
x++;
return( BLAS_dnormi( &n, x ) );
}
REAL BLAS_CALLMODEL my_dnormi( int *_n, REAL *x )
{
/* ===============================================================
dnormi returns the infinity-norm of the vector x.
=============================================================== */
int j;
register REAL hold, absval;
int n = *_n;
x--;
hold = 0.0;
/* for(j = 1; j <= n; j++) */
for(j = n; j > 0; j--) {
absval = fabs(x[j]);
hold = MAX( hold, absval );
}
return( hold );
}
/* ************************************************************************ */
/* Subvector and submatrix access routines (Fortran compatibility) */
/* ************************************************************************ */
#ifndef UseMacroVector
int subvec( int item)
{
return( item-1 );
}
#endif
int submat( int nrowb, int row, int col)
{
return( nrowb*(col-1) + subvec(row) );
}
int posmat( int nrowb, int row, int col)
{
return( submat(nrowb, row, col)+BLAS_BASE );
}
/* ************************************************************************ */
/* Randomization functions */
/* ************************************************************************ */
void randomseed(int seeds[])
/* Simply create some default seed values */
{
seeds[1] = 123456;
seeds[2] = 234567;
seeds[3] = 345678;
}
void randomdens( int n, REAL *x, REAL r1, REAL r2, REAL densty, int *seeds )
{
/* ------------------------------------------------------------------
random generates a vector x[*] of random numbers
in the range (r1, r2) with (approximate) specified density.
seeds[*] must be initialized before the first call.
------------------------------------------------------------------ */
int i;
REAL *y;
y = (REAL *) malloc(sizeof(*y) * (n+1));
ddrand( n, x, 1, seeds );
ddrand( n, y, 1, seeds );
for (i = 1; i<=n; i++) {
if (y[i] < densty)
x[i] = r1 + (r2 - r1) * x[i];
else
x[i] = 0.0;
}
free(y);
}
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
void ddrand( int n, REAL *x, int incx, int *seeds )
{
/* ------------------------------------------------------------------
ddrand fills a vector x with uniformly distributed random numbers
in the interval (0, 1) using a method due to Wichman and Hill.
seeds[1..3] should be set to integer values
between 1 and 30000 before the first entry.
Integer arithmetic up to 30323 is required.
Blatantly copied from Wichman and Hill 19-January-1987.
14-Feb-94. Original version.
30 Jun 1999. seeds stored in an array.
30 Jun 1999. This version of ddrand.
------------------------------------------------------------------ */
int ix, xix;
if (n < 1) return;
for (ix = 1; ix<=1+(n-1)*incx; ix=ix+incx) {
seeds[1] = 171*(seeds[1] % 177) - 2*(seeds[1]/177);
seeds[2] = 172*(seeds[2] % 176) - 35*(seeds[2]/176);
seeds[3] = 170*(seeds[3] % 178) - 63*(seeds[3]/178);
if (seeds[1] < 0) seeds[1] = seeds[1] + 30269;
if (seeds[2] < 0) seeds[2] = seeds[2] + 30307;
if (seeds[3] < 0) seeds[3] = seeds[3] + 30323;
x[ix] = ((REAL) seeds[1])/30269.0 +
((REAL) seeds[2])/30307.0 +
((REAL) seeds[3])/30323.0;
xix = (int) x[ix];
x[ix] = fabs(x[ix] - xix);
}
}

132
gtsam/3rdparty/lp_solve_5.5/shared/myblas.h vendored
View File

@ -1,132 +0,0 @@
#ifndef HEADER_myblas
#define HEADER_myblas
/* ************************************************************************ */
/* BLAS function interface with local and external loadable versions */
/* Author: Kjell Eikland */
/* Version: Initial version spring 2004 */
/* Licence: LGPL */
/* ************************************************************************ */
/* Changes: 19 September 2004 Moved function pointer variable */
/* declarations from myblas.h to myblas.c */
/* to avoid linker problems with the Mac. */
/* 20 April 2005 Modified all double types to REAL to self- */
/* adjust to global settings. Note that BLAS */
/* as of now does not have double double. */
/* ************************************************************************ */
#define BLAS_BASE 1
#define UseMacroVector
#if defined LoadableBlasLib
# if LoadableBlasLib == 0
# undef LoadableBlasLib
# endif
#else
# define LoadableBlasLib
#endif
/* ************************************************************************ */
/* Include necessary libraries */
/* ************************************************************************ */
#include "commonlib.h"
#ifdef LoadableBlasLib
#ifdef WIN32
#include <windows.h>
#else
#include <dlfcn.h>
#endif
#endif
#ifdef __cplusplus
extern "C" {
#endif
/* ************************************************************************ */
/* BLAS functions */
/* ************************************************************************ */
#ifndef BLAS_CALLMODEL
#ifdef WIN32
# define BLAS_CALLMODEL _cdecl
#else
# define BLAS_CALLMODEL
#endif
#endif
typedef void (BLAS_CALLMODEL BLAS_dscal_func) (int *n, REAL *da, REAL *dx, int *incx);
typedef void (BLAS_CALLMODEL BLAS_dcopy_func) (int *n, REAL *dx, int *incx, REAL *dy, int *incy);
typedef void (BLAS_CALLMODEL BLAS_daxpy_func) (int *n, REAL *da, REAL *dx, int *incx, REAL *dy, int *incy);
typedef void (BLAS_CALLMODEL BLAS_dswap_func) (int *n, REAL *dx, int *incx, REAL *dy, int *incy);
typedef double (BLAS_CALLMODEL BLAS_ddot_func) (int *n, REAL *dx, int *incx, REAL *dy, int *incy);
typedef int (BLAS_CALLMODEL BLAS_idamax_func)(int *n, REAL *x, int *is);
typedef void (BLAS_CALLMODEL BLAS_dload_func) (int *n, REAL *da, REAL *dx, int *incx);
typedef double (BLAS_CALLMODEL BLAS_dnormi_func)(int *n, REAL *x);
#ifndef __WINAPI
# ifdef WIN32
# define __WINAPI WINAPI
# else
# define __WINAPI
# endif
#endif
void init_BLAS(void);
MYBOOL is_nativeBLAS(void);
MYBOOL load_BLAS(char *libname);
MYBOOL unload_BLAS(void);
/* ************************************************************************ */
/* User-callable BLAS definitions (C base 1) */
/* ************************************************************************ */
void dscal ( int n, REAL da, REAL *dx, int incx );
void dcopy ( int n, REAL *dx, int incx, REAL *dy, int incy );
void daxpy ( int n, REAL da, REAL *dx, int incx, REAL *dy, int incy );
void dswap ( int n, REAL *dx, int incx, REAL *dy, int incy );
REAL ddot ( int n, REAL *dx, int incx, REAL *dy, int incy );
int idamax( int n, REAL *x, int is );
void dload ( int n, REAL da, REAL *dx, int incx );
REAL dnormi( int n, REAL *x );
/* ************************************************************************ */
/* Locally implemented BLAS functions (C base 0) */
/* ************************************************************************ */
void BLAS_CALLMODEL my_dscal ( int *n, REAL *da, REAL *dx, int *incx );
void BLAS_CALLMODEL my_dcopy ( int *n, REAL *dx, int *incx, REAL *dy, int *incy );
void BLAS_CALLMODEL my_daxpy ( int *n, REAL *da, REAL *dx, int *incx, REAL *dy, int *incy );
void BLAS_CALLMODEL my_dswap ( int *n, REAL *dx, int *incx, REAL *dy, int *incy );
REAL BLAS_CALLMODEL my_ddot ( int *n, REAL *dx, int *incx, REAL *dy, int *incy );
int BLAS_CALLMODEL my_idamax( int *n, REAL *x, int *is );
void BLAS_CALLMODEL my_dload ( int *n, REAL *da, REAL *dx, int *incx );
REAL BLAS_CALLMODEL my_dnormi( int *n, REAL *x );
/* ************************************************************************ */
/* Subvector and submatrix access routines (Fortran compatibility) */
/* ************************************************************************ */
#ifdef UseMacroVector
#define subvec(item) (item - 1)
#else
int subvec( int item );
#endif
int submat( int nrowb, int row, int col );
int posmat( int nrowb, int row, int col );
/* ************************************************************************ */
/* Randomization functions */
/* ************************************************************************ */
void randomseed(int *seeds);
void randomdens( int n, REAL *x, REAL r1, REAL r2, REAL densty, int *seeds);
void ddrand( int n, REAL *x, int incx, int *seeds );
#ifdef __cplusplus
}
#endif
#endif