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