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