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