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;
22: PetscInt its,m,n;
23: PetscReal norm;
24: PetscMPIInt size,rank;
25: PetscScalar one = 1.0;
27: PetscInitialize(&argc,&args,(char*)0,help);
28: MPI_Comm_size(PETSC_COMM_WORLD,&size);
29: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
31: PetscOptionsGetBool(NULL,NULL,"-partition",&partition,NULL);
32: PetscOptionsGetBool(NULL,NULL,"-displayIS",&displayIS,NULL);
33: PetscOptionsGetBool(NULL,NULL,"-displayMat",&displayMat,NULL);
35: /* Determine file from which we read the matrix.*/
36: PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg);
39: /* - - - - - - - - - - - - - - - - - - - - - - - -
40: Load system
41: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
42: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
43: MatCreate(PETSC_COMM_WORLD,&A);
44: MatLoad(A,fd);
45: PetscViewerDestroy(&fd);
46: MatGetLocalSize(A,&m,&n);
49: /* Create rhs vector of all ones */
50: VecCreate(PETSC_COMM_WORLD,&b);
51: VecSetSizes(b,m,PETSC_DECIDE);
52: VecSetFromOptions(b);
53: VecSet(b,one);
55: VecDuplicate(b,&x);
56: VecDuplicate(b,&u);
57: VecSet(x,0.0);
59: /* - - - - - - - - - - - - - - - - - - - - - - - -
60: Test partition
61: - - - - - - - - - - - - - - - - - - - - - - - - - */
62: if (partition) {
63: MatPartitioning mpart;
64: IS mis,nis,is;
65: PetscInt *count;
66: Mat BB;
68: if (displayMat) {
69: PetscPrintf(PETSC_COMM_WORLD,"Before partitioning/reordering, A:\n");
70: MatView(A,PETSC_VIEWER_DRAW_WORLD);
71: }
73: PetscMalloc1(size,&count);
74: MatPartitioningCreate(PETSC_COMM_WORLD, &mpart);
75: MatPartitioningSetAdjacency(mpart, A);
76: /* MatPartitioningSetVertexWeights(mpart, weight); */
77: MatPartitioningSetFromOptions(mpart);
78: MatPartitioningApply(mpart, &mis);
79: MatPartitioningDestroy(&mpart);
80: if (displayIS) {
81: PetscPrintf(PETSC_COMM_WORLD,"mis, new processor assignment:\n");
82: ISView(mis,PETSC_VIEWER_STDOUT_WORLD);
83: }
85: ISPartitioningToNumbering(mis,&nis);
86: if (displayIS) {
87: PetscPrintf(PETSC_COMM_WORLD,"nis:\n");
88: ISView(nis,PETSC_VIEWER_STDOUT_WORLD);
89: }
91: ISPartitioningCount(mis,size,count);
92: ISDestroy(&mis);
93: if (displayIS && rank == 0) {
94: PetscInt i;
95: PetscPrintf(PETSC_COMM_SELF,"[ %d ] count:\n",rank);
96: for (i=0; i<size; i++) PetscPrintf(PETSC_COMM_WORLD," %d",count[i]);
97: PetscPrintf(PETSC_COMM_WORLD,"\n");
98: }
100: ISInvertPermutation(nis, count[rank], &is);
101: PetscFree(count);
102: ISDestroy(&nis);
103: ISSort(is);
104: if (displayIS) {
105: PetscPrintf(PETSC_COMM_WORLD,"inverse of nis - maps new local rows to old global rows:\n");
106: ISView(is,PETSC_VIEWER_STDOUT_WORLD);
107: }
109: MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&BB);
110: if (displayMat) {
111: PetscPrintf(PETSC_COMM_WORLD,"After partitioning/reordering, A:\n");
112: MatView(BB,PETSC_VIEWER_DRAW_WORLD);
113: }
115: /* need to move the vector also */
116: ISDestroy(&is);
117: MatDestroy(&A);
118: A = BB;
119: }
121: /* Create linear solver; set operators; set runtime options.*/
122: KSPCreate(PETSC_COMM_WORLD,&ksp);
123: KSPSetOperators(ksp,A,A);
124: KSPSetFromOptions(ksp);
126: /* - - - - - - - - - - - - - - - - - - - - - - - -
127: Solve system
128: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
129: KSPSolve(ksp,b,x);
130: KSPGetIterationNumber(ksp,&its);
132: /* Check error */
133: MatMult(A,x,u);
134: VecAXPY(u,-1.0,b);
135: VecNorm(u,NORM_2,&norm);
136: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
137: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);
138: flg = PETSC_FALSE;
139: PetscOptionsGetBool(NULL,NULL, "-ksp_reason", &flg,NULL);
140: if (flg) {
141: KSPConvergedReason reason;
142: KSPGetConvergedReason(ksp,&reason);
143: PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %D\n", reason);
144: }
146: /* Free work space.*/
147: MatDestroy(&A)); PetscCall(VecDestroy(&b);
148: VecDestroy(&u)); PetscCall(VecDestroy(&x);
149: KSPDestroy(&ksp);
151: PetscFinalize();
152: return 0;
153: }
155: /*TEST
157: test:
158: args: -f ${DATAFILESPATH}/matrices/small -partition -mat_partitioning_type parmetis
159: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) parmetis
160: output_file: output/ex31.out
161: nsize: 3
163: TEST*/