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;
20: PetscInitialize(&argc, &args, (char *)0, help);
21: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
23: /* Determine file from which we read the matrix A */
24: PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg);
27: /* Load matrix A */
28: PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd);
29: MatCreate(PETSC_COMM_WORLD, &A);
30: MatLoad(A, fd);
31: PetscViewerDestroy(&fd);
33: /* Print (for testing only) */
34: PetscOptionsHasName(NULL, NULL, "-view_mats", &viewmats);
35: if (viewmats) {
36: if (rank == 0) printf("A_aij:\n");
37: MatView(A, 0);
38: }
40: /* Test MatTransposeMatMult_aij_aij() */
41: MatTransposeMatMult(A, A, MAT_INITIAL_MATRIX, fill, &C);
42: if (viewmats) {
43: if (rank == 0) printf("\nC = A_aij^T * A_aij:\n");
44: MatView(C, 0);
45: }
46: MatDestroy(&C);
47: MatGetLocalSize(A, &m, &n);
49: /* create a dense matrix Bdense */
50: MatCreate(PETSC_COMM_WORLD, &Bdense);
51: MatSetSizes(Bdense, m, PETSC_DECIDE, PETSC_DECIDE, BN);
52: MatSetType(Bdense, MATDENSE);
53: MatSetFromOptions(Bdense);
54: MatSetUp(Bdense);
55: MatGetOwnershipRange(Bdense, &rstart, &rend);
57: PetscMalloc3(m, &rows, BN, &cols, m * BN, &array);
58: for (i = 0; i < m; i++) rows[i] = rstart + i;
59: PetscRandomCreate(PETSC_COMM_WORLD, &rand);
60: PetscRandomSetFromOptions(rand);
61: for (j = 0; j < BN; j++) {
62: cols[j] = j;
63: for (i = 0; i < m; i++) {
64: PetscRandomGetValue(rand, &rval);
65: array[m * j + i] = rval;
66: }
67: }
68: MatSetValues(Bdense, m, rows, BN, cols, array, INSERT_VALUES);
69: MatAssemblyBegin(Bdense, MAT_FINAL_ASSEMBLY);
70: MatAssemblyEnd(Bdense, MAT_FINAL_ASSEMBLY);
71: PetscRandomDestroy(&rand);
72: PetscFree3(rows, cols, array);
73: if (viewmats) {
74: if (rank == 0) printf("\nBdense:\n");
75: MatView(Bdense, 0);
76: }
78: /* Test MatTransposeMatMult_aij_dense() */
79: MatTransposeMatMult(A, Bdense, MAT_INITIAL_MATRIX, fill, &C);
80: MatTransposeMatMult(A, Bdense, MAT_REUSE_MATRIX, fill, &C);
81: if (viewmats) {
82: if (rank == 0) printf("\nC=A^T*Bdense:\n");
83: MatView(C, 0);
84: }
86: /* Check accuracy */
87: MatCreate(PETSC_COMM_WORLD, &Cdense);
88: MatSetSizes(Cdense, n, PETSC_DECIDE, PETSC_DECIDE, BN);
89: MatSetType(Cdense, MATDENSE);
90: MatSetFromOptions(Cdense);
91: MatSetUp(Cdense);
92: MatAssemblyBegin(Cdense, MAT_FINAL_ASSEMBLY);
93: MatAssemblyEnd(Cdense, MAT_FINAL_ASSEMBLY);
95: MPI_Comm_size(PETSC_COMM_WORLD, &size);
96: if (size == 1) {
97: VecCreateSeqWithArray(PETSC_COMM_SELF, 1, m, NULL, &x);
98: VecCreateSeqWithArray(PETSC_COMM_SELF, 1, n, NULL, &y);
99: } else {
100: VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, m, PETSC_DECIDE, NULL, &x);
101: VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, PETSC_DECIDE, NULL, &y);
102: }
104: /* Cdense[:,j] = A^T * Bdense[:,j] */
105: MatDenseGetArray(Bdense, &Barray);
106: MatDenseGetArray(Cdense, &Carray);
107: for (j = 0; j < BN; j++) {
108: VecPlaceArray(x, Barray);
109: VecPlaceArray(y, Carray);
111: MatMultTranspose(A, x, y);
113: VecResetArray(x);
114: VecResetArray(y);
115: Barray += m;
116: Carray += n;
117: }
118: MatDenseRestoreArray(Bdense, &Barray);
119: MatDenseRestoreArray(Cdense, &Carray);
120: if (viewmats) {
121: if (rank == 0) printf("\nCdense:\n");
122: MatView(Cdense, 0);
123: }
125: MatEqual(C, Cdense, &flg);
126: if (!flg) {
127: if (rank == 0) printf(" C != Cdense\n");
128: }
130: /* Free data structures */
131: MatDestroy(&A);
132: MatDestroy(&C);
133: MatDestroy(&Bdense);
134: MatDestroy(&Cdense);
135: VecDestroy(&x);
136: VecDestroy(&y);
137: 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*/