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