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