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;

 21:   PetscInitialize(&argc, &args, (char *)0, help);
 22:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 23:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 24:   PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), 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++) MatSetValues(B, 1, &i, 1, &i, &one, INSERT_VALUES);
 43:   MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
 44:   MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);

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

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

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

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

 69: /*TEST

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

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

 81: TEST*/