/////////////////////////////////////////////////////////////////////////////////
////
////  CRS sparse matrices manipulation routines
////  Copyright (C) 2004-2008 Manolis Lourakis (lourakis at ics forth gr)
////  Institute of Computer Science, Foundation for Research & Technology - Hellas
////  Heraklion, Crete, Greece.
////
////  This program is free software; you can redistribute it and/or modify
////  it under the terms of the GNU General Public License as published by
////  the Free Software Foundation; either version 2 of the License, or
////  (at your option) any later version.
////
////  This program is distributed in the hope that it will be useful,
////  but WITHOUT ANY WARRANTY; without even the implied warranty of
////  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
////  GNU General Public License for more details.
////
///////////////////////////////////////////////////////////////////////////////////

#include <stdio.h>
#include <stdlib.h>

#include "sba.h"

/* allocate a sparse CRS matrix */
void sba_crsm_alloc(struct sba_crsm *sm, int nr, int nc, int nnz) {
	int msz;

	sm->nr = nr;
	sm->nc = nc;
	sm->nnz = nnz;
	msz = 2 * nnz + nr + 1;
	sm->val = (int *)malloc(msz * sizeof(int)); /* required memory is allocated in a single step */
	if (!sm->val) {
		fprintf(stderr, "SBA: memory allocation request failed in sba_crsm_alloc() [nr=%d, nc=%d, nnz=%d]\n", nr, nc,
				nnz);
		exit(1);
	}
	sm->colidx = sm->val + nnz;
	sm->rowptr = sm->colidx + nnz;
}

/* free a sparse CRS matrix */
void sba_crsm_free(struct sba_crsm *sm) {
	sm->nr = sm->nc = sm->nnz = -1;
	free(sm->val);
	sm->val = sm->colidx = sm->rowptr = NULL;
}

/* returns the index of the (i, j) element. No bounds checking! */
int sba_crsm_elmidx(struct sba_crsm *sm, int i, int j) {
	register int low, high, mid, diff;

	low = sm->rowptr[i];
	high = sm->rowptr[i + 1] - 1;

	/* binary search for finding the element at column j */
	while (low <= high) {
		/* following early termination test seems to actually slow down the search */
		// if(j<sm->colidx[low] || j>sm->colidx[high]) return -1; /* not found */

		/* mid=low+((high-low)>>1) ensures no index overflows */
		mid = (low + high) >> 1; //(low+high)/2;
		diff = j - sm->colidx[mid];
		if (diff < 0)
			high = mid - 1;
		else if (diff > 0)
			low = mid + 1;
		else
			return mid;
	}

	return -1; /* not found */
}

/* similarly to sba_crsm_elmidx() above, returns the index of the (i, j) element using the
 * fact that the index of element (i, jp) was previously found to be jpidx. This can be
 * slightly faster than sba_crsm_elmidx(). No bounds checking!
 */
int sba_crsm_elmidxp(struct sba_crsm *sm, int i, int j, int jp, int jpidx) {
	register int low, high, mid, diff;

	diff = j - jp;
	if (diff > 0) {
		low = jpidx + 1;
		high = sm->rowptr[i + 1] - 1;
	} else if (diff == 0)
		return jpidx;
	else { /* diff<0 */
		low = sm->rowptr[i];
		high = jpidx - 1;
	}

	/* binary search for finding the element at column j */
	while (low <= high) {
		/* following early termination test seems to actually slow down the search */
		// if(j<sm->colidx[low] || j>sm->colidx[high]) return -1; /* not found */

		/* mid=low+((high-low)>>1) ensures no index overflows */
		mid = (low + high) >> 1; //(low+high)/2;
		diff = j - sm->colidx[mid];
		if (diff < 0)
			high = mid - 1;
		else if (diff > 0)
			low = mid + 1;
		else
			return mid;
	}

	return -1; /* not found */
}

/* returns the number of nonzero elements in row i and
 * fills up the vidxs and jidxs arrays with the val and column
 * indexes of the elements found, respectively.
 * vidxs and jidxs are assumed preallocated and of max. size sm->nc
 */
int sba_crsm_row_elmidxs(struct sba_crsm *sm, int i, int *vidxs, int *jidxs) {
	register int j, k;

	for (j = sm->rowptr[i], k = 0; j < sm->rowptr[i + 1]; ++j, ++k) {
		vidxs[k] = j;
		jidxs[k] = sm->colidx[j];
	}

	return k;
}

/* returns the number of nonzero elements in col j and
 * fills up the vidxs and iidxs arrays with the val and row
 * indexes of the elements found, respectively.
 * vidxs and iidxs are assumed preallocated and of max. size sm->nr
 */
