static char help[] = "Tests sequential and parallel MatMatMult() and MatAXPY(...,SUBSET_NONZERO_PATTERN) \n\
Input arguments are:\n\
  -f <input_file>  : file to load\n\n";
/* e.g., mpiexec -n 3 ./ex113 -f <file> */

#include "petscmat.h"
#include<iostream>

void save_sparse_matrix(Mat& mat, const char *fname)
{
int ierr=0;
PetscViewer petsc_viewer;
ierr = PetscViewerCreate (PETSC_COMM_WORLD, &petsc_viewer);
ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,fname,&petsc_viewer);
ierr = PetscViewerSetFormat (petsc_viewer, PETSC_VIEWER_ASCII_MATLAB);
ierr = MatView (mat, petsc_viewer);
ierr = PetscViewerDestroy (petsc_viewer);
}




#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **args)
{
  Mat            A,A1,A2,Mtmp,dstMat;
  PetscViewer    viewer;
  PetscErrorCode ierr;
  PetscReal      fill=4.0;
  char           file[128]; 
  PetscTruth     flg;

  PetscInitialize(&argc,&args,(char *)0,help);
  
  /*  Load the matrix A */
  //ierr = PetscOptionsGetString(PETSC_NULL,"-f",file,127,&flg);CHKERRQ(ierr);
  //if (!flg) SETERRQ(1,"Must indicate a file name for matrix A with the -f option.");
  //sparse 
  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"F_bin.txt",FILE_MODE_READ,&viewer);CHKERRQ(ierr);
  ierr = MatLoad(viewer,MATMPIAIJ,&A1);CHKERRQ(ierr); 
  ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr);             

  //dense
  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"tempA_bin.txt",FILE_MODE_READ,&viewer);CHKERRQ(ierr);
  ierr = MatLoad(viewer,MATMPIDENSE,&A2);CHKERRQ(ierr); 
  ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr);             

  //ierr = MatDuplicate(A,MAT_COPY_VALUES,&A1);CHKERRQ(ierr);
  //ierr = MatDuplicate(A,MAT_COPY_VALUES,&A2);CHKERRQ(ierr);

  ierr = MatMatMult(A1,A2,MAT_INITIAL_MATRIX,fill,&A);CHKERRQ(ierr);
  //ierr = MatMatMult(A,Mtmp,MAT_INITIAL_MATRIX,fill,&dstMat);CHKERRQ(ierr); 
  //ierr = MatDestroy(Mtmp);CHKERRQ(ierr);
  std::cout<<"Finishing MatMatMult"<<std::endl;

 
  //char *fname="A.txt";
  //save_sparse_matrix(A,fname); 

  MatTranspose(A,MAT_REUSE_MATRIX,&A);

  std::cout<<"Finishing MatTranspose"<<std::endl;

  Mat subA;
  int row_g,col_g;
  MatGetSize(A,&row_g,&col_g);
  IS isrow,iscol;
  int rstart,rend; 
  ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
  ierr = ISCreateStride(PETSC_COMM_WORLD,(rend-rstart)/2,rstart,2,&isrow);CHKERRQ(ierr);
  /* list ALL the columns we want */
  ierr = ISCreateStride(PETSC_COMM_WORLD,col_g/2,3,3,&iscol);CHKERRQ(ierr);
  ierr = MatGetSubMatrix(A,isrow,iscol,PETSC_DECIDE,MAT_INITIAL_MATRIX,&subA);CHKERRQ(ierr);

  std::cout<<"Finishing MatGetSubMatrix"<<std::endl;

  //fname="subA.txt";
  //save_sparse_matrix(subA,fname);

  ierr = MatDestroy(A);CHKERRQ(ierr); 
  ierr = MatDestroy(A1);CHKERRQ(ierr);
  ierr = MatDestroy(A2);CHKERRQ(ierr); 
  MatDestroy(subA); 
  ierr = PetscFinalize();CHKERRQ(ierr);
  return 0;
}