#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "commonlib.h"
#include "myblas.h"
#include "sparselib.h"


sparseMatrix *createMatrix(int dimLimit, int lenLimit, int initVectors)
{
  int          initsize;
  sparseMatrix *matrix;

  if(initVectors < 0)
    initVectors = 0;
  if(initVectors == 0)
    initsize = MIN(INITIALSIZE, dimLimit);
  else
    initsize = MAX(INITIALSIZE, initVectors);

  CALLOC(matrix, 1);
  matrix->limit = dimLimit;
  matrix->limitVector = lenLimit;
  resizeMatrix(matrix, initsize);
  while(initVectors > 0) {
    initVectors--;
    appendMatrix(matrix, createVector(lenLimit, 2));
  }
  return(matrix);
}


void resizeMatrix(sparseMatrix *matrix, int newSize)
{
  int oldSize;

  if(matrix == NULL)
    oldSize = 0;
  else
    oldSize = matrix->size;
  while(oldSize>newSize) {
	  oldSize--;
	  freeVector(matrix->list[oldSize]);
    return;
  }
  REALLOC(matrix->list, newSize);
  while(oldSize<newSize) {
	  matrix->list[oldSize] = NULL;
	  oldSize++;
  }
  if(newSize>0)
    matrix->size = newSize;
}

int appendMatrix(sparseMatrix *matrix, sparseVector *newVector)
{
   if(matrix->count == matrix->size)
     resizeMatrix(matrix, matrix->size + 10);
   matrix->list[matrix->count] = newVector;
   matrix->count++;
   putDiagonalIndex(newVector, matrix->count);
   return(matrix->count);
}


int NZcountMatrix(sparseMatrix *matrix)
{
  int i, nz;

  nz = 0;
  for(i = 0; i < matrix->count; i++)
    nz += matrix->list[i]->count;

  return( nz );
}


void freeMatrix(sparseMatrix *matrix)
{
  resizeMatrix(matrix, 0);
  FREE(matrix);
}


void printMatrix(int n, sparseMatrix *matrix, int modulo, MYBOOL showEmpty)
{
   int i;
   for(i = 1; i<=matrix->count; i++) 
     if(matrix->list[i-1] != NULL) { 
       if(showEmpty || matrix->list[i-1]->count>0)
         printVector(n, matrix->list[i-1], modulo);
   }
}


sparseVector *createVector(int dimLimit, int initSize)
{
  sparseVector *newitem;
  CALLOC(newitem, 1);
  newitem->limit = dimLimit;
  initSize = resizeVector(newitem, initSize);
  return(newitem);
}


sparseVector *cloneVector(sparseVector *sparse)
{
  sparseVector *hold;
  hold = createVector(sparse->limit, sparse->count);
  hold->count = sparse->count;
  MEMCOPY(&hold->value[0], &sparse->value[0], (sparse->count+1));
  MEMCOPY(&hold->index[0], &sparse->index[0], (sparse->count+1));
  return(hold);
}

int redimensionVector(sparseVector *sparse, int newDim)
{
  int olddim, i;

  olddim = sparse->limit;
  sparse->limit = newDim;
  if(lastIndex(sparse)>newDim) {
    i = sparse->count;
    while(i>0 && sparse->index[i]>newDim) i--;
    sparse->count = i;
    resizeVector(sparse, sparse->count);
  }
  return(olddim);
}


int resizeVector(sparseVector *sparse, int newSize)
{
  int oldsize;

  oldsize = sparse->size;
  REALLOC(sparse->value, (newSize+1));
  REALLOC(sparse->index, (newSize+1));
  sparse->size = newSize;
  return(oldsize);
}


void moveVector(sparseVector *sparse, int destPos, int sourcePos, int itemCount)
{
  int i;
  
  if(itemCount <= 0 || sourcePos == destPos)
    return;

#if defined DOFASTMATH
  if(TRUE) {
    MEMMOVE(&sparse->value[destPos], &sparse->value[sourcePos], itemCount);
    MEMMOVE(&sparse->index[destPos], &sparse->index[sourcePos], itemCount);
  }
  else {
    int    *idxPtr1, *idxPtr2;
    double *valPtr1, *valPtr2;

    for(i = 1, idxPtr1 = sparse->index+destPos, idxPtr2 = sparse->index+sourcePos,
               valPtr1 = sparse->value+destPos, valPtr2 = sparse->value+sourcePos; 
        i<=itemCount; i++, idxPtr1++, idxPtr2++, valPtr1++, valPtr2++) {
      *idxPtr1 = *idxPtr2;
      *valPtr1 = *valPtr2;
    }
  }
#else
  for(i = 1; i<=itemCount; i++) {
    sparse->value[destPos] = sparse->value[sourcePos];
    sparse->index[destPos] = sparse->index[sourcePos];
    destPos++;
    sourcePos++;
  }
#endif
}


