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