1 static char help[] = "Test partition. Reads a PETSc matrix and vector from a file and solves a linear system.\n\
2 This Input parameters include\n\
3 -f <input_file> : file to load \n\
4 -partition -mat_partitioning_view \n\\n";
5
6 #include <petscksp.h>
7
main(int argc,char ** args)8 int main(int argc, char **args)
9 {
10 KSP ksp; /* linear solver context */
11 Mat A; /* matrix */
12 Vec x, b, u; /* approx solution, RHS, exact solution */
13 PetscViewer fd; /* viewer */
14 char file[PETSC_MAX_PATH_LEN]; /* input file name */
15 PetscBool flg, partition = PETSC_FALSE, displayIS = PETSC_FALSE, displayMat = PETSC_FALSE;
16 PetscInt its, m, n;
17 PetscReal norm;
18 PetscMPIInt size, rank;
19 PetscScalar one = 1.0;
20
21 PetscFunctionBeginUser;
22 PetscCall(PetscInitialize(&argc, &args, NULL, help));
23 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
24 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
25
26 PetscCall(PetscOptionsGetBool(NULL, NULL, "-partition", &partition, NULL));
27 PetscCall(PetscOptionsGetBool(NULL, NULL, "-displayIS", &displayIS, NULL));
28 PetscCall(PetscOptionsGetBool(NULL, NULL, "-displayMat", &displayMat, NULL));
29
30 /* Determine file from which we read the matrix.*/
31 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg));
32 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "Must indicate binary file with the -f option");
33
34 /* - - - - - - - - - - - - - - - - - - - - - - - -
35 Load system
36 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
37 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
38 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
39 PetscCall(MatLoad(A, fd));
40 PetscCall(PetscViewerDestroy(&fd));
41 PetscCall(MatGetLocalSize(A, &m, &n));
42 PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")", m, n);
43
44 /* Create rhs vector of all ones */
45 PetscCall(VecCreate(PETSC_COMM_WORLD, &b));
46 PetscCall(VecSetSizes(b, m, PETSC_DECIDE));
47 PetscCall(VecSetFromOptions(b));
48 PetscCall(VecSet(b, one));
49
50 PetscCall(VecDuplicate(b, &x));
51 PetscCall(VecDuplicate(b, &u));
52 PetscCall(VecSet(x, 0.0));
53
54 /* - - - - - - - - - - - - - - - - - - - - - - - -
55 Test partition
56 - - - - - - - - - - - - - - - - - - - - - - - - - */
57 if (partition) {
58 MatPartitioning mpart;
59 IS mis, nis, is;
60 PetscInt *count;
61 Mat BB;
62
63 if (displayMat) {
64 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Before partitioning/reordering, A:\n"));
65 PetscCall(MatView(A, PETSC_VIEWER_DRAW_WORLD));
66 }
67
68 PetscCall(PetscMalloc1(size, &count));
69 PetscCall(MatPartitioningCreate(PETSC_COMM_WORLD, &mpart));
70 PetscCall(MatPartitioningSetAdjacency(mpart, A));
71 /* PetscCall(MatPartitioningSetVertexWeights(mpart, weight)); */
72 PetscCall(MatPartitioningSetFromOptions(mpart));
73 PetscCall(MatPartitioningApply(mpart, &mis));
74 PetscCall(MatPartitioningDestroy(&mpart));
75 if (displayIS) {
76 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "mis, new processor assignment:\n"));
77 PetscCall(ISView(mis, PETSC_VIEWER_STDOUT_WORLD));
78 }
79
80 PetscCall(ISPartitioningToNumbering(mis, &nis));
81 if (displayIS) {
82 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "nis:\n"));
83 PetscCall(ISView(nis, PETSC_VIEWER_STDOUT_WORLD));
84 }
85
86 PetscCall(ISPartitioningCount(mis, size, count));
87 PetscCall(ISDestroy(&mis));
88 if (displayIS && rank == 0) {
89 PetscInt i;
90 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[ %d ] count:\n", rank));
91 for (i = 0; i < size; i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT, count[i]));
92 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n"));
93 }
94
95 PetscCall(ISInvertPermutation(nis, count[rank], &is));
96 PetscCall(PetscFree(count));
97 PetscCall(ISDestroy(&nis));
98 PetscCall(ISSort(is));
99 if (displayIS) {
100 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "inverse of nis - maps new local rows to old global rows:\n"));
101 PetscCall(ISView(is, PETSC_VIEWER_STDOUT_WORLD));
102 }
103
104 PetscCall(MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, &BB));
105 if (displayMat) {
106 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "After partitioning/reordering, A:\n"));
107 PetscCall(MatView(BB, PETSC_VIEWER_DRAW_WORLD));
108 }
109
110 /* need to move the vector also */
111 PetscCall(ISDestroy(&is));
112 PetscCall(MatDestroy(&A));
113 A = BB;
114 }
115
116 /* Create linear solver; set operators; set runtime options.*/
117 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
118 PetscCall(KSPSetOperators(ksp, A, A));
119 PetscCall(KSPSetFromOptions(ksp));
120
121 /* - - - - - - - - - - - - - - - - - - - - - - - -
122 Solve system
123 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
124 PetscCall(KSPSolve(ksp, b, x));
125 PetscCall(KSPGetIterationNumber(ksp, &its));
126
127 /* Check error */
128 PetscCall(MatMult(A, x, u));
129 PetscCall(VecAXPY(u, -1.0, b));
130 PetscCall(VecNorm(u, NORM_2, &norm));
131 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of iterations = %3" PetscInt_FMT "\n", its));
132 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Residual norm %g\n", (double)norm));
133 flg = PETSC_FALSE;
134 PetscCall(PetscOptionsGetBool(NULL, NULL, "-ksp_reason", &flg, NULL));
135 if (flg) {
136 KSPConvergedReason reason;
137 PetscCall(KSPGetConvergedReason(ksp, &reason));
138 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "KSPConvergedReason: %s\n", KSPConvergedReasons[reason]));
139 }
140
141 /* Free work space.*/
142 PetscCall(MatDestroy(&A));
143 PetscCall(VecDestroy(&b));
144 PetscCall(VecDestroy(&u));
145 PetscCall(VecDestroy(&x));
146 PetscCall(KSPDestroy(&ksp));
147
148 PetscCall(PetscFinalize());
149 return 0;
150 }
151
152 /*TEST
153
154 test:
155 args: -f ${DATAFILESPATH}/matrices/small -partition -mat_partitioning_type parmetis
156 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) parmetis
157 output_file: output/ex31.out
158 nsize: 3
159
160 TEST*/
161