Commits

Shlomi Fish committed ceff58b

Add the shared sources.

Comments (0)

Files changed (6)

shared/commonlib.c

+
+#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, "%s", 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, "%s", 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, "%s", 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

shared/commonlib.h

+#ifndef HEADER_commonlib
+#define HEADER_commonlib
+
+#include <stdlib.h>
+#include <stdio.h>
+#ifdef WIN32
+  #include <windows.h>
+#endif
+
+/* 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 18446744073709551615ll
+  #else
+    #define MAXUINT64 18446744073709551615l
+  #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 */
+/*
+*   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;
+    int x;
+
+    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++)
+    {
+        x = 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];
+
+}
+/*
+*   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
+
+#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;
+#if !defined DOFASTMATH
+  int      m, mp1;
+#endif
+  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;
+#if !defined DOFASTMATH
+  int      m, mp1;
+#endif
+  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;
+#if !defined DOFASTMATH
+  int      ix, m, mp1;
+#endif
+  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)