Actual source code: ex163.c
2: static char help[] = "Tests MatTransposeMatMult() on MatLoad() matrix \n\n";
4: #include <petscmat.h>
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;
19: PetscFunctionBeginUser;
20: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
21: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
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");
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));
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: }
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));
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));
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: }
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: }
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));
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: }
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));
111: PetscCall(MatMultTranspose(A, x, y));
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: }
125: PetscCall(MatEqual(C, Cdense, &flg));
126: if (!flg) {
127: if (rank == 0) printf(" C != Cdense\n");
128: }
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: }
141: /*TEST
143: test:
144: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
145: args: -f ${DATAFILESPATH}/matrices/small
146: output_file: output/ex163.out
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
155: TEST*/