1 static char help[] = "Creates a matrix, inserts some values, and tests MatCreateSubMatrices() and MatZeroEntries().\n\n";
2
3 #include <petscmat.h>
4
main(int argc,char ** argv)5 int main(int argc, char **argv)
6 {
7 Mat mat, submat, submat1, *submatrices;
8 PetscInt m = 10, n = 10, i = 4, tmp, rstart, rend;
9 IS irow, icol;
10 PetscScalar value = 1.0;
11 PetscViewer sviewer;
12 PetscBool allA = PETSC_FALSE;
13
14 PetscFunctionBeginUser;
15 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
16 PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_COMMON));
17 PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF, PETSC_VIEWER_ASCII_COMMON));
18
19 PetscCall(MatCreate(PETSC_COMM_WORLD, &mat));
20 PetscCall(MatSetSizes(mat, PETSC_DECIDE, PETSC_DECIDE, m, n));
21 PetscCall(MatSetFromOptions(mat));
22 PetscCall(MatSetUp(mat));
23 PetscCall(MatGetOwnershipRange(mat, &rstart, &rend));
24 for (i = rstart; i < rend; i++) {
25 value = (PetscReal)i + 1;
26 tmp = i % 5;
27 PetscCall(MatSetValues(mat, 1, &tmp, 1, &i, &value, INSERT_VALUES));
28 }
29 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
30 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
31 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "Original matrix\n"));
32 PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
33
34 /* Test MatCreateSubMatrix_XXX_All(), i.e., submatrix = A */
35 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_all", &allA, NULL));
36 if (allA) {
37 PetscCall(ISCreateStride(PETSC_COMM_SELF, m, 0, 1, &irow));
38 PetscCall(ISCreateStride(PETSC_COMM_SELF, n, 0, 1, &icol));
39 PetscCall(MatCreateSubMatrices(mat, 1, &irow, &icol, MAT_INITIAL_MATRIX, &submatrices));
40 PetscCall(MatCreateSubMatrices(mat, 1, &irow, &icol, MAT_REUSE_MATRIX, &submatrices));
41 submat = *submatrices;
42
43 /* sviewer will cause the submatrices (one per processor) to be printed in the correct order */
44 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "\nSubmatrices with all\n"));
45 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "--------------------\n"));
46 PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
47 PetscCall(MatView(submat, sviewer));
48 PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
49
50 PetscCall(ISDestroy(&irow));
51 PetscCall(ISDestroy(&icol));
52
53 /* test getting a reference on a submat */
54 PetscCall(PetscObjectReference((PetscObject)submat));
55 PetscCall(MatDestroySubMatrices(1, &submatrices));
56 PetscCall(MatDestroy(&submat));
57 }
58
59 /* Form submatrix with rows 2-4 and columns 4-8 */
60 PetscCall(ISCreateStride(PETSC_COMM_SELF, 3, 2, 1, &irow));
61 PetscCall(ISCreateStride(PETSC_COMM_SELF, 5, 4, 1, &icol));
62 PetscCall(MatCreateSubMatrices(mat, 1, &irow, &icol, MAT_INITIAL_MATRIX, &submatrices));
63 submat = *submatrices;
64
65 /* Test reuse submatrices */
66 PetscCall(MatCreateSubMatrices(mat, 1, &irow, &icol, MAT_REUSE_MATRIX, &submatrices));
67
68 /* sviewer will cause the submatrices (one per processor) to be printed in the correct order */
69 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "\nSubmatrices\n"));
70 PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
71 PetscCall(MatView(submat, sviewer));
72 PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
73 PetscCall(PetscObjectReference((PetscObject)submat));
74 PetscCall(MatDestroySubMatrices(1, &submatrices));
75 PetscCall(MatDestroy(&submat));
76
77 /* Form submatrix with rows 2-4 and all columns */
78 PetscCall(ISDestroy(&icol));
79 PetscCall(ISCreateStride(PETSC_COMM_SELF, 10, 0, 1, &icol));
80 PetscCall(MatCreateSubMatrices(mat, 1, &irow, &icol, MAT_INITIAL_MATRIX, &submatrices));
81 PetscCall(MatCreateSubMatrices(mat, 1, &irow, &icol, MAT_REUSE_MATRIX, &submatrices));
82 submat = *submatrices;
83
84 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "\nSubmatrices with allcolumns\n"));
85 PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
86 PetscCall(MatView(submat, sviewer));
87 PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
88
89 /* Test MatDuplicate */
90 PetscCall(MatDuplicate(submat, MAT_COPY_VALUES, &submat1));
91 PetscCall(MatDestroy(&submat1));
92
93 /* Zero the original matrix */
94 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "Original zeroed matrix\n"));
95 PetscCall(MatZeroEntries(mat));
96 PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));
97
98 PetscCall(ISDestroy(&irow));
99 PetscCall(ISDestroy(&icol));
100 PetscCall(PetscObjectReference((PetscObject)submat));
101 PetscCall(MatDestroySubMatrices(1, &submatrices));
102 PetscCall(MatDestroy(&submat));
103 PetscCall(MatDestroy(&mat));
104 PetscCall(PetscFinalize());
105 return 0;
106 }
107
108 /*TEST
109
110 test:
111 args: -mat_type aij
112
113 test:
114 suffix: 2
115 args: -mat_type dense
116
117 test:
118 suffix: 3
119 nsize: 3
120 args: -mat_type aij
121
122 test:
123 suffix: 4
124 nsize: 3
125 args: -mat_type dense
126
127 test:
128 suffix: 5
129 nsize: 3
130 args: -mat_type aij -test_all
131
132 TEST*/
133