1c4762a1bSJed Brown static char help[] = "Test MatTransposeMatMult() \n\n";
2c4762a1bSJed Brown
3c4762a1bSJed Brown /* Example:
4c4762a1bSJed Brown mpiexec -n 8 ./ex209 -f <inputfile> (e.g., inputfile=ceres_solver_iteration_001_A.petsc)
5c4762a1bSJed Brown */
6c4762a1bSJed Brown
7c4762a1bSJed Brown #include <petscmat.h>
8c4762a1bSJed Brown
main(int argc,char ** args)9d71ae5a4SJacob Faibussowitsch int main(int argc, char **args)
10d71ae5a4SJacob Faibussowitsch {
11c4762a1bSJed Brown Mat A, C, AtA, B;
12c4762a1bSJed Brown PetscViewer fd;
13c4762a1bSJed Brown char file[PETSC_MAX_PATH_LEN];
14c4762a1bSJed Brown PetscReal fill = 4.0;
15c4762a1bSJed Brown PetscMPIInt rank, size;
16c4762a1bSJed Brown PetscBool equal;
17c4762a1bSJed Brown PetscInt i, m, n, rstart, rend;
18c4762a1bSJed Brown PetscScalar one = 1.0;
19c4762a1bSJed Brown
20327415f7SBarry Smith PetscFunctionBeginUser;
21*c8025a54SPierre Jolivet PetscCall(PetscInitialize(&argc, &args, NULL, help));
229566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
239566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
249566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), NULL));
25c4762a1bSJed Brown
26c4762a1bSJed Brown /* Load matrix A */
279566063dSJacob Faibussowitsch PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
289566063dSJacob Faibussowitsch PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
299566063dSJacob Faibussowitsch PetscCall(MatSetType(A, MATAIJ));
309566063dSJacob Faibussowitsch PetscCall(MatLoad(A, fd));
319566063dSJacob Faibussowitsch PetscCall(PetscViewerDestroy(&fd));
32c4762a1bSJed Brown
33c4762a1bSJed Brown /* Create identity matrix B */
349566063dSJacob Faibussowitsch PetscCall(MatGetLocalSize(A, &m, &n));
359566063dSJacob Faibussowitsch PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
369566063dSJacob Faibussowitsch PetscCall(MatSetSizes(B, m, m, PETSC_DECIDE, PETSC_DECIDE));
379566063dSJacob Faibussowitsch PetscCall(MatSetType(B, MATAIJ));
389566063dSJacob Faibussowitsch PetscCall(MatSetFromOptions(B));
399566063dSJacob Faibussowitsch PetscCall(MatSetUp(B));
40c4762a1bSJed Brown
419566063dSJacob Faibussowitsch PetscCall(MatGetOwnershipRange(B, &rstart, &rend));
4248a46eb9SPierre Jolivet for (i = rstart; i < rend; i++) PetscCall(MatSetValues(B, 1, &i, 1, &i, &one, INSERT_VALUES));
439566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
449566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
45c4762a1bSJed Brown
46c4762a1bSJed Brown /* Compute AtA = A^T*B*A, B = identity matrix */
479566063dSJacob Faibussowitsch PetscCall(MatPtAP(B, A, MAT_INITIAL_MATRIX, fill, &AtA));
489566063dSJacob Faibussowitsch PetscCall(MatPtAP(B, A, MAT_REUSE_MATRIX, fill, &AtA));
49dd400576SPatrick Sanan if (rank == 0) printf("C = A^T*B*A is done...\n");
509566063dSJacob Faibussowitsch PetscCall(MatDestroy(&B));
51c4762a1bSJed Brown
52c4762a1bSJed Brown /* Compute C = A^T*A */
539566063dSJacob Faibussowitsch PetscCall(MatTransposeMatMult(A, A, MAT_INITIAL_MATRIX, fill, &C));
54dd400576SPatrick Sanan if (rank == 0) printf("C = A^T*A is done...\n");
559566063dSJacob Faibussowitsch PetscCall(MatTransposeMatMult(A, A, MAT_REUSE_MATRIX, fill, &C));
56dd400576SPatrick Sanan if (rank == 0) printf("REUSE C = A^T*A is done...\n");
57c4762a1bSJed Brown
58c4762a1bSJed Brown /* Compare C and AtA */
599566063dSJacob Faibussowitsch PetscCall(MatMultEqual(C, AtA, 20, &equal));
6028b400f6SJacob Faibussowitsch PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "A^T*A != At*A");
619566063dSJacob Faibussowitsch PetscCall(MatDestroy(&AtA));
62c4762a1bSJed Brown
639566063dSJacob Faibussowitsch PetscCall(MatDestroy(&C));
649566063dSJacob Faibussowitsch PetscCall(MatDestroy(&A));
659566063dSJacob Faibussowitsch PetscCall(PetscFinalize());
66b122ec5aSJacob Faibussowitsch return 0;
67c4762a1bSJed Brown }
68c4762a1bSJed Brown
69c4762a1bSJed Brown /*TEST
70c4762a1bSJed Brown
71c4762a1bSJed Brown test:
72dfd57a17SPierre Jolivet requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
73c4762a1bSJed Brown args: -f ${DATAFILESPATH}/matrices/arco1
74c4762a1bSJed Brown
75c4762a1bSJed Brown test:
76c4762a1bSJed Brown suffix: 2
77c4762a1bSJed Brown nsize: 4
78dfd57a17SPierre Jolivet requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
79c4762a1bSJed Brown args: -f ${DATAFILESPATH}/matrices/arco1 -matptap_via nonscalable
80c4762a1bSJed Brown
81c4762a1bSJed Brown TEST*/
82