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;
23: PetscInitialize(&argc, &args, (char *)0, help);
24: MPI_Comm_size(PETSC_COMM_WORLD, &size);
25: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
27: PetscOptionsGetBool(NULL, NULL, "-partition", &partition, NULL);
28: PetscOptionsGetBool(NULL, NULL, "-displayIS", &displayIS, NULL);
29: PetscOptionsGetBool(NULL, NULL, "-displayMat", &displayMat, NULL);
31: /* Determine file from which we read the matrix.*/
32: PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg);
35: /* - - - - - - - - - - - - - - - - - - - - - - - -
36: Load system
37: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
38: PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd);
39: MatCreate(PETSC_COMM_WORLD, &A);
40: MatLoad(A, fd);
41: PetscViewerDestroy(&fd);
42: MatGetLocalSize(A, &m, &n);
45: /* Create rhs vector of all ones */
46: VecCreate(PETSC_COMM_WORLD, &b);
47: VecSetSizes(b, m, PETSC_DECIDE);
48: VecSetFromOptions(b);
49: VecSet(b, one);
51: VecDuplicate(b, &x);
52: VecDuplicate(b, &u);
53: 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: PetscPrintf(PETSC_COMM_WORLD, "Before partitioning/reordering, A:\n");
66: MatView(A, PETSC_VIEWER_DRAW_WORLD);
67: }
69: PetscMalloc1(size, &count);
70: MatPartitioningCreate(PETSC_COMM_WORLD, &mpart);
71: MatPartitioningSetAdjacency(mpart, A);
72: /* MatPartitioningSetVertexWeights(mpart, weight); */
73: MatPartitioningSetFromOptions(mpart);
74: MatPartitioningApply(mpart, &mis);
75: MatPartitioningDestroy(&mpart);
76: if (displayIS) {
77: PetscPrintf(PETSC_COMM_WORLD, "mis, new processor assignment:\n");
78: ISView(mis, PETSC_VIEWER_STDOUT_WORLD);
79: }
81: ISPartitioningToNumbering(mis, &nis);
82: if (displayIS) {
83: PetscPrintf(PETSC_COMM_WORLD, "nis:\n");
84: ISView(nis, PETSC_VIEWER_STDOUT_WORLD);
85: }
87: ISPartitioningCount(mis, size, count);
88: ISDestroy(&mis);
89: if (displayIS && rank == 0) {
90: PetscInt i;
91: PetscPrintf(PETSC_COMM_SELF, "[ %d ] count:\n", rank);
92: for (i = 0; i < size; i++) PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT, count[i]);
93: PetscPrintf(PETSC_COMM_WORLD, "\n");
94: }
96: ISInvertPermutation(nis, count[rank], &is);
97: PetscFree(count);
98: ISDestroy(&nis);
99: ISSort(is);
100: if (displayIS) {
101: PetscPrintf(PETSC_COMM_WORLD, "inverse of nis - maps new local rows to old global rows:\n");
102: ISView(is, PETSC_VIEWER_STDOUT_WORLD);
103: }
105: MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, &BB);
106: if (displayMat) {
107: PetscPrintf(PETSC_COMM_WORLD, "After partitioning/reordering, A:\n");
108: MatView(BB, PETSC_VIEWER_DRAW_WORLD);
109: }
111: /* need to move the vector also */
112: ISDestroy(&is);
113: MatDestroy(&A);
114: A = BB;
115: }
117: /* Create linear solver; set operators; set runtime options.*/
118: KSPCreate(PETSC_COMM_WORLD, &ksp);
119: KSPSetOperators(ksp, A, A);
120: KSPSetFromOptions(ksp);
122: /* - - - - - - - - - - - - - - - - - - - - - - - -
123: Solve system
124: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125: KSPSolve(ksp, b, x);
126: KSPGetIterationNumber(ksp, &its);
128: /* Check error */
129: MatMult(A, x, u);
130: VecAXPY(u, -1.0, b);
131: VecNorm(u, NORM_2, &norm);
132: PetscPrintf(PETSC_COMM_WORLD, "Number of iterations = %3" PetscInt_FMT "\n", its);
133: PetscPrintf(PETSC_COMM_WORLD, "Residual norm %g\n", (double)norm);
134: flg = PETSC_FALSE;
135: PetscOptionsGetBool(NULL, NULL, "-ksp_reason", &flg, NULL);
136: if (flg) {
137: KSPConvergedReason reason;
138: KSPGetConvergedReason(ksp, &reason);
139: PetscPrintf(PETSC_COMM_WORLD, "KSPConvergedReason: %s\n", KSPConvergedReasons[reason]);
140: }
142: /* Free work space.*/
143: MatDestroy(&A);
144: VecDestroy(&b);
145: VecDestroy(&u);
146: VecDestroy(&x);
147: KSPDestroy(&ksp);
149: 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*/