void rotateVector(sparseVector *sparse, int startPos, int chainSize, int stepDelta)
{
/*  int idxHold; */
/*  double valHold; */

}


void swapVector(sparseVector *sparse1, sparseVector *sparse2)
{
  int n, m, *idx;
  REAL *val;

  n = sparse1->count;
  sparse1->count = sparse2->count;
  sparse2->count = n;

  n = sparse1->size;
  sparse1->size = sparse2->size;
  sparse2->size = n;

  n = sparse1->limit;
  sparse1->limit = sparse2->limit;
  sparse2->limit = n;

  idx = sparse1->index;
  sparse1->index = sparse2->index;
  sparse2->index = idx;

  val = sparse1->value;
  sparse1->value = sparse2->value;
  sparse2->value = val;

  n = getDiagonalIndex(sparse1);
  m = getDiagonalIndex(sparse2);
  putDiagonalIndex(sparse1, m);
  putDiagonalIndex(sparse2, n);

}


void freeVector(sparseVector *sparse)
{
	if(sparse != NULL) {
    FREE(sparse->value);
    FREE(sparse->index);
    FREE(sparse);
  }
}


MYBOOL verifyVector(sparseVector *sparse)
{
  int i, n, k1, k2, kd; 
  int  err = 0;
  double vd;

  n = sparse->count;
  kd = sparse->index[0];
  vd = sparse->value[0];
  if(n <= 1)
    return(TRUE);
  k1 = 0;
  k2 = sparse->index[1];
  if(k2 == kd && sparse->value[1] != vd) 
    err = 2;

  for(i = 2; i <= n && err == 0; i++) {
    k1 = k2;
    k2 = sparse->index[i];
    if(k1 >= k2) err = 1;
    if(k2 == kd && sparse->value[i] != vd) err = 2;
  }
  if(err == 0)
    return(TRUE);
  else if(err == 1)
    printf("Invalid sparse vector index order");
  else if(err == 2)
    printf("Invalid sparse vector diagonal value");
  return(FALSE);
}


int firstIndex(sparseVector *sparse)
{
  return(sparse->index[1]);
}


int lastIndex(sparseVector *sparse)
{
  return(sparse->index[sparse->count]);
}


int getDiagonalIndex(sparseVector *sparse)
{
  return(sparse->index[0]);
}


int putDiagonalIndex(sparseVector *sparse, int index)
{
  int oldindex;
  oldindex = sparse->index[0];
  if(index > 0) {
    sparse->index[0] = 0; /* Must temporarily set to zero to force vector search in getItem */
    sparse->value[0] = getItem(sparse, index);
  }
  else
    sparse->value[0] = 0;
  sparse->index[0] = index;
  return(oldindex);
}


MYBOOL putDiagonal(sparseVector *sparse, REAL value)
{
  if(sparse->index[0]>0) {
    putItem(sparse, sparse->index[0], value); 
    return(TRUE);
  }
  else
    return(FALSE);
}


REAL getDiagonal(sparseVector *sparse)
{
   return(sparse->value[0]);
}


REAL getItem(sparseVector *sparse, int targetIndex)
{
  /* First check if we want the diagonal element */ 
  if(targetIndex == sparse->index[0])
    return(sparse->value[0]);

  /* If not, search for the variable's position in the index list */
  targetIndex = findIndex(targetIndex, sparse->index, sparse->count, BLAS_BASE);
  if(targetIndex < 0)
    return(0);
  else
    return(sparse->value[targetIndex]);
}


REAL addtoItem(sparseVector *sparse, int targetIndex, REAL value)
{
  int idx;

  if(targetIndex > 0) 
    idx = findIndex(targetIndex, sparse->index, sparse->count, BLAS_BASE);
  else {
    idx = -targetIndex;
    if(idx > sparse->count)
      /* Index error; ignore item */
      return(0.0);
  }

  if(idx <=0 )
    value = putItem(sparse, targetIndex, value);
  else {
    value += sparse->value[idx];
    putItem(sparse, -idx, value);
  }
  return(value);
}


