1 static char help[] = "Tests MatCreateSubmatrix() in parallel.";
2
3 #include <petscmat.h>
4
main(int argc,char ** args)5 int main(int argc, char **args)
6 {
7 Mat C, A;
8 PetscInt i, j, m = 3, n = 2, rstart, rend;
9 PetscMPIInt size, rank;
10 PetscScalar v;
11 IS isrow, iscol;
12 PetscBool test_matmatmult = PETSC_FALSE;
13
14 PetscFunctionBeginUser;
15 PetscCall(PetscInitialize(&argc, &args, NULL, help));
16 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matmatmult", &test_matmatmult, NULL));
17
18 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
19 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
20 n = 2 * size;
21
22 PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
23 PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n));
24 PetscCall(MatSetFromOptions(C));
25 PetscCall(MatSetUp(C));
26
27 /*
28 This is JUST to generate a nice test matrix, all processors fill up
29 the entire matrix. This is not something one would ever do in practice.
30 */
31 PetscCall(MatGetOwnershipRange(C, &rstart, &rend));
32 for (i = rstart; i < rend; i++) {
33 for (j = 0; j < m * n; j++) {
34 v = i + j + 1;
35 PetscCall(MatSetValues(C, 1, &i, 1, &j, &v, INSERT_VALUES));
36 }
37 }
38
39 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
40 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
41 PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_COMMON));
42 PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
43
44 /*
45 Generate a new matrix consisting of every second row and column of
46 the original matrix
47 */
48 PetscCall(MatGetOwnershipRange(C, &rstart, &rend));
49 /* Create parallel IS with the rows we want on THIS processor */
50 PetscCall(ISCreateStride(PETSC_COMM_WORLD, (rend - rstart) / 2, rstart, 2, &isrow));
51 /* Create parallel IS with the rows we want on THIS processor (same as rows for now) */
52 PetscCall(ISCreateStride(PETSC_COMM_WORLD, (rend - rstart) / 2, rstart, 2, &iscol));
53
54 PetscCall(MatCreateSubMatrix(C, isrow, iscol, MAT_INITIAL_MATRIX, &A));
55 PetscCall(MatCreateSubMatrix(C, isrow, iscol, MAT_REUSE_MATRIX, &A));
56 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
57
58 if (test_matmatmult) {
59 PetscCall(MatDestroy(&C));
60 PetscCall(MatMatMult(A, A, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C));
61 PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
62 }
63
64 PetscCall(ISDestroy(&isrow));
65 PetscCall(ISDestroy(&iscol));
66 PetscCall(MatDestroy(&A));
67 PetscCall(MatDestroy(&C));
68 PetscCall(PetscFinalize());
69 return 0;
70 }
71
72 /*TEST
73
74 test:
75
76 test:
77 suffix: 2
78 nsize: 3
79
80 test:
81 suffix: 2_baij
82 nsize: 3
83 args: -mat_type baij
84
85 test:
86 suffix: 2_sbaij
87 nsize: 3
88 args: -mat_type sbaij
89
90 test:
91 suffix: baij
92 args: -mat_type baij
93 output_file: output/ex59_1_baij.out
94
95 test:
96 suffix: sbaij
97 args: -mat_type sbaij
98 output_file: output/ex59_1_sbaij.out
99
100 test:
101 suffix: kok
102 nsize: 3
103 requires: kokkos_kernels
104 args: -mat_type aijkokkos -test_matmatmult
105 filter: grep -v -i type
106 output_file: output/ex59_kok.out
107
108 TEST*/
109