static char help[] = "Test PtAP \n\
Reads PETSc matrix A and P, then comput Pt*A*P \n\
Input parameters include\n\
  -fA <input_file> -fP <input_file>: second files to load (projection) \n\n";

#include "petscmat.h"

#undef WRITEFILE
#undef __FUNCT__
#define __FUNCT__ "main"
PetscInt main(PetscInt argc,char **args)
{
  Mat            A,P,C;    
  PetscViewer    fd;               
  char           file[2][PETSC_MAX_PATH_LEN];     
  PetscTruth     flg; 
  PetscErrorCode ierr;
  PetscReal      fill=2.0;

  PetscInitialize(&argc,&args,(char *)0,help);
#ifdef WRITEFILE
 {
    PetscViewer    viewer;
    ierr = PetscPrintf(PETSC_COMM_WORLD,"writing matrix A in binary to A.dat ...\n");CHKERRQ(ierr);
    ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"A.dat",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
    ierr = MatView(A,viewer);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr);

    ierr = PetscPrintf(PETSC_COMM_WORLD,"writing matrix P in binary to P.dat ...\n");CHKERRQ(ierr);
    ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"P.dat",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
    ierr = MatView(P,viewer);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr);
 }
#endif

  /* read the two matrices, A (square) and P (projection) */
  ierr = PetscOptionsGetString(PETSC_NULL,"-fA",file[0],PETSC_MAX_PATH_LEN-1,&flg);CHKERRQ(ierr);
  if (!flg) SETERRQ(PETSC_ERR_USER,"Must indicate binary file with the -fA options");
  ierr = PetscOptionsGetString(PETSC_NULL,"-fP",file[1],PETSC_MAX_PATH_LEN-1,&flg);CHKERRQ(ierr);
  if (!flg) SETERRQ(PETSC_ERR_USER,"Must indicate binary file with the -fP options");
   
  /* Load matrices */
  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&fd);CHKERRQ(ierr);
  ierr  = MatLoad(fd,MATAIJ,&A);CHKERRQ(ierr);
  ierr = PetscViewerDestroy(fd);CHKERRQ(ierr); 
  //ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr);
    
  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[1],FILE_MODE_READ,&fd);CHKERRQ(ierr);
  ierr = MatLoad(fd,MATAIJ,&P);CHKERRQ(ierr);
  ierr = PetscViewerDestroy(fd);CHKERRQ(ierr);
   
  ierr = MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);CHKERRQ(ierr);

  ierr = MatDestroy(C);CHKERRQ(ierr);
  ierr = MatDestroy(P);CHKERRQ(ierr);
  ierr = MatDestroy(A);CHKERRQ(ierr);
  ierr = PetscFinalize();CHKERRQ(ierr);
  return 0;
}