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));
53: /* - - - - - - - - - - - - - - - - - - - - - - - -
54: Test partition
55: - - - - - - - - - - - - - - - - - - - - - - - - - */
56: if (partition) {
57: MatPartitioning mpart;
58: IS mis, nis, is;
59: PetscInt *count;
60: Mat BB;
62: if (displayMat) {
63: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Before partitioning/reordering, A:\n"));
64: PetscCall(MatView(A, PETSC_VIEWER_DRAW_WORLD));
65: }
67: PetscCall(PetscMalloc1(size, &count));
68: PetscCall(MatPartitioningCreate(PETSC_COMM_WORLD, &mpart));
69: PetscCall(MatPartitioningSetAdjacency(mpart, A));
70: /* PetscCall(MatPartitioningSetVertexWeights(mpart, weight)); */
71: PetscCall(MatPartitioningSetFromOptions(mpart));
72: PetscCall(MatPartitioningApply(mpart, &mis));
73: PetscCall(MatPartitioningDestroy(&mpart));
74: if (displayIS) {
75: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "mis, new processor assignment:\n"));
76: PetscCall(ISView(mis, PETSC_VIEWER_STDOUT_WORLD));
77: }
79: PetscCall(ISPartitioningToNumbering(mis, &nis));
80: if (displayIS) {
81: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "nis:\n"));
82: PetscCall(ISView(nis, PETSC_VIEWER_STDOUT_WORLD));
83: }
85: PetscCall(ISPartitioningCount(mis, size, count));
86: PetscCall(ISDestroy(&mis));
87: if (displayIS && rank == 0) {
88: PetscCall(PetscPrintf(PETSC_COMM_SELF, "[ %d ] count:\n", rank));
89: for (PetscInt i = 0; i < size; i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT, count[i]));
90: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n"));
91: }
93: PetscCall(ISInvertPermutation(nis, count[rank], &is));
94: PetscCall(PetscFree(count));
95: PetscCall(ISDestroy(&nis));
96: PetscCall(ISSort(is));
97: if (displayIS) {
98: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "inverse of nis - maps new local rows to old global rows:\n"));
99: PetscCall(ISView(is, PETSC_VIEWER_STDOUT_WORLD));
100: }
102: PetscCall(MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, &BB));
103: if (displayMat) {
104: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "After partitioning/reordering, A:\n"));
105: PetscCall(MatView(BB, PETSC_VIEWER_DRAW_WORLD));
106: }
108: /* need to move the vector also */
109: PetscCall(ISDestroy(&is));
110: PetscCall(MatDestroy(&A));
111: A = BB;
112: }
114: /* Create linear solver; set operators; set runtime options.*/
115: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
116: PetscCall(KSPSetOperators(ksp, A, A));
117: PetscCall(KSPSetFromOptions(ksp));
119: /* - - - - - - - - - - - - - - - - - - - - - - - -
120: Solve system
121: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
122: PetscCall(KSPSolve(ksp, b, x));
123: PetscCall(KSPGetIterationNumber(ksp, &its));
125: /* Check error */
126: PetscCall(MatMult(A, x, u));
127: PetscCall(VecAXPY(u, -1.0, b));
128: PetscCall(VecNorm(u, NORM_2, &norm));
129: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of iterations = %3" PetscInt_FMT "\n", its));
130: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Residual norm %g\n", (double)norm));
131: flg = PETSC_FALSE;
132: PetscCall(PetscOptionsGetBool(NULL, NULL, "-ksp_reason", &flg, NULL));
133: if (flg) {
134: KSPConvergedReason reason;
135: PetscCall(KSPGetConvergedReason(ksp, &reason));
136: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "KSPConvergedReason: %s\n", KSPConvergedReasons[reason]));
137: }
139: /* Free work space.*/
140: PetscCall(MatDestroy(&A));
141: PetscCall(VecDestroy(&b));
142: PetscCall(VecDestroy(&u));
143: PetscCall(VecDestroy(&x));
144: PetscCall(KSPDestroy(&ksp));
146: PetscCall(PetscFinalize());
147: return 0;
148: }
150: /*TEST
152: test:
153: args: -f ${DATAFILESPATH}/matrices/small -partition -mat_partitioning_type parmetis
154: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) parmetis
155: output_file: output/ex31.out
156: nsize: 3
158: TEST*/