/*************************************************************************/
/*                                                                       */
/*  Copyright (c) 1994 Stanford University                               */
/*                                                                       */
/*  All rights reserved.                                                 */
/*                                                                       */
/*  Permission is given to use, copy, and modify this software for any   */
/*  non-commercial purpose as long as this copyright notice is not       */
/*  removed.  All other uses, including redistribution in whole or in    */
/*  part, are forbidden without prior written permission.                */
/*                                                                       */
/*  This software is provided with absolutely no warranty and no         */
/*  support.                                                             */
/*                                                                       */
/*************************************************************************/

EXTERN_ENV

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "matrix.h"

#define Error(m) { printf(m); exit(0); }
#define AddMember(set, new) { long s, n; s = set; n = new; link[n] = link[s]; link[s] = n; }

long maxm;

SMatrix NewMatrix(long n, long m, long nz)
{
  SMatrix M;

  M.n = n; M.m = m;
  M.col = (long *) MyMalloc((n+1)*sizeof(long), DISTRIBUTED);
  M.startrow = (long *) MyMalloc((n+1)*sizeof(long), DISTRIBUTED);
  M.row = (long *) MyMalloc((m+n)*sizeof(long), DISTRIBUTED);
  if (nz) {
	M.nz = (double *) MyMalloc((m+n)*sizeof(double), DISTRIBUTED);
  } else {
	M.nz = NULL;
  }

  if (!M.col || !M.row || (nz && !M.nz)) {
    printf("NewMatrix %ld %ld: Out of memory\n", n, m);
    exit(0);
  }

  return(M);
}


void FreeMatrix(SMatrix M)
{
  MyFree(M.col);
  MyFree(M.startrow);
  MyFree(M.row);
  if (M.nz)
    MyFree((void *)M.nz);
}


double *NewVector(long n)
{
  double *v;

  v = (double *) MyMalloc(n*sizeof(double), DISTRIBUTED);

  if (!v && n) {
    printf("Out of memory: NewVector(%ld)\n", n);
    exit(0);
  }

  return(v);
}


double Value(long i, long j)
{
  if (i == j) {
	return((double) maxm+0.1);
  } else {
	return(-1.0);
  }
}


SMatrix ReadSparse(char *name, char *probName)
{
	FILE *fp;
	long n, m, i, j;
	long n_rows, tmp;
	long numer_lines;
	long colnum, colsize, rownum, rowsize;
	char buf[100], type[4];
	SMatrix M, F;

	if (!name || name[0] == 0) {
		fp = stdin;
	} else {
		fp = fopen(name, "r");
	}

	if (!fp) {
		Error("Error opening file\n");
	}

	fscanf(fp, "%72c", buf);

	fscanf(fp, "%8c", probName);
	probName[8] = 0;
	DumpLine(fp);

	for (i=0; i<5; i++) {
	  fscanf(fp, "%14c", buf);
	  sscanf(buf, "%ld", &tmp);
	  if (i == 3)
	    numer_lines = tmp;
	}
	DumpLine(fp);

	fscanf(fp, "%3c", type);
	type[3] = 0;
	if (!(type[0] != 'C' && type[1] == 'S' && type[2] == 'A')) {
	  fprintf(stderr, "Wrong type: %s\n", type);
	  exit(0);
	}

	fscanf(fp, "%11c", buf); /* pad */

	fscanf(fp, "%14c", buf); sscanf(buf, "%ld", &n_rows);

	fscanf(fp, "%14c", buf); sscanf(buf, "%ld", &n);

	fscanf(fp, "%14c", buf); sscanf(buf, "%ld", &m);

	fscanf(fp, "%14c", buf); sscanf(buf, "%ld", &tmp);
	if (tmp != 0)
	  printf("This is not an assembled matrix!\n");
	if (n_rows != n)
	  printf("Matrix is not symmetric\n");
	DumpLine(fp);

	fscanf(fp, "%16c", buf);
	ParseIntFormat(buf, &colnum, &colsize);
	fscanf(fp, "%16c", buf);
	ParseIntFormat(buf, &rownum, &rowsize);
	fscanf(fp, "%20c", buf);
	fscanf(fp, "%20c", buf);
		
	DumpLine(fp); /* format statement */

	M = NewMatrix(n, m, 0);

	ReadVector(fp, n+1, M.col, colnum, colsize);

	ReadVector(fp, m, M.row, rownum, rowsize);

	for (i=0; i<numer_lines; i++) /* dump numeric values */
	  DumpLine(fp);

	for (i=0; i<n; i++)
		ISort(M, i);

	for (i=0; i<=n; i++)
	  M.startrow[i] = M.col[i];

	fclose(fp);

	F = LowerToFull(M);

	maxm = 0;
	for (i=0; i<n; i++)
	  if (F.col[i+1]-F.col[i] > maxm)
	    maxm = F.col[i+1]-F.col[i];

	if (F.nz) {
	  for (j=0; j<n; j++)
	    for (i=F.col[j]; i<F.col[j+1]; i++)
	      F.nz[i] = Value(F.row[i], j);
	}

	FreeMatrix(M);

	return(F);
}

