#include <stdlib.h>
#include "petscmat.h"
#include "/home/hpcin/soft/petsc-3.0.0-p8/src/mat/impls/dense/mpi/mpidense.c"

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **argv)
{

  Mat            A,B,C;
  PetscScalar    *a,*b,*c;
  PetscErrorCode ierr;
  PetscInt       rows,columns, i,j;
  PetscScalar    matrixElementValue; 
  PetscMPIInt    rank,size;

  PetscInitialize(&argc,&argv,0,0);
  MPI_Comm_size(PETSC_COMM_WORLD,&size);
  MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

  FILE *archive;

  archive = fopen("matrix.txt", "r"); //open archive
 
  fscanf (archive, "%d", &rows); //reading the sizes of matrix
  fscanf (archive, "%d", &columns);
 
  printf("rows: %d\n",rows);
  printf("columns: %d\n",columns); 

  PetscPrintf(PETSC_COMM_WORLD,"Number of processors = %d, rank = %d\n",size,rank);
 
  // Allocate array 
  ierr = PetscMalloc(rows*columns*sizeof(PetscScalar),&a);CHKERRQ(ierr);
  ierr = PetscMalloc(rows*columns*sizeof(PetscScalar),&b);CHKERRQ(ierr);
  ierr = PetscMalloc(rows*columns*sizeof(PetscScalar),&c);CHKERRQ(ierr);
 
  //Setting values
  for (i=0; i<rows; i++) {
    for (j=0; j<columns; j++) {
      fscanf(archive, "%lf", &matrixElementValue);
      a[i+j*rows] = matrixElementValue; 
      b[i+j*rows] = a[i+j*rows]; 
      c[i+j*rows] = 0;
    }
  }

  fclose(archive); //closing

  ierr = MatCreateMPIDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, rows, columns, a,&A);CHKERRQ(ierr);     
  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);   

  ierr = MatCreateMPIDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, rows, columns, b,&B);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  
  ierr = MatCreateMPIDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, rows, columns, c,&C);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(C,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(C,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);

  ierr = MatMatMultNumeric_MPIDense_MPIDense(A,B,C);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  

  printf("MatMatMult C:\n"); MatView(C,0);
 

  ierr = PetscFree(a);CHKERRQ(ierr);
  ierr = PetscFree(b);CHKERRQ(ierr);
  ierr = PetscFree(c);CHKERRQ(ierr);

  ierr = MatDestroy(A);CHKERRQ(ierr);
  ierr = MatDestroy(B);CHKERRQ(ierr);
  ierr = MatDestroy(C);CHKERRQ(ierr);

 
  ierr = PetscFinalize();CHKERRQ(ierr);
  return 0;
}