REAL putItem(sparseVector *sparse, int targetIndex, REAL value)
{
  REAL last = 0.0;
  int  posIndex;

  if(targetIndex < 0) {
    posIndex = -targetIndex;
    if(posIndex > sparse->count)
      return(last);
    targetIndex = sparse->index[posIndex];
  }
  else
    posIndex = findIndex(targetIndex, sparse->index, sparse->count, BLAS_BASE);

  if(fabs(value) < MACHINEPREC)
    value = 0;

  if(targetIndex == sparse->index[0]) 
    sparse->value[0] = value;

  if(posIndex < 0) {
    if(value != 0) {
      if(sparse->count == sparse->size)
        resizeVector(sparse, sparse->size + RESIZEDELTA);
      posIndex = -posIndex;
      sparse->count++;
      if(posIndex < sparse->count) 
        moveVector(sparse, posIndex+1, posIndex, sparse->count-posIndex);
      sparse->value[posIndex] = value;
      sparse->index[posIndex] = targetIndex;
    }
  }
  else {
    if(value == 0) {
      last = sparse->value[posIndex];
      if(sparse->count > posIndex) 
        moveVector(sparse, posIndex, posIndex+1, sparse->count-posIndex);
      sparse->count--;
    }
    else {
      sparse->value[posIndex] = value;
      sparse->index[posIndex] = targetIndex;
    }
  }

#ifdef DEBUG_SPARSELIB
  verifyVector(sparse);
#endif

  return(last);
}


void swapItems(sparseVector *sparse, int firstIndex, int secondIndex)
{
  int i,j,ki,kj;
  REAL hold;

  if(firstIndex == secondIndex)
    return;
  if(firstIndex > secondIndex) {
    i = firstIndex;
    firstIndex = secondIndex;
    secondIndex = i;
  }

  if(FALSE) {
    i = 1;
    ki = 0;
    while(i <= sparse->count && (ki = sparse->index[i])<firstIndex) i++;
    j = i;
    kj = 0;
    while(j <= sparse->count && (kj = sparse->index[j])<secondIndex) j++;
  }
  else {
    i = findIndex(firstIndex, sparse->index, sparse->count, BLAS_BASE);
    if(i < 0)
      i = -i;
    j = findIndex(secondIndex, sparse->index, sparse->count, BLAS_BASE);
    if(j < 0)
      j = -j;
  }

  if(i > sparse->count)
    ki = 0;
  else
    ki = sparse->index[i];
  if(j > sparse->count)
    kj = 0;
  else
    kj = sparse->index[j];

  if(ki == firstIndex && kj == secondIndex) {   /* Found both -> swap in place */
    hold = sparse->value[i];
    sparse->value[i] = sparse->value[j];
    sparse->value[j] = hold;

    if(sparse->index[0] == firstIndex)
      sparse->value[0] = sparse->value[i];
    else if(sparse->index[0] == secondIndex)
      sparse->value[0] = sparse->value[j];
  }
  else if(ki == firstIndex) {                   /* Found first, but not the second -> shift left */
    j--;
    if(i < j) {
      hold = sparse->value[i];
      moveVector(sparse, i, i+1, j-i);
      sparse->value[j] = hold;
    }
    sparse->index[j] = secondIndex;

    if(sparse->index[0] == firstIndex)
      sparse->value[0] = 0;
    else if(sparse->index[0] == secondIndex)
      sparse->value[0] = sparse->value[j];

  }
  else if(kj == secondIndex) {                  /* Found second, but not the first -> shift right */
    if(i < j) {
      hold = sparse->value[j];
      moveVector(sparse, i+1, i, j-i);
      sparse->value[i] = hold;
    }
    sparse->index[i] = firstIndex;

    if(sparse->index[0] == firstIndex)
      sparse->value[0] = sparse->value[i];
    else if(sparse->index[0] == secondIndex)
      sparse->value[0] = 0;
  }

#ifdef DEBUG_SPARSELIB
  verifyVector(sparse);
#endif

}


