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*/