int sba_crsm_col_elmidxs(struct sba_crsm *sm, int j, int *vidxs, int *iidxs) {
	register int *nextrowptr = sm->rowptr + 1;
	register int i, l;
	register int low, high, mid, diff;

	for (i = l = 0; i < sm->nr; ++i) {
		low = sm->rowptr[i];
		high = nextrowptr[i] - 1;

		/* binary search attempting to find an element at column j */
		while (low <= high) {
			// if(j<sm->colidx[low] || j>sm->colidx[high]) break; /* not found */

			mid = (low + high) >> 1; //(low+high)/2;
			diff = j - sm->colidx[mid];
			if (diff < 0)
				high = mid - 1;
			else if (diff > 0)
				low = mid + 1;
			else { /* found */
				vidxs[l] = mid;
				iidxs[l++] = i;
				break;
			}
		}
	}

	return l;
}

/* a more straighforward (but slower) implementation of the above function */
/***
int sba_crsm_col_elmidxs(struct sba_crsm *sm, int j, int *vidxs, int *iidxs)
{
register int i, k, l;

  for(i=l=0; i<sm->nr; ++i)
	for(k=sm->rowptr[i]; k<sm->rowptr[i+1]; ++k)
	  if(sm->colidx[k]==j){
		vidxs[l]=k;
		iidxs[l++]=i;
	  }

  return l;
}
***/

#if 0
/* returns 1 if there exists a row i having columns j and k,
 * i.e. a row i s.t. elements (i, j) and (i, k) are nonzero;
 * 0 otherwise
 */ 
int sba_crsm_common_row(struct sba_crsm *sm, int j, int k)
{
register int i, low, high, mid, diff;

  if(j==k) return 1;

  for(i=0; i<sm->nr; ++i){
    low=sm->rowptr[i];
    high=sm->rowptr[i+1]-1;
    if(j<sm->colidx[low] || j>sm->colidx[high] || /* j not found */
       k<sm->colidx[low] || k>sm->colidx[high])   /* k not found */
      continue;

    /* binary search for finding the element at column j */
    while(low<=high){
      mid=(low+high)>>1; //(low+high)/2;
      diff=j-sm->colidx[mid];
      if(diff<0)
        high=mid-1;
      else if(diff>0)
        low=mid+1;
      else
        goto jfound;
    }

    continue; /* j not found */

jfound:
    if(j>k){
      low=sm->rowptr[i];
      high=mid-1;
    }
    else{
      low=mid+1;
      high=sm->rowptr[i+1]-1;
    }

    if(k<sm->colidx[low] || k>sm->colidx[high]) continue; /* k not found */

    /* binary search for finding the element at column k */
    while(low<=high){
      mid=(low+high)>>1; //(low+high)/2;
      diff=k-sm->colidx[mid];
      if(diff<0)
        high=mid-1;
      else if(diff>0)
        low=mid+1;
      else /* found */
        return 1;
    }
  }

  return 0;
}
#endif

#if 0

/* sample code using the above routines */

main()
{
int mat[7][6]={
    {10, 0, 0, 0, -2, 0},
    {3,  9, 0, 0,  0, 3},
    {0,  7, 8, 7,  0, 0},
    {3,  0, 8, 7,  5, 0},
    {0,  8, 0, 9,  9, 13},
    {0,  4, 0, 0,  2, -1},
    {3,  7, 0, 9,  2, 0}
};

struct sba_crsm sm;
int i, j, k, l;
int vidxs[7], /* max(6, 7) */
    jidxs[6], iidxs[7];


  sba_crsm_build(&sm, mat[0], 7, 6);
  sba_crsm_print(&sm, stdout);

  for(i=0; i<7; ++i){
    for(j=0; j<6; ++j)
      printf("%3d ", ((k=sba_crsm_elmidx(&sm, i, j))!=-1)? sm.val[k] : 0);
    printf("\n");
  }

  for(i=0; i<7; ++i){
    k=sba_crsm_row_elmidxs(&sm, i, vidxs, jidxs);
    printf("row %d\n", i);
    for(l=0; l<k; ++l){
      j=jidxs[l];
      printf("%d %d  ", j, sm.val[vidxs[l]]); 
    }
    printf("\n");
  }

  for(j=0; j<6; ++j){
    k=sba_crsm_col_elmidxs(&sm, j, vidxs, iidxs);
    printf("col %d\n", j);
    for(l=0; l<k; ++l){
      i=iidxs[l];
      printf("%d %d  ", i, sm.val[vidxs[l]]); 
    }
    printf("\n");
  }

  sba_crsm_free(&sm);
}
#endif