void clearVector(sparseVector *sparse, int indexStart, int indexEnd)
{
  int i;

  i = sparse->count;
  if(i==0) return;

  if(indexStart<=0)
    indexStart=sparse->index[1];
  if(indexEnd<=0)
    indexEnd=sparse->index[i];

  if(indexStart>indexEnd) return;

  if(sparse->index[0]>=indexStart && sparse->index[0]<=indexEnd) {
    sparse->value[0] = 0;
  }
  if(indexStart<=sparse->index[1] && indexEnd>=sparse->index[i]) 
    sparse->count = 0;
  else {
    while(i>0 && sparse->index[i]>indexEnd) i--;
    indexEnd = i;
    while(i>0 && sparse->index[i]>=indexStart) i--;
    indexStart = i+1;
    if(indexEnd>=indexStart) {
      i = sparse->count-indexEnd;
      moveVector(sparse, indexStart, indexEnd+1, i);
      sparse->count -= indexEnd-indexStart+1;
    }
  }

#ifdef DEBUG_SPARSELIB
  verifyVector(sparse);
#endif

}


int getVector(sparseVector *sparse, REAL *dense, int indexStart, int indexEnd, MYBOOL doClear)
{
  int i,k;

  i = 1;
  while(i<=sparse->count && sparse->index[i]<indexStart) i++;

  while(i<=sparse->count && (k=sparse->index[i])<=indexEnd) {
    while(indexStart<k) {
      dense[indexStart] = 0;
      indexStart++;
    }
    dense[indexStart] = sparse->value[i];
    indexStart++;
    i++;
  }

  while(indexStart<=indexEnd) {
    dense[indexStart] = 0;
    indexStart++;
  }

  k = sparse->count;
  if(doClear) {
    sparse->count = 0;
    sparse->value[0] = 0;
  }
  return(k);
}

void putVector(sparseVector *sparse, REAL *dense, int indexStart, int indexEnd)
{
  int i,n;

  n = sparse->count;
  if(indexStart<=0)
    indexStart=sparse->index[1];
  if(indexEnd<=0)
    indexEnd=sparse->index[n];

  if(n==0 || sparse->index[n]<indexStart) {
    i = sparse->index[0];
    if(i>=indexStart && i<=indexEnd)
      sparse->value[0] = 0;
    for(i = indexStart; i<=indexEnd; i++) {
      if(dense[i] == 0) continue;
      if(sparse->size == sparse->count)
        resizeVector(sparse, sparse->size + RESIZEDELTA);
      sparse->count++;
      sparse->value[sparse->count] = dense[i];
      sparse->index[sparse->count] = i;
      if(i == sparse->index[0]) 
        sparse->value[0] = dense[i];
    }
  }
  else {
    while(indexStart <= indexEnd) {
      putItem(sparse, indexStart, dense[indexStart]);
      indexStart++;
    }
  }

#ifdef DEBUG_SPARSELIB
  verifyVector(sparse);
#endif

}


void fillVector(sparseVector *sparse, int count, REAL value)
{
  int i;

  if(sparse->count > 0) 
    clearVector(sparse, 0, 0);
  for(i = 1; i<=count; i++)
    putItem(sparse, i, value);
}


REAL dotVector(sparseVector *sparse, REAL *dense, int indexStart, int indexEnd)
{
  int  i, n;
  long REAL sum;
  
  n = sparse->count;
  sum = 0;

  if(n > 0) {
    if(indexStart<=0)
      indexStart=sparse->index[1];
    if(indexEnd<=0)
      indexEnd=sparse->index[n];

    if(indexStart > 1) {
      i = findIndex(indexStart, sparse->index, sparse->count, BLAS_BASE);
      if(i < 0) {
        i = -i;
        if(i > n) 
          return(sum);
      }
    }
    else
      i = 1;

    /* CPU intensive loop; provide alternative evaluation models */
#if defined DOFASTMATH
    {
      /* Do fast pointer arithmetic */
      int  *indexptr;
      REAL *valueptr;
/*      for(i = 1, indexptr = sparse->index + 1;
          i <= n && (*indexptr) < indexStart; i++, indexptr++); */
      indexptr = sparse->index + i;
      for(valueptr = sparse->value + i;
          i <= n && (*indexptr) <= indexEnd;  i++, indexptr++, valueptr++) 
        sum += (*valueptr) * dense[(*indexptr)];
    }
#else
    {
      /* Do traditional indexed access */
      int k;
/*      i = 1; */
/*      while(i<=n && sparse->index[i]<indexStart) i++; */
      while(i<=n && (k = sparse->index[i])<=indexEnd) {
        sum += sparse->value[i] * dense[k];
        i++;
      }
    }
#endif    
  }

  return(sum);
}


