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

  6: #include <petscmat.h>

  8: #undef WRITEFILE
 11: PetscInt main(PetscInt argc,char **args)
 12: {
 13:   Mat            A,P,C;
 14:   PetscViewer    fd;
 15:   char           file[2][PETSC_MAX_PATH_LEN];
 16:   PetscBool      flg;
 18:   PetscReal      fill=2.0;

 20:   PetscInitialize(&argc,&args,(char*)0,help);
 21: #if defined(WRITEFILE)
 22:   {
 23:     PetscViewer viewer;
 24:     PetscPrintf(PETSC_COMM_WORLD,"writing matrix A in binary to A.dat ...\n");
 25:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,"A.dat",FILE_MODE_WRITE,&viewer);
 26:     MatView(A,viewer);
 27:     PetscViewerDestroy(&viewer);

 29:     PetscPrintf(PETSC_COMM_WORLD,"writing matrix P in binary to P.dat ...\n");
 30:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,"P.dat",FILE_MODE_WRITE,&viewer);
 31:     MatView(P,viewer);
 32:     PetscViewerDestroy(&viewer);
 33:   }
 34: #endif

 36:   /* read the two matrices, A (square) and P (projection) */
 37:   PetscOptionsGetString(NULL,"-fA",file[0],PETSC_MAX_PATH_LEN,&flg);
 38:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Must indicate binary file with the -fA options");
 39:   PetscOptionsGetString(NULL,"-fP",file[1],PETSC_MAX_PATH_LEN,&flg);
 40:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Must indicate binary file with the -fP options");

 42:   /* Load matrices */
 43:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&fd);
 44:   MatCreate(PETSC_COMM_WORLD,&A);
 45:   MatLoad(A,fd);
 46:   PetscViewerDestroy(&fd);
 47:   /* MatGetSize(A,&m,&n); */

 49:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[1],FILE_MODE_READ,&fd);
 50:   MatCreate(PETSC_COMM_WORLD,&P);
 51:   MatLoad(P,fd);
 52:   PetscViewerDestroy(&fd);

 54:   MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);

 56:   MatDestroy(&C);
 57:   MatDestroy(&P);
 58:   MatDestroy(&A);
 59:   PetscFinalize();
 60:   return 0;
 61: }