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