1 2 static char help[] = "Tests MatTransposeMatMult() on MatLoad() matrix \n\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc, char **args) 7 { 8 Mat A, C, Bdense, Cdense; 9 PetscViewer fd; /* viewer */ 10 char file[PETSC_MAX_PATH_LEN]; /* input file name */ 11 PetscBool flg, viewmats = PETSC_FALSE; 12 PetscMPIInt rank, size; 13 PetscReal fill = 1.0; 14 PetscInt m, n, i, j, BN = 10, rstart, rend, *rows, *cols; 15 PetscScalar *Barray, *Carray, rval, *array; 16 Vec x, y; 17 PetscRandom rand; 18 19 PetscFunctionBeginUser; 20 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 21 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 22 23 /* Determine file from which we read the matrix A */ 24 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg)); 25 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option"); 26 27 /* Load matrix A */ 28 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd)); 29 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 30 PetscCall(MatLoad(A, fd)); 31 PetscCall(PetscViewerDestroy(&fd)); 32 33 /* Print (for testing only) */ 34 PetscCall(PetscOptionsHasName(NULL, NULL, "-view_mats", &viewmats)); 35 if (viewmats) { 36 if (rank == 0) printf("A_aij:\n"); 37 PetscCall(MatView(A, 0)); 38 } 39 40 /* Test MatTransposeMatMult_aij_aij() */ 41 PetscCall(MatTransposeMatMult(A, A, MAT_INITIAL_MATRIX, fill, &C)); 42 if (viewmats) { 43 if (rank == 0) printf("\nC = A_aij^T * A_aij:\n"); 44 PetscCall(MatView(C, 0)); 45 } 46 PetscCall(MatDestroy(&C)); 47 PetscCall(MatGetLocalSize(A, &m, &n)); 48 49 /* create a dense matrix Bdense */ 50 PetscCall(MatCreate(PETSC_COMM_WORLD, &Bdense)); 51 PetscCall(MatSetSizes(Bdense, m, PETSC_DECIDE, PETSC_DECIDE, BN)); 52 PetscCall(MatSetType(Bdense, MATDENSE)); 53 PetscCall(MatSetFromOptions(Bdense)); 54 PetscCall(MatSetUp(Bdense)); 55 PetscCall(MatGetOwnershipRange(Bdense, &rstart, &rend)); 56 57 PetscCall(PetscMalloc3(m, &rows, BN, &cols, m * BN, &array)); 58 for (i = 0; i < m; i++) rows[i] = rstart + i; 59 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand)); 60 PetscCall(PetscRandomSetFromOptions(rand)); 61 for (j = 0; j < BN; j++) { 62 cols[j] = j; 63 for (i = 0; i < m; i++) { 64 PetscCall(PetscRandomGetValue(rand, &rval)); 65 array[m * j + i] = rval; 66 } 67 } 68 PetscCall(MatSetValues(Bdense, m, rows, BN, cols, array, INSERT_VALUES)); 69 PetscCall(MatAssemblyBegin(Bdense, MAT_FINAL_ASSEMBLY)); 70 PetscCall(MatAssemblyEnd(Bdense, MAT_FINAL_ASSEMBLY)); 71 PetscCall(PetscRandomDestroy(&rand)); 72 PetscCall(PetscFree3(rows, cols, array)); 73 if (viewmats) { 74 if (rank == 0) printf("\nBdense:\n"); 75 PetscCall(MatView(Bdense, 0)); 76 } 77 78 /* Test MatTransposeMatMult_aij_dense() */ 79 PetscCall(MatTransposeMatMult(A, Bdense, MAT_INITIAL_MATRIX, fill, &C)); 80 PetscCall(MatTransposeMatMult(A, Bdense, MAT_REUSE_MATRIX, fill, &C)); 81 if (viewmats) { 82 if (rank == 0) printf("\nC=A^T*Bdense:\n"); 83 PetscCall(MatView(C, 0)); 84 } 85 86 /* Check accuracy */ 87 PetscCall(MatCreate(PETSC_COMM_WORLD, &Cdense)); 88 PetscCall(MatSetSizes(Cdense, n, PETSC_DECIDE, PETSC_DECIDE, BN)); 89 PetscCall(MatSetType(Cdense, MATDENSE)); 90 PetscCall(MatSetFromOptions(Cdense)); 91 PetscCall(MatSetUp(Cdense)); 92 PetscCall(MatAssemblyBegin(Cdense, MAT_FINAL_ASSEMBLY)); 93 PetscCall(MatAssemblyEnd(Cdense, MAT_FINAL_ASSEMBLY)); 94 95 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 96 if (size == 1) { 97 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, m, NULL, &x)); 98 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, n, NULL, &y)); 99 } else { 100 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, m, PETSC_DECIDE, NULL, &x)); 101 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, PETSC_DECIDE, NULL, &y)); 102 } 103 104 /* Cdense[:,j] = A^T * Bdense[:,j] */ 105 PetscCall(MatDenseGetArray(Bdense, &Barray)); 106 PetscCall(MatDenseGetArray(Cdense, &Carray)); 107 for (j = 0; j < BN; j++) { 108 PetscCall(VecPlaceArray(x, Barray)); 109 PetscCall(VecPlaceArray(y, Carray)); 110 111 PetscCall(MatMultTranspose(A, x, y)); 112 113 PetscCall(VecResetArray(x)); 114 PetscCall(VecResetArray(y)); 115 Barray += m; 116 Carray += n; 117 } 118 PetscCall(MatDenseRestoreArray(Bdense, &Barray)); 119 PetscCall(MatDenseRestoreArray(Cdense, &Carray)); 120 if (viewmats) { 121 if (rank == 0) printf("\nCdense:\n"); 122 PetscCall(MatView(Cdense, 0)); 123 } 124 125 PetscCall(MatEqual(C, Cdense, &flg)); 126 if (!flg) { 127 if (rank == 0) printf(" C != Cdense\n"); 128 } 129 130 /* Free data structures */ 131 PetscCall(MatDestroy(&A)); 132 PetscCall(MatDestroy(&C)); 133 PetscCall(MatDestroy(&Bdense)); 134 PetscCall(MatDestroy(&Cdense)); 135 PetscCall(VecDestroy(&x)); 136 PetscCall(VecDestroy(&y)); 137 PetscCall(PetscFinalize()); 138 return 0; 139 } 140 141 /*TEST 142 143 test: 144 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 145 args: -f ${DATAFILESPATH}/matrices/small 146 output_file: output/ex163.out 147 148 test: 149 suffix: 2 150 nsize: 3 151 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 152 args: -f ${DATAFILESPATH}/matrices/small 153 output_file: output/ex163.out 154 155 TEST*/ 156