/*
 * Some basic functions for mtx reading and formating
 * 
 * Author: Petros Anastasiadis(panastas@cslab.ece.ntua.gr) 
 */
#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#include <errno.h>
#include "matrix_op.h"
#include "util.h"
#include "input.h"

int mtx_read(int ** I, int ** cooCol, double ** cooVal, int * n, int * m, int * n_z, char * name)
{
	
	char c;
	char *type, *format, *var_type, *symmetry, *string=NULL;
	FILE *fp ;
	size_t len=0;

	if ((fp=fopen(name, "r"))==NULL) error("Invalid File");
	else if ((strstr(name, "mtx"))==NULL) error("Invalid File Type (*.mtx)");

	getline(&string, &len, fp);
	strtok(string," ");
	type = strtok(NULL," ");
	format = strtok(NULL," ");
	var_type = strtok(NULL," ");
	symmetry = strtok(NULL,"\n");
	//printf("type=%s, format=%s, var_type=%s, ", type, format, var_type);
	if (strcmp(type,"matrix")){
		printf("type=%s unsupported...terminating\n\n\n\n\n\n\n\n\n\n\n\n", type);
		exit(1);
	}
	if (strcmp(format,"coordinate")){
		printf("format=%s unsupported...terminating\n\n\n\n\n\n\n\n\n\n\n\n", format);
		exit(1);
	}
	if (strcmp(var_type,"integer") && strcmp(var_type,"real")){
		printf("Var_type=%s unsupported...terminating\n\n\n\n\n\n\n\n\n\n\n\n", var_type);
		exit(1);
	}
	while((c=getc(fp))=='%') while( (c=getc(fp))!='\n') ; 
	ungetc(c, fp);
	int k, lines = 0, sym_k=0;
	fscanf(fp,"%d %d %d", n, m, &lines);
	//printf("n=%d, m=%d, lines=%d, ", *n, *m, lines);
	
	*n_z = 0;
	if (!strcmp(symmetry,"symmetric")){
		get_nz_symmetric(n_z, name);
		//printf("symmetry=symmetric\n");
	}
	else if (!strcmp(symmetry,"general")) {
		*n_z=lines;
		//printf("symmetry=general\n");
	}
	else {
		printf("Invalid symmetry value:%s\n", symmetry); 
		return 0; 
	}
	//printf("n_z=%d\n", *n_z);
	*I = (int *) malloc(*n_z * sizeof(int));
	*cooCol = (int *) malloc(*n_z * sizeof(int));
	*cooVal = (double *) malloc(*n_z * sizeof(double));
	double dum;
	if ( !*I || !*cooCol || !*cooVal ) return 0;
	if (!strcmp(symmetry,"symmetric")){
		for (k = 0; k < lines; k++) {
			fscanf(fp,"%d %d %lf", &((*I)[sym_k]), &((*cooCol)[sym_k]), &dum);
			(*cooVal)[sym_k]=(double) dum;
			(*I)[sym_k]--;
			(*cooCol)[sym_k]--;
			sym_k++;
			if ((*I)[sym_k-1] != (*cooCol)[sym_k-1]) {
				(*I)[sym_k] = (*cooCol)[sym_k-1];
				(*cooCol)[sym_k] = (*I)[sym_k-1];
				(*cooVal)[sym_k] = (*cooVal)[sym_k-1];
				sym_k++;
			}
		}
		if (sym_k!=*n_z){
			printf("Error in symmetric read: sym_k=%d n_z=%d\n", sym_k, *n_z);
			return 0;
		}
	}
	else if (!strcmp(symmetry,"general")) for (k = 0; k < lines; k++){
		fscanf(fp,"%d %d %lf", &((*I)[k]), &((*cooCol)[k]), &dum);
		(*cooVal)[k] = (double) dum;
		(*I)[k]--;
		(*cooCol)[k]--;
	}
	quickSort( *I, *cooCol, *cooVal, 0, *n_z-1);
	fclose(fp);
	return 1;
}

void get_nz_symmetric( int * n_z, char* name)
{
	char c;
	FILE *fp ;
	if ((fp=fopen(name, "r"))==NULL){
		printf("Problem in symmetric read pass\n");
		exit(1);
	}

	while((c=getc(fp))=='%') while( (c=getc(fp))!='\n') ; 
	ungetc(c, fp);
	int k, i, j, n, m, lines;
	double x;
	fscanf(fp,"%d %d %d", &n, &m, &lines);
	for (k = 0; k < lines; k++){
		fscanf(fp,"%d %d %lf", &i, &j, &x);
		(*n_z)++;
		if(i!=j) (*n_z)++;
	}
}

void csr_transform(double ** A, int n, int m, int n_z, double  *csrValA, int *csrRowPtrA, int *csrColIndA)
{
	int i,j,k=0;
	for (i = 0; i < n; i++){
		csrRowPtrA[i]=k;
		for (j = 0; j < m; j++){
			if (A[i][j]!=0.0){
				csrValA[k]=A[i][j];
				csrColIndA[k]= j;
				k++;
			}
		}
	}
	csrRowPtrA[i]=k;
	if (k!=n_z) printf("Error at non zeroes: %d\n", k-n_z);
	return;
}			

void quickSort( int *a, int * b, double * c, int l, int r)
{
	int j;
	if( l < r ) 
	{	// divide and conquer		
		j = partition( a, b, c, l, r);
		quickSort( a, b, c, l, j-1);
		quickSort( a, b, c, j+1, r);
	}
}

int partition( int *a, int * b, double * c, int l, int r) 
{
	int pivot, i, j, t;
	double t1;
	pivot = a[l];
	i = l; j = r+1;
		
	while(1)
	{
		do ++i; while( a[i] <= pivot && i <= r );
   		do --j; while( a[j] > pivot );
   		if( i >= j ) break;
   		t = a[i]; a[i] = a[j]; a[j] = t;
		t = b[i]; b[i] = b[j]; b[j] = t;
		t1 = c[i]; c[i] = c[j]; c[j] = t1;
   	}
   	t = a[l]; a[l] = a[j]; a[j] = t;
	t = b[l]; b[l] = b[j]; b[j] = t;
	t1 = c[l]; c[l] = c[j]; c[j] = t1;
   	return j;
}