Actual source code: ex31.c
petsc-3.3-p6 2013-02-11
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>
17: int main(int argc,char **args)
18: {
19: KSP ksp; /* linear solver context */
20: Mat A; /* matrix */
21: Vec x,b,u; /* approx solution, RHS, exact solution */
22: PetscViewer fd; /* viewer */
23: char file[PETSC_MAX_PATH_LEN]; /* input file name */
24: PetscBool flg,partition=PETSC_FALSE,displayIS=PETSC_FALSE,displayMat=PETSC_FALSE;
26: PetscInt its,m,n;
27: PetscReal norm;
28: PetscMPIInt size,rank;
29: PetscScalar one = 1.0;
31: PetscInitialize(&argc,&args,(char *)0,help);
32: MPI_Comm_size(PETSC_COMM_WORLD,&size);
33: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
34:
35: PetscOptionsGetBool(PETSC_NULL,"-partition",&partition,PETSC_NULL);
36: PetscOptionsGetBool(PETSC_NULL,"-displayIS",&displayIS,PETSC_NULL);
37: PetscOptionsGetBool(PETSC_NULL,"-displayMat",&displayMat,PETSC_NULL);
39: /* Determine file from which we read the matrix.*/
40: PetscOptionsGetString(PETSC_NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);
41: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f option");
43: /* - - - - - - - - - - - - - - - - - - - - - - - -
44: Load system
45: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
46: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
47: MatCreate(PETSC_COMM_WORLD,&A);
48: MatLoad(A,fd);
49: PetscViewerDestroy(&fd);
50: MatGetLocalSize(A,&m,&n);
51: if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);
53: /* Create rhs vector of all ones */
54: VecCreate(PETSC_COMM_WORLD,&b);
55: VecSetSizes(b,m,PETSC_DECIDE);
56: VecSetFromOptions(b);
57: VecSet(b,one);
59: VecDuplicate(b,&x);
60: VecDuplicate(b,&u);
61: VecSet(x,0.0);
63: /* - - - - - - - - - - - - - - - - - - - - - - - -
64: Test partition
65: - - - - - - - - - - - - - - - - - - - - - - - - - */
66: if (partition){
67: MatPartitioning mpart;
68: IS mis,nis,is;
69: PetscInt *count;
70: Mat BB;
71:
72: if (displayMat){
73: if (!rank) printf("Before partitioning/reordering, A:\n");
74: MatView(A,PETSC_VIEWER_DRAW_WORLD);
75: }
76:
77: PetscMalloc(size*sizeof(PetscInt),&count);
78: MatPartitioningCreate(PETSC_COMM_WORLD, &mpart);
79: MatPartitioningSetAdjacency(mpart, A);
80: /* MatPartitioningSetVertexWeights(mpart, weight); */
81: MatPartitioningSetFromOptions(mpart);
82: MatPartitioningApply(mpart, &mis);
83: MatPartitioningDestroy(&mpart);
84: if (displayIS){
85: PetscPrintf(PETSC_COMM_WORLD,"mis, new processor assignment:\n");
86: ISView(mis,PETSC_VIEWER_STDOUT_WORLD);
87: }
89: ISPartitioningToNumbering(mis,&nis);
90: if (displayIS){
91: PetscPrintf(PETSC_COMM_WORLD,"nis:\n");
92: ISView(nis,PETSC_VIEWER_STDOUT_WORLD);
93: }
95: ISPartitioningCount(mis,size,count);
96: ISDestroy(&mis);
97: if (displayIS && !rank ) {
98: PetscInt i;
99: printf("[ %d ] count:\n",rank);
100: for (i=0; i<size; i++){
101: printf(" %d",count[i]);
102: }
103: printf("\n");
104: }
105:
106: ISInvertPermutation(nis, count[rank], &is);
107: PetscFree(count);
108: ISDestroy(&nis);
109: ISSort(is);
110: if (displayIS){
111: PetscPrintf(PETSC_COMM_WORLD,"inverse of nis - maps new local rows to old global rows:\n");
112: ISView(is,PETSC_VIEWER_STDOUT_WORLD);
113: }
115: MatGetSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&BB);
116: if (displayMat){
117: if (!rank) printf("After partitioning/reordering, A:\n");
118: MatView(BB,PETSC_VIEWER_DRAW_WORLD);
119: }
121: /* need to move the vector also */
122: ISDestroy(&is);
123: MatDestroy(&A);
124: A = BB;
125: }
127: /* Create linear solver; set operators; set runtime options.*/
128: KSPCreate(PETSC_COMM_WORLD,&ksp);
129: KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);
130: KSPSetFromOptions(ksp);
132: /* - - - - - - - - - - - - - - - - - - - - - - - -
133: Solve system
134: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
135: KSPSolve(ksp,b,x);
136: KSPGetIterationNumber(ksp,&its);
138: /* Check error */
139: MatMult(A,x,u);
140: VecAXPY(u,-1.0,b);
141: VecNorm(u,NORM_2,&norm);
142: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
143: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm);
144: flg = PETSC_FALSE;
145: PetscOptionsGetBool(PETSC_NULL, "-ksp_reason", &flg,PETSC_NULL);
146: if (flg){
147: KSPConvergedReason reason;
148: KSPGetConvergedReason(ksp,&reason);
149: PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %D\n", reason);
150: }
152: /* Free work space.*/
153: MatDestroy(&A); VecDestroy(&b);
154: VecDestroy(&u); VecDestroy(&x);
155: KSPDestroy(&ksp);
157: PetscFinalize();
158: return 0;
159: }