void daxpyVector1(sparseVector *sparse, REAL scalar, REAL *dense, int indexStart, int indexEnd)
{
  int i, n;

  if(scalar == 0) return;

  n = sparse->count;
  if(indexStart<=0)
    indexStart=sparse->index[1];
  if(indexEnd<=0)
    indexEnd=sparse->index[n];

  /* CPU intensive loop; provide alternative evaluation models */
#if defined DOFASTMATH
  {
    /* Do fast pointer arithmetic */
    int    *indexptr;
    REAL *valueptr;
    for(i = 1, indexptr = sparse->index + 1;
        i <= n && (*indexptr) < indexStart; i++, indexptr++);
    for(valueptr = sparse->value + i;
        i <= n && (*indexptr) <= indexEnd;  i++, indexptr++, valueptr++) 
      dense[(*indexptr)] += (*valueptr) * scalar;
  }
#else
  {
    /* Do traditional indexed access */
    int k;
    for(i = 1; i<= n; i++) {
      k = sparse->index[i];
      if(k<indexStart) continue;
      if(k>indexEnd) break;
      dense[k] += sparse->value[i] * scalar;
    }
  }
#endif  
}
void daxpyVector2(REAL *dense, REAL scalar, sparseVector *sparse, int indexStart, int indexEnd)
{
  sparseVector *hold;

  hold = createVector(sparse->limit, sparse->count);
  putDiagonalIndex(hold, getDiagonalIndex(sparse));
  putVector(hold, dense, indexStart, indexEnd);
  daxpyVector3(hold, scalar, sparse, indexStart, indexEnd);
  freeVector(hold);
}
void daxpyVector3(sparseVector *sparse1, REAL scalar, sparseVector *sparse2, int indexStart, int indexEnd)
{
  int i1, i2, k, p1, p2, c1, c2;
  sparseVector *hold;

  if(sparse1->count == 0) return;

 /* Spool to start positions */
  i1 = 1;
  c1 = sparse1->count;
  while(i1 <= c1 && sparse1->index[i1] < indexStart) i1++;
  if(i1 <= c1)
    p1 = sparse1->index[i1];
  else
    p1 = indexEnd+1;

  i2 = 1;
  c2 = sparse2->count;
  while(i2 <= c2 && sparse2->index[i2] < indexStart) i2++;
  if(i2 <= c2)
    p2 = sparse2->index[i2];
  else
    p2 = indexEnd+1;

 /* Create a temporary vector */
  k = c1+c2;
  if(k > 0) {
    hold = createVector(MAX(sparse1->limit, sparse2->limit), k);
    putDiagonalIndex(hold, getDiagonalIndex(sparse2));
  }
  else
    hold = sparse2;

 /* Loop over all items in both vectors */
  while((i1 <= c1 && p1 <= indexEnd) || 
        (i2 <= c2 && p2 <= indexEnd)) {

    k = 0;

   /* Add/spool exclusive right-vector items */
    while(i2 <= c2 && p2 < p1) {
      if(hold != sparse2)
        putItem(hold, p2, sparse2->value[i2]);
      i2++;
      if(i2 <= c2)
        p2 = sparse2->index[i2];
      else
        p2 = indexEnd+1;
      k++;
    }
   /* Add equal-indexed items */
    while(i1 <= c1 && i2 <= c2 && p1 == p2) {
/*      if(hold != sparse2) */
        putItem(hold, p1, scalar*sparse1->value[i1]+sparse2->value[i2]);
/*      else
          addtoItem(sparse2, -i2, scalar*sparse1->value[i1]); */
      i1++;
      if(i1 <= c1)
        p1 = sparse1->index[i1];
      else
        p1 = indexEnd+1;
      i2++;
      if(i2 <= c2)
        p2 = sparse2->index[i2];
      else
        p2 = indexEnd+1;
      k++;
    }
   /* Add exclusive left-vector items */
    while(i1 <= c1 && p1 < p2) {
      putItem(hold, p1, scalar*sparse1->value[i1]);
/*      if(hold == sparse2) c2++; */
      i1++;
      if(i1 <= c1)
        p1 = sparse1->index[i1];
      else
        p1 = indexEnd+1;
      k++;
    }

    if(k == 0) break;
  }

/*  if(hold != sparse2) */
  {
    swapVector(hold, sparse2);
    freeVector(hold);
  }

#ifdef DEBUG_SPARSELIB
  verifyVector(sparse2);
#endif

}


