1: static char help[] = "Test MatPtAP,  MatMatMatMult\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
  9: int main(int argc,char **args)
 10: {
 11:   Mat            A,P,C,R,RAP;
 12:   PetscViewer    fd;
 13:   char           file[2][PETSC_MAX_PATH_LEN];
 14:   PetscBool      flg,testPtAP=PETSC_TRUE,testRARt=PETSC_TRUE;
 16:   PetscReal      fill=2.0,norm;

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

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

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

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

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

 54:   PetscOptionsGetBool(NULL,NULL,"-testPtAP",&testPtAP,NULL);
 55:   PetscOptionsGetBool(NULL,NULL,"-testRARt",&testRARt,NULL);

 57:   MatTranspose(P,MAT_INITIAL_MATRIX,&R);

 59:   if (testPtAP) {
 60:     MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);
 61:     MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);

 63:     /* Check PtAP = RAP */
 64:     MatMatMatMult(R,A,P,MAT_INITIAL_MATRIX,2.0,&RAP);
 65:     MatAXPY(C,-1.0,RAP,DIFFERENT_NONZERO_PATTERN);
 66:     MatNorm(C,NORM_FROBENIUS,&norm);
 67:     if (norm > 1.e-14) {PetscPrintf(PETSC_COMM_SELF,"norm(PtAP - RAP)= %g\n",norm);}
 68:     MatDestroy(&C);
 69:     MatDestroy(&RAP);
 70:   }

 72:   if (testRARt) {
 73:     MatRARt(A,R,MAT_INITIAL_MATRIX,fill,&C);
 74:     MatRARt(A,R,MAT_REUSE_MATRIX,fill,&C);

 76:     /* Check RARt = RAP */
 77:     MatMatMatMult(R,A,P,MAT_INITIAL_MATRIX,2.0,&RAP);
 78:     MatAXPY(C,-1.0,RAP,DIFFERENT_NONZERO_PATTERN);
 79:     MatNorm(C,NORM_FROBENIUS,&norm);
 80:     if (norm > 1.e-14) {PetscPrintf(PETSC_COMM_SELF,"norm(RARt - RAP)= %g\n",norm);}
 81:     MatDestroy(&C);
 82:     MatDestroy(&RAP);
 83:   }

 85:   MatDestroy(&R);
 86:   MatDestroy(&P);
 87:   MatDestroy(&A);
 88:   PetscFinalize();
 89:   return ierr;
 90: }