Actual source code: ex4.c
1: static char help[] = "Creates a matrix, inserts some values, and tests MatCreateSubMatrices() and MatZeroEntries().\n\n";
3: #include <petscmat.h>
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;
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));
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));
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;
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));
50: PetscCall(ISDestroy(&irow));
51: PetscCall(ISDestroy(&icol));
53: /* test getting a reference on a submat */
54: PetscCall(PetscObjectReference((PetscObject)submat));
55: PetscCall(MatDestroySubMatrices(1, &submatrices));
56: PetscCall(MatDestroy(&submat));
57: }
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;
65: /* Test reuse submatrices */
66: PetscCall(MatCreateSubMatrices(mat, 1, &irow, &icol, MAT_REUSE_MATRIX, &submatrices));
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));
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;
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));
89: /* Test MatDuplicate */
90: PetscCall(MatDuplicate(submat, MAT_COPY_VALUES, &submat1));
91: PetscCall(MatDestroy(&submat1));
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));
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: }
108: /*TEST
110: test:
111: args: -mat_type aij
113: test:
114: suffix: 2
115: args: -mat_type dense
117: test:
118: suffix: 3
119: nsize: 3
120: args: -mat_type aij
122: test:
123: suffix: 4
124: nsize: 3
125: args: -mat_type dense
127: test:
128: suffix: 5
129: nsize: 3
130: args: -mat_type aij -test_all
132: TEST*/