Actual source code: ex209.c

  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];
 14:   PetscReal      fill = 4.0;
 15:   PetscMPIInt    rank,size;
 16:   PetscBool      equal;
 17:   PetscInt       i,m,n,rstart,rend;
 18:   PetscScalar    one=1.0;

 20:   PetscInitialize(&argc,&args,(char*)0,help);
 21:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 22:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 23:   PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),NULL);

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

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

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

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

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

 59:   /* Compare C and AtA */
 60:   MatMultEqual(C,AtA,20,&equal);
 62:   MatDestroy(&AtA);

 64:   MatDestroy(&C);
 65:   MatDestroy(&A);
 66:   PetscFinalize();
 67:   return 0;
 68: }

 70: /*TEST

 72:    test:
 73:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
 74:       args: -f ${DATAFILESPATH}/matrices/arco1

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

 82: TEST*/