Actual source code: ex73.c
2: static char help[] = "Reads a PETSc matrix from a file partitions it\n\n";
4: /*T
5: Concepts: partitioning
6: Processors: n
7: T*/
9: /*
10: Include "petscmat.h" so that we can use matrices. Note that this file
11: automatically includes:
12: petscsys.h - base PETSc routines petscvec.h - vectors
13: petscmat.h - matrices
14: petscis.h - index sets
15: petscviewer.h - viewers
17: Example of usage:
18: mpiexec -n 3 ex73 -f <matfile> -mat_partitioning_type parmetis/scotch -viewer_binary_skip_info -nox
19: */
20: #include <petscksp.h>
24: int main(int argc,char **args)
25: {
26: const MatType mtype = MATMPIAIJ; /* matrix format */
27: Mat A,B; /* matrix */
28: PetscViewer fd; /* viewer */
29: char file[PETSC_MAX_PATH_LEN]; /* input file name */
30: PetscBool flg,viewMats,viewIS;
31: PetscInt ierr,*nlocal,m,n;
32: PetscMPIInt rank,size;
33: MatPartitioning part;
34: IS is,isn;
35: Vec xin, xout;
36: VecScatter scat;
38: PetscInitialize(&argc,&args,(char *)0,help);
39: MPI_Comm_size(PETSC_COMM_WORLD,&size);
40: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
41: PetscOptionsHasName(PETSC_NULL, "-view_mats", &viewMats);
42: PetscOptionsHasName(PETSC_NULL, "-view_is", &viewIS);
44: /*
45: Determine file from which we read the matrix
46: */
47: PetscOptionsGetString(PETSC_NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);
49: /*
50: Open binary file. Note that we use FILE_MODE_READ to indicate
51: reading from this file.
52: */
53: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
55: /*
56: Load the matrix and vector; then destroy the viewer.
57: */
58: MatCreate(PETSC_COMM_WORLD,&A);
59: MatSetType(A,mtype);
60: MatLoad(A,fd);
61: VecCreate(PETSC_COMM_WORLD,&xin);
62: VecLoad(xin,fd);
63: PetscViewerDestroy(&fd);
64: if (viewMats){
65: if (!rank) printf("Original matrix:\n");
66: MatView(A,PETSC_VIEWER_DRAW_WORLD);
67: }
69: /* Partition the graph of the matrix */
70: MatPartitioningCreate(PETSC_COMM_WORLD,&part);
71: MatPartitioningSetAdjacency(part,A);
72: MatPartitioningSetFromOptions(part);
74: /* get new processor owner number of each vertex */
75: MatPartitioningApply(part,&is);
76: if (viewIS){
77: if (!rank) printf("IS1 - new processor ownership:\n");
78: ISView(is,PETSC_VIEWER_STDOUT_WORLD);
79: }
81: /* get new global number of each old global number */
82: ISPartitioningToNumbering(is,&isn);
83: if (viewIS){
84: if (!rank) printf("IS2 - new global numbering:\n");
85: ISView(isn,PETSC_VIEWER_STDOUT_WORLD);
86: }
88: /* get number of new vertices for each processor */
89: PetscMalloc(size*sizeof(PetscInt),&nlocal);
90: ISPartitioningCount(is,size,nlocal);
91: ISDestroy(&is);
93: /* get old global number of each new global number */
94: ISInvertPermutation(isn,nlocal[rank],&is);
95: PetscFree(nlocal);
96: ISDestroy(&isn);
97: MatPartitioningDestroy(&part);
98: if (viewIS){
99: if (!rank) printf("IS3=inv(IS2) - old global number of each new global number:\n");
100: ISView(is,PETSC_VIEWER_STDOUT_WORLD);
101: }
103: /* move the matrix rows to the new processes they have been assigned to by the permutation */
104: ISSort(is);
105: MatGetSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&B);
106: MatDestroy(&A);
108: /* move the vector rows to the new processes they have been assigned to */
109: MatGetLocalSize(B,&m,&n);
110: VecCreateMPI(PETSC_COMM_WORLD,m,PETSC_DECIDE,&xout);
111: VecScatterCreate(xin,PETSC_NULL,xout,is,&scat);
112: VecScatterBegin(scat,xin,xout,INSERT_VALUES,SCATTER_FORWARD);
113: VecScatterEnd(scat,xin,xout,INSERT_VALUES,SCATTER_FORWARD);
114: VecScatterDestroy(&scat);
115: ISDestroy(&is);
116: if (viewMats){
117: if (!rank) printf("Partitioned matrix:\n");
118: MatView(B,PETSC_VIEWER_DRAW_WORLD);
119: }
121: {
122: PetscInt rstart,i,*nzd,*nzo,nzl,nzmax = 0,*ncols,nrow,j;
123: Mat J;
124: const PetscInt *cols;
125: const PetscScalar *vals;
126: PetscScalar *nvals;
127:
128: MatGetOwnershipRange(B,&rstart,PETSC_NULL);
129: PetscMalloc(2*m*sizeof(PetscInt),&nzd);
130: PetscMemzero(nzd,2*m*sizeof(PetscInt));
131: PetscMalloc(2*m*sizeof(PetscInt),&nzo);
132: PetscMemzero(nzo,2*m*sizeof(PetscInt));
133: for (i=0; i<m; i++) {
134: MatGetRow(B,i+rstart,&nzl,&cols,PETSC_NULL);
135: for (j=0; j<nzl; j++) {
136: if (cols[j] >= rstart && cols[j] < rstart+n) {nzd[2*i] += 2; nzd[2*i+1] += 2;}
137: else {nzo[2*i] += 2; nzo[2*i+1] += 2;}
138: }
139: nzmax = PetscMax(nzmax,nzd[2*i]+nzo[2*i]);
140: MatRestoreRow(B,i+rstart,&nzl,&cols,PETSC_NULL);
141: }
142: MatCreateMPIAIJ(PETSC_COMM_WORLD,2*m,2*m,PETSC_DECIDE,PETSC_DECIDE,0,nzd,0,nzo,&J);
143: PetscInfo(0,"Created empty Jacobian matrix\n");
144: PetscFree(nzd);
145: PetscFree(nzo);
146: PetscMalloc2(nzmax,PetscInt,&ncols,nzmax,PetscScalar,&nvals);
147: PetscMemzero(nvals,nzmax*sizeof(PetscScalar));
148: for (i=0; i<m; i++) {
149: MatGetRow(B,i+rstart,&nzl,&cols,&vals);
150: for (j=0; j<nzl; j++) {
151: ncols[2*j] = 2*cols[j];
152: ncols[2*j+1] = 2*cols[j]+1;
153: }
154: nrow = 2*(i+rstart);
155: MatSetValues(J,1,&nrow,2*nzl,ncols,nvals,INSERT_VALUES);
156: nrow = 2*(i+rstart) + 1;
157: MatSetValues(J,1,&nrow,2*nzl,ncols,nvals,INSERT_VALUES);
158: MatRestoreRow(B,i+rstart,&nzl,&cols,&vals);
159: }
160: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
161: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
162: if (viewMats){
163: if (!rank) printf("Jacobian matrix structure:\n");
164: MatView(J,PETSC_VIEWER_DRAW_WORLD);
165: }
166: MatDestroy(&J);
167: PetscFree2(ncols,nvals);
168: }
169:
170: /*
171: Free work space. All PETSc objects should be destroyed when they
172: are no longer needed.
173: */
174: MatDestroy(&B);
175: VecDestroy(&xin);
176: VecDestroy(&xout);
177: PetscFinalize();
178: return 0;
179: }