#include "petscmat.h"
#define MAXLINE 1000

#undef __FUNCT__
#define __FUNCT__ "ReadMatFromFile"
/*
   ReadMatFromFile - Reads matrix data stored in ASCII i j value format and creates an AIJ matrix
   Input
.  File pointer *fp
   Output
.  SeqAIJ matrix

   Notes:
   Format of the matrix stored in the ASCII file
   Line 1 - Matrix Name (ignored)
   Line 2 - nrows ncols nz where nz = total number of nonzeros in the matrix
   Line 3 onwards matrix data in i j value format.

   If there are several matrices stored in the same file, then the file pointer should
   point to the matrix name
*/
PetscErrorCode ReadMatFromFile(FILE * fp,Mat * Mptr)
{
  Mat            A;
  char           line[MAXLINE]; 
  PetscErrorCode ierr;
  int            m,n,nz,*col=0,*row=0; /* these are fscaned so kept as int */
  PetscInt       i,col_i,row_i,*nnz,shift = 1;
  PetscScalar    val_i,*val = 0;
  
  PetscFunctionBegin;
  fgets(line,MAXLINE,fp); /* Skip first line  */
  /* Read rows,cols,number of nonzeros */
  fscanf(fp,"%d %d %d\n",&m,&n,&nz);
  
  ierr = PetscMalloc(m*sizeof(PetscInt),&nnz);CHKERRQ(ierr);
  ierr = PetscMemzero(nnz,m*sizeof(PetscInt));CHKERRQ(ierr);
  ierr = PetscMalloc(nz*sizeof(PetscReal),&val);CHKERRQ(ierr);
  ierr = PetscMalloc(nz*sizeof(PetscInt),&row);CHKERRQ(ierr);
  ierr = PetscMalloc(nz*sizeof(PetscInt),&col);CHKERRQ(ierr);
  for (i=0; i<nz; i++) {
     fscanf(fp,"%d %d %lf\n",row+i,col+i,(double*)val+i);
     row[i] -= shift; col[i] -= shift;  /* set index set starts at 0 */
      nnz[row[i]]++;
    }

  ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,m,n,0,nnz,&A);CHKERRQ(ierr);
  ierr = PetscFree(nnz);CHKERRQ(ierr);
  
  for (i=0; i<nz; i++) {
     row_i =(int)row[i]; col_i =(int)col[i]; val_i = val[i];
     ierr = MatSetValues(A,1,&row_i,1,&col_i,&val_i,INSERT_VALUES);CHKERRQ(ierr);
  }
  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = PetscFree(col);CHKERRQ(ierr);
  ierr = PetscFree(val);CHKERRQ(ierr);
  ierr = PetscFree(row);CHKERRQ(ierr);

  *Mptr = A; /* Point Mptr to A */
   
  PetscFunctionReturn(0);
}