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 Mat A, C, AtA, B; 11 PetscViewer fd; 12 char file[PETSC_MAX_PATH_LEN]; 13 PetscReal fill = 4.0; 14 PetscMPIInt rank, size; 15 PetscBool equal; 16 PetscInt i, m, n, rstart, rend; 17 PetscScalar one = 1.0; 18 19 PetscFunctionBeginUser; 20 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 21 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 22 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 23 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), NULL)); 24 25 /* Load matrix A */ 26 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd)); 27 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 28 PetscCall(MatSetType(A, MATAIJ)); 29 PetscCall(MatLoad(A, fd)); 30 PetscCall(PetscViewerDestroy(&fd)); 31 32 /* Create identity matrix B */ 33 PetscCall(MatGetLocalSize(A, &m, &n)); 34 PetscCall(MatCreate(PETSC_COMM_WORLD, &B)); 35 PetscCall(MatSetSizes(B, m, m, PETSC_DECIDE, PETSC_DECIDE)); 36 PetscCall(MatSetType(B, MATAIJ)); 37 PetscCall(MatSetFromOptions(B)); 38 PetscCall(MatSetUp(B)); 39 40 PetscCall(MatGetOwnershipRange(B, &rstart, &rend)); 41 for (i = rstart; i < rend; i++) PetscCall(MatSetValues(B, 1, &i, 1, &i, &one, INSERT_VALUES)); 42 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 43 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 44 45 /* Compute AtA = A^T*B*A, B = identity matrix */ 46 PetscCall(MatPtAP(B, A, MAT_INITIAL_MATRIX, fill, &AtA)); 47 PetscCall(MatPtAP(B, A, MAT_REUSE_MATRIX, fill, &AtA)); 48 if (rank == 0) printf("C = A^T*B*A is done...\n"); 49 PetscCall(MatDestroy(&B)); 50 51 /* Compute C = A^T*A */ 52 PetscCall(MatTransposeMatMult(A, A, MAT_INITIAL_MATRIX, fill, &C)); 53 if (rank == 0) printf("C = A^T*A is done...\n"); 54 PetscCall(MatTransposeMatMult(A, A, MAT_REUSE_MATRIX, fill, &C)); 55 if (rank == 0) printf("REUSE C = A^T*A is done...\n"); 56 57 /* Compare C and AtA */ 58 PetscCall(MatMultEqual(C, AtA, 20, &equal)); 59 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "A^T*A != At*A"); 60 PetscCall(MatDestroy(&AtA)); 61 62 PetscCall(MatDestroy(&C)); 63 PetscCall(MatDestroy(&A)); 64 PetscCall(PetscFinalize()); 65 return 0; 66 } 67 68 /*TEST 69 70 test: 71 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 72 args: -f ${DATAFILESPATH}/matrices/arco1 73 74 test: 75 suffix: 2 76 nsize: 4 77 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 78 args: -f ${DATAFILESPATH}/matrices/arco1 -matptap_via nonscalable 79 80 TEST*/ 81