Actual source code: ex31.c
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";
6: #include <petscksp.h>
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;
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));
26: PetscCall(PetscOptionsGetBool(NULL, NULL, "-partition", &partition, NULL));
27: PetscCall(PetscOptionsGetBool(NULL, NULL, "-displayIS", &displayIS, NULL));
28: PetscCall(PetscOptionsGetBool(NULL, NULL, "-displayMat", &displayMat, NULL));
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");
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);
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));
50: PetscCall(VecDuplicate(b, &x));
51: PetscCall(VecDuplicate(b, &u));
52: PetscCall(VecSet(x, 0.0));
54: /* - - - - - - - - - - - - - - - - - - - - - - - -
55: Test partition
56: - - - - - - - - - - - - - - - - - - - - - - - - - */
57: if (partition) {
58: MatPartitioning mpart;
59: IS mis, nis, is;
60: PetscInt *count;
61: Mat BB;
63: if (displayMat) {
64: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Before partitioning/reordering, A:\n"));
65: PetscCall(MatView(A, PETSC_VIEWER_DRAW_WORLD));
66: }
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: }
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: }
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: }
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: }
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: }
110: /* need to move the vector also */
111: PetscCall(ISDestroy(&is));
112: PetscCall(MatDestroy(&A));
113: A = BB;
114: }
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));
121: /* - - - - - - - - - - - - - - - - - - - - - - - -
122: Solve system
123: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
124: PetscCall(KSPSolve(ksp, b, x));
125: PetscCall(KSPGetIterationNumber(ksp, &its));
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: }
141: /* Free work space.*/
142: PetscCall(MatDestroy(&A));
143: PetscCall(VecDestroy(&b));
144: PetscCall(VecDestroy(&u));
145: PetscCall(VecDestroy(&x));
146: PetscCall(KSPDestroy(&ksp));
148: PetscCall(PetscFinalize());
149: return 0;
150: }
152: /*TEST
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
160: TEST*/