Actual source code: ex33.c
1: static char help[] = "Test memory scalability of MatMatMult() for AIJ and DENSE matrices. \n\
2: Modified from the code contributed by Ian Lin <iancclin@umich.edu> \n\n";
4: /*
5: Example:
6: mpiexec -n <np> ./ex33 -mem_view -matproduct_batch_size <Bbn>
7: */
9: #include <petsc.h>
11: PetscErrorCode Print_memory(PetscLogDouble mem)
12: {
13: double max_mem, min_mem;
15: PetscFunctionBeginUser;
16: PetscCallMPI(MPI_Reduce(&mem, &max_mem, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD));
17: PetscCallMPI(MPI_Reduce(&mem, &min_mem, 1, MPI_DOUBLE, MPI_MIN, 0, MPI_COMM_WORLD));
18: max_mem = max_mem / 1024.0 / 1024.0;
19: min_mem = min_mem / 1024.0 / 1024.0;
20: PetscCall(PetscPrintf(MPI_COMM_WORLD, " max and min memory across all processors %.4f Mb, %.4f Mb.\n", (double)max_mem, (double)min_mem));
21: PetscFunctionReturn(PETSC_SUCCESS);
22: }
24: /*
25: Illustrate how to use MPI derived data types.
26: It would save memory significantly. See MatMPIDenseScatter()
27: */
28: PetscErrorCode TestMPIDerivedDataType(void)
29: {
30: MPI_Datatype type1, type2, rtype1, rtype2;
31: PetscScalar buffer[24]; /* An array of 4 rows, 6 cols */
32: MPI_Status status;
33: PetscMPIInt rank, size, disp[2];
35: PetscFunctionBeginUser;
36: PetscCallMPI(MPI_Comm_size(MPI_COMM_WORLD, &size));
37: PetscCheck(size >= 2, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, "Must use at least 2 processors");
38: PetscCallMPI(MPI_Comm_rank(MPI_COMM_WORLD, &rank));
40: if (rank == 0) {
41: /* proc[0] sends 2 rows to proc[1] */
42: for (PetscInt i = 0; i < 24; i++) buffer[i] = (PetscScalar)i;
44: disp[0] = 0;
45: disp[1] = 2;
46: PetscCallMPI(MPI_Type_create_indexed_block(2, 1, (const PetscMPIInt *)disp, MPIU_SCALAR, &type1));
47: /* one column has 4 entries */
48: PetscCallMPI(MPI_Type_create_resized(type1, 0, 4 * sizeof(PetscScalar), &type2));
49: PetscCallMPI(MPI_Type_commit(&type2));
50: PetscCallMPI(MPI_Send(buffer, 6, type2, 1, 123, MPI_COMM_WORLD));
52: } else if (rank == 1) {
53: /* proc[1] receives 2 rows from proc[0], and put them into contiguous rows, starting at the row 1 (disp[0]) */
54: for (PetscInt i = 0; i < 24; i++) buffer[i] = 0.0;
56: disp[0] = 1;
57: PetscCallMPI(MPI_Type_create_indexed_block(1, 2, (const PetscMPIInt *)disp, MPIU_SCALAR, &rtype1));
58: PetscCallMPI(MPI_Type_create_resized(rtype1, 0, 4 * sizeof(PetscScalar), &rtype2));
60: PetscCallMPI(MPI_Type_commit(&rtype2));
61: PetscCallMPI(MPI_Recv(buffer, 6, rtype2, 0, 123, MPI_COMM_WORLD, &status));
62: for (PetscInt i = 0; i < 4; i++) {
63: for (PetscInt j = 0; j < 6; j++) PetscCall(PetscPrintf(MPI_COMM_SELF, " %g", (double)PetscRealPart(buffer[i + j * 4])));
64: PetscCall(PetscPrintf(MPI_COMM_SELF, "\n"));
65: }
66: }
68: if (rank == 0) {
69: PetscCallMPI(MPI_Type_free(&type1));
70: PetscCallMPI(MPI_Type_free(&type2));
71: } else if (rank == 1) {
72: PetscCallMPI(MPI_Type_free(&rtype1));
73: PetscCallMPI(MPI_Type_free(&rtype2));
74: }
75: PetscCallMPI(MPI_Barrier(MPI_COMM_WORLD));
76: PetscFunctionReturn(PETSC_SUCCESS);
77: }
79: int main(int argc, char **args)
80: {
81: PetscInt mA = 2700, nX = 80, nz = 40;
82: /* PetscInt mA=6,nX=5,nz=2; //small test */
83: PetscLogDouble mem;
84: Mat A, X, Y;
85: PetscBool flg = PETSC_FALSE;
87: PetscFunctionBeginUser;
88: PetscCall(PetscInitialize(&argc, &args, NULL, help));
89: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_mpiderivedtype", &flg, NULL));
90: if (flg) {
91: PetscCall(TestMPIDerivedDataType());
92: PetscCall(PetscFinalize());
93: return 0;
94: }
96: PetscCall(PetscOptionsGetBool(NULL, NULL, "-mem_view", &flg, NULL));
97: PetscCall(PetscMemoryGetCurrentUsage(&mem));
98: if (flg) {
99: PetscCall(PetscPrintf(MPI_COMM_WORLD, "Before start,"));
100: PetscCall(Print_memory(mem));
101: }
103: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, mA, mA, nz, NULL, nz, NULL, &A));
104: PetscCall(MatSetRandom(A, NULL));
105: PetscCall(PetscMemoryGetCurrentUsage(&mem));
106: if (flg) {
107: PetscCall(PetscPrintf(MPI_COMM_WORLD, "After creating A,"));
108: PetscCall(Print_memory(mem));
109: }
111: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, mA, nX, NULL, &X));
112: PetscCall(MatSetRandom(X, NULL));
113: PetscCall(PetscMemoryGetCurrentUsage(&mem));
114: if (flg) {
115: PetscCall(PetscPrintf(MPI_COMM_WORLD, "After creating X,"));
116: PetscCall(Print_memory(mem));
117: }
119: PetscCall(MatMatMult(A, X, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Y));
120: PetscCall(PetscMemoryGetCurrentUsage(&mem));
121: if (flg) {
122: PetscCall(PetscPrintf(MPI_COMM_WORLD, "After MatMatMult,"));
123: PetscCall(Print_memory(mem));
124: }
126: /* Test reuse */
127: PetscCall(MatMatMult(A, X, MAT_REUSE_MATRIX, PETSC_DETERMINE, &Y));
128: PetscCall(PetscMemoryGetCurrentUsage(&mem));
129: if (flg) {
130: PetscCall(PetscPrintf(MPI_COMM_WORLD, "After reuse MatMatMult,"));
131: PetscCall(Print_memory(mem));
132: }
134: /* Check accuracy */
135: PetscCall(MatMatMultEqual(A, X, Y, 10, &flg));
136: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "Error in MatMatMult()");
138: PetscCall(MatDestroy(&A));
139: PetscCall(MatDestroy(&X));
140: PetscCall(MatDestroy(&Y));
142: PetscCall(PetscFinalize());
143: return 0;
144: }
146: /*TEST
148: test:
149: suffix: 1
150: nsize: 4
151: output_file: output/empty.out
153: test:
154: suffix: 2
155: nsize: 8
156: output_file: output/empty.out
158: test:
159: suffix: 3
160: nsize: 2
161: args: -test_mpiderivedtype
162: output_file: output/ex33_3.out
164: TEST*/