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