void DumpLine(FILE *fp)
{
	long c;

	while ((c = fgetc(fp)) != '\n')
		;
}

void ParseIntFormat(char *buf, long *num, long *size)
{
  char *tmp;

  tmp = buf;

  while (*tmp++ != '(')
    ;
  sscanf(tmp, "%ld", num);

  while (*tmp++ != 'I')
    ;
  sscanf(tmp, "%ld", size);
}


void ReadVector(FILE *fp, long n, long *where, long perline, long persize)
{
  long i, j, item;
  char tmp, buf[100];

  i = 0;
  while (i < n) {
    fgets(buf, 100, fp);    /* read a line at a time */
    for (j=0; j<perline && i<n; j++) {
      tmp = buf[(j+1)*persize]; buf[(j+1)*persize] = 0;  /* null terminate */
      item = atoi(&buf[j*persize]);
      buf[(j+1)*persize] = tmp;
      where[i++] = item-1;
    }
  }
}

SMatrix LowerToFull(SMatrix L)
{
  SMatrix M;
  long *link, *first;
  long i, j, nextj, ind = 0;

  link = (long *) malloc(L.n*sizeof(long));
  first = (long *) malloc(L.n*sizeof(long));

  for (i=0; i<L.n; i++) {
    link[i] = first[i] = -1;
  }

  M = NewMatrix(L.n, 2*(L.m-L.n)+L.n, 0);

  for (i=0; i<L.n; i++) {
    M.col[i] = ind;

    for (j=L.col[i]; j<L.col[i+1]; j++) {
      if (L.row[j] >= i) {
	M.row[ind++] = L.row[j];
      }
    }

    j = link[i];
    while (j != -1) {
      nextj = link[j];
      M.row[ind++] = j;
      first[j]++;
      if (first[j] < L.col[j+1]) {
	AddMember(L.row[first[j]], j);
      }
      j = nextj;
    }

    first[i] = L.col[i];
    if (L.row[first[i]] == i) {
      first[i]++;
    } else {
      fprintf(stderr, "Missing diagonal: %ld: ", i);
      for (j=L.col[i]; j<L.col[i+1]; j++) {
	fprintf(stderr, "%ld ", L.row[j]);
      }
      fprintf(stderr, "\n");
    }

    if (first[i] < L.col[i+1]) {
      AddMember(L.row[first[i]], i);
    }
  }

  M.col[M.n] = ind;

  for (i=0; i<=L.n; i++) {
    M.startrow[i] = M.col[i];
  }

  if (ind != M.m) {
    printf("Lost some\n");
  }

  free(link); free(first);

  return(M);
}


void ISort(SMatrix M, long k)
{
	long hi, lo;
	long i, j, tmp;

	hi = M.col[k+1];
	lo = M.col[k];

	for (i=lo; i<hi; i++) {
		tmp = M.row[i];
		j = i;
		while (M.row[j-1] > tmp && j > lo) {
			M.row[j] = M.row[j-1];
			j--;
			}
		M.row[j] = tmp;
		}

}