1 static char help[] = "Test MatTransposeMatMult() \n\n"; 2 3 /* Example: 4 mpiexec -n 8 ./ex209 -f <inputfile> (e.g., inputfile=ceres_solver_iteration_001_A.petsc) 5 */ 6 7 #include <petscmat.h> 8 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; 19 20 PetscFunctionBeginUser; 21 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 22 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 23 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 24 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), NULL)); 25 26 /* Load matrix A */ 27 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd)); 28 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 29 PetscCall(MatSetType(A, MATAIJ)); 30 PetscCall(MatLoad(A, fd)); 31 PetscCall(PetscViewerDestroy(&fd)); 32 33 /* Create identity matrix B */ 34 PetscCall(MatGetLocalSize(A, &m, &n)); 35 PetscCall(MatCreate(PETSC_COMM_WORLD, &B)); 36 PetscCall(MatSetSizes(B, m, m, PETSC_DECIDE, PETSC_DECIDE)); 37 PetscCall(MatSetType(B, MATAIJ)); 38 PetscCall(MatSetFromOptions(B)); 39 PetscCall(MatSetUp(B)); 40 41 PetscCall(MatGetOwnershipRange(B, &rstart, &rend)); 42 for (i = rstart; i < rend; i++) PetscCall(MatSetValues(B, 1, &i, 1, &i, &one, INSERT_VALUES)); 43 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 44 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 45 46 /* Compute AtA = A^T*B*A, B = identity matrix */ 47 PetscCall(MatPtAP(B, A, MAT_INITIAL_MATRIX, fill, &AtA)); 48 PetscCall(MatPtAP(B, A, MAT_REUSE_MATRIX, fill, &AtA)); 49 if (rank == 0) printf("C = A^T*B*A is done...\n"); 50 PetscCall(MatDestroy(&B)); 51 52 /* Compute C = A^T*A */ 53 PetscCall(MatTransposeMatMult(A, A, MAT_INITIAL_MATRIX, fill, &C)); 54 if (rank == 0) printf("C = A^T*A is done...\n"); 55 PetscCall(MatTransposeMatMult(A, A, MAT_REUSE_MATRIX, fill, &C)); 56 if (rank == 0) printf("REUSE C = A^T*A is done...\n"); 57 58 /* Compare C and AtA */ 59 PetscCall(MatMultEqual(C, AtA, 20, &equal)); 60 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "A^T*A != At*A"); 61 PetscCall(MatDestroy(&AtA)); 62 63 PetscCall(MatDestroy(&C)); 64 PetscCall(MatDestroy(&A)); 65 PetscCall(PetscFinalize()); 66 return 0; 67 } 68 69 /*TEST 70 71 test: 72 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 73 args: -f ${DATAFILESPATH}/matrices/arco1 74 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 80 81 TEST*/ 82