1: static char help[] = "Test MatTransposeMatMult() \n\n";

  3: /* Example:
  4:   mpiexec -n 8 ./ex209 -f <inputfile> (e.g., inputfile=ceres_solver_iteration_001_A.petsc)
  5: */

  7:  #include <petscmat.h>

  9: int main(int argc,char **args)
 10: {
 11:   Mat            A,C,AtA,B;
 12:   PetscViewer    fd;
 13:   char           file[PETSC_MAX_PATH_LEN];
 15:   PetscReal      fill = 4.0;
 16:   PetscMPIInt    rank,size;
 17:   PetscBool      equal;
 18:   PetscInt       i,m,n,rstart,rend;
 19:   PetscScalar    one=1.0;

 21:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 22:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 23:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 24:   PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,NULL);

 26:   /* Load matrix A */
 27:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
 28:   MatCreate(PETSC_COMM_WORLD,&A);
 29:   MatSetType(A,MATAIJ);
 30:   MatLoad(A,fd);
 31:   PetscViewerDestroy(&fd);

 33:   /* Create identity matrix B */
 34:   MatGetLocalSize(A,&m,&n);
 35:   MatCreate(PETSC_COMM_WORLD,&B);
 36:   MatSetSizes(B,m,m,PETSC_DECIDE,PETSC_DECIDE);
 37:   MatSetType(B,MATAIJ);
 38:   MatSetFromOptions(B);
 39:   MatSetUp(B);

 41:   MatGetOwnershipRange(B,&rstart,&rend);
 42:   for (i=rstart; i<rend; i++) {
 43:     MatSetValues(B,1,&i,1,&i,&one,INSERT_VALUES);
 44:   }
 45:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 46:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

 48:   /* Compute AtA = A^T*B*A, B = identity matrix */
 49:   MatPtAP(B,A,MAT_INITIAL_MATRIX,fill,&AtA);
 50:   MatPtAP(B,A,MAT_REUSE_MATRIX,fill,&AtA);
 51:   if (!rank) printf("C = A^T*B*A is done...\n");
 52:   MatDestroy(&B);

 54:   /* Compute C = A^T*A */
 55:   MatTransposeMatMult(A,A,MAT_INITIAL_MATRIX,fill,&C);
 56:   if (!rank) printf("C = A^T*A is done...\n");
 57:   MatTransposeMatMult(A,A,MAT_REUSE_MATRIX,fill,&C);
 58:   if (!rank) printf("REUSE C = A^T*A is done...\n");

 60:   /* Compare C and AtA */
 61:   MatMultEqual(C,AtA,20,&equal);
 62:   if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"A^T*A != At*A");
 63:   MatDestroy(&AtA);

 65:   MatDestroy(&C);
 66:   MatDestroy(&A);
 67:   PetscFinalize();
 68:   return ierr;
 69: }


 72: /*TEST

 74:    test:
 75:       requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
 76:       args: -f ${DATAFILESPATH}/matrices/arco1

 78:    test:
 79:       suffix: 2
 80:       nsize: 4
 81:       requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
 82:       args: -f ${DATAFILESPATH}/matrices/arco1 -matptap_via nonscalable

 84: TEST*/