void dswapVector1(sparseVector *sparse, REAL *dense, int indexStart, int indexEnd)
{
  int i, d, n;
  REAL *x;

  if(indexStart <= 0)
    indexStart = 1;
  n = lastIndex(sparse);
  if(indexEnd <= 0) 
    indexEnd = n;
  CALLOC(x, (MAX(indexEnd,n)+1));

  getVector(sparse, x, indexStart, n, FALSE);
  d = getDiagonalIndex(sparse);
  clearVector(sparse, indexStart, n);
  for(i = indexStart; i<=indexEnd; i++) {
    if(dense[i] != 0)
      putItem(sparse, i, dense[i]);
  }
  for(i = indexEnd+1; i<=n; i++) {
    if(x[i] != 0)
      putItem(sparse, i, x[i]);
  }
  MEMCOPY(&dense[indexStart], &x[indexStart], (indexEnd-indexStart+1));

#ifdef DEBUG_SPARSELIB
  verifyVector(sparse);
#endif
  
  FREE(x);
}
void dswapVector2(REAL *dense, sparseVector *sparse, int indexStart, int indexEnd)
{
  dswapVector1(sparse, dense, indexStart, indexEnd);
}


void dswapVector3(sparseVector *sparse1, sparseVector *sparse2, int indexStart, int indexEnd)
{

  REAL *dense1, *dense2;

  if(indexStart<=0)
    indexStart = 1;
  if(indexEnd<=0)
    indexEnd = MAX(lastIndex(sparse1), lastIndex(sparse2));

  if(indexStart <= firstIndex(sparse1) && indexStart <= firstIndex(sparse2) && 
     indexEnd >= lastIndex(sparse1) && indexEnd >= lastIndex(sparse2)) {
    swapVector(sparse1, sparse2);
  }
  else {

    CALLOC(dense1, (indexEnd+1));
    CALLOC(dense2, (indexEnd+1));
    getVector(sparse1, dense1, indexStart, indexEnd, TRUE);
    getVector(sparse2, dense2, indexStart, indexEnd, TRUE);
    clearVector(sparse1, indexStart, indexEnd);
    clearVector(sparse2, indexStart, indexEnd);
    putVector(sparse1, dense2, indexStart, indexEnd);
    putVector(sparse2, dense1, indexStart, indexEnd);
    FREE(dense1);
    FREE(dense2);
  }
}


int idamaxVector(sparseVector *sparse, int is, REAL *maxValue)
{
  int    i, n, imax;
  REAL xmax;

  n = sparse->count;
  imax = 1;
  if(n == 0)
    xmax = 0;
  else {
    xmax = fabs(sparse->value[imax]);

    /* CPU intensive loop; provide alternative evaluation models */
#if defined DOFASTMATH
    {
      /* Do fast pointer arithmetic */
      int  *indexptr;
      REAL *valueptr;
      for(i = 1, indexptr = sparse->index + 1;
          i <= n && (*indexptr) <= is; i++, indexptr++);
      for(valueptr = sparse->value + i;
          i <= n; i++, indexptr++, valueptr++) {
	      if((*valueptr)>xmax) {
		      xmax = (*valueptr);
		      imax = (*indexptr);
        }
      }
    }
#else
    {
      REAL xtest;
      /* Do traditional indexed access */
      i = 1;
      while(i <= n && sparse->index[i] <= is) i++;
      for(; i<=n; i++) {
        xtest = fabs(sparse->value[i]);
	      if(xtest>xmax) {
		      xmax = xtest;
		      imax = sparse->index[i];
        }
      }
    }
#endif    
  }
  if(maxValue != NULL)
    (*maxValue) = sparse->index[imax];
  return(imax);
}


void printVector(int n, sparseVector *sparse, int modulo )
{
  int i,j,k;

  if(sparse == NULL) return;

  if (modulo <= 0) modulo = 5;
  for (i = 1, j = 1; j<=n; i++, j++) {
    if(i<=sparse->count)
      k = sparse->index[i];
    else
      k = n+1;
    while (j < k) {
      if(mod(j, modulo) == 1) 
        printf("\n%2d:%12g", j, 0.0);
      else
        printf(" %2d:%12g", j, 0.0);
      j++;
    }
    if(k<=n) {
      if(mod(j, modulo) == 1) 
        printf("\n%2d:%12g", k, sparse->value[i]);
      else
        printf(" %2d:%12g", k, sparse->value[i]);
    }
  }
  if(mod(j, modulo) != 0) printf("\n");
}