Actual source code: ex73.c
petsc-3.6.1 2015-08-06
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 <petscmat.h>
24: int main(int argc,char **args)
25: {
26: 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,viewVecs;
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(NULL, "-view_mats", &viewMats);
42: PetscOptionsHasName(NULL, "-view_is", &viewIS);
43: PetscOptionsHasName(NULL, "-view_vecs", &viewVecs);
45: /*
46: Determine file from which we read the matrix
47: */
48: PetscOptionsGetString(NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);
50: /*
51: Open binary file. Note that we use FILE_MODE_READ to indicate
52: reading from this file.
53: */
54: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
56: /*
57: Load the matrix and vector; then destroy the viewer.
58: */
59: MatCreate(PETSC_COMM_WORLD,&A);
60: MatSetType(A,mtype);
61: MatLoad(A,fd);
62: VecCreate(PETSC_COMM_WORLD,&xin);
63: VecLoad(xin,fd);
64: PetscViewerDestroy(&fd);
65: if (viewMats) {
66: if (!rank) printf("Original matrix:\n");
67: MatView(A,PETSC_VIEWER_DRAW_WORLD);
68: }
69: if (viewVecs) {
70: if (!rank) printf("Original vector:\n");
71: VecView(xin,PETSC_VIEWER_STDOUT_WORLD);
72: }
74: /* Partition the graph of the matrix */
75: MatPartitioningCreate(PETSC_COMM_WORLD,&part);
76: MatPartitioningSetAdjacency(part,A);
77: MatPartitioningSetFromOptions(part);
79: /* get new processor owner number of each vertex */
80: MatPartitioningApply(part,&is);
81: if (viewIS) {
82: if (!rank) printf("IS1 - new processor ownership:\n");
83: ISView(is,PETSC_VIEWER_STDOUT_WORLD);
84: }
86: /* get new global number of each old global number */
87: ISPartitioningToNumbering(is,&isn);
88: if (viewIS) {
89: if (!rank) printf("IS2 - new global numbering:\n");
90: ISView(isn,PETSC_VIEWER_STDOUT_WORLD);
91: }
93: /* get number of new vertices for each processor */
94: PetscMalloc1(size,&nlocal);
95: ISPartitioningCount(is,size,nlocal);
96: ISDestroy(&is);
98: /* get old global number of each new global number */
99: ISInvertPermutation(isn,nlocal[rank],&is);
100: PetscFree(nlocal);
101: ISDestroy(&isn);
102: MatPartitioningDestroy(&part);
103: if (viewIS) {
104: if (!rank) printf("IS3=inv(IS2) - old global number of each new global number:\n");
105: ISView(is,PETSC_VIEWER_STDOUT_WORLD);
106: }
108: /* move the matrix rows to the new processes they have been assigned to by the permutation */
109: ISSort(is);
110: MatGetSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&B);
111: MatDestroy(&A);
113: /* move the vector rows to the new processes they have been assigned to */
114: MatGetLocalSize(B,&m,&n);
115: VecCreateMPI(PETSC_COMM_WORLD,m,PETSC_DECIDE,&xout);
116: VecScatterCreate(xin,is,xout,NULL,&scat);
117: VecScatterBegin(scat,xin,xout,INSERT_VALUES,SCATTER_FORWARD);
118: VecScatterEnd(scat,xin,xout,INSERT_VALUES,SCATTER_FORWARD);
119: VecScatterDestroy(&scat);
120: ISDestroy(&is);
121: if (viewMats) {
122: if (!rank) printf("Partitioned matrix:\n");
123: MatView(B,PETSC_VIEWER_DRAW_WORLD);
124: }
125: if (viewVecs) {
126: if (!rank) printf("Mapped vector:\n");
127: VecView(xout,PETSC_VIEWER_STDOUT_WORLD);
128: }
130: {
131: PetscInt rstart,i,*nzd,*nzo,nzl,nzmax = 0,*ncols,nrow,j;
132: Mat J;
133: const PetscInt *cols;
134: const PetscScalar *vals;
135: PetscScalar *nvals;
137: MatGetOwnershipRange(B,&rstart,NULL);
138: PetscCalloc1(2*m,&nzd);
139: PetscCalloc1(2*m,&nzo);
140: for (i=0; i<m; i++) {
141: MatGetRow(B,i+rstart,&nzl,&cols,NULL);
142: for (j=0; j<nzl; j++) {
143: if (cols[j] >= rstart && cols[j] < rstart+n) {
144: nzd[2*i] += 2;
145: nzd[2*i+1] += 2;
146: } else {
147: nzo[2*i] += 2;
148: nzo[2*i+1] += 2;
149: }
150: }
151: nzmax = PetscMax(nzmax,nzd[2*i]+nzo[2*i]);
152: MatRestoreRow(B,i+rstart,&nzl,&cols,NULL);
153: }
154: MatCreateAIJ(PETSC_COMM_WORLD,2*m,2*m,PETSC_DECIDE,PETSC_DECIDE,0,nzd,0,nzo,&J);
155: PetscInfo(0,"Created empty Jacobian matrix\n");
156: PetscFree(nzd);
157: PetscFree(nzo);
158: PetscMalloc2(nzmax,&ncols,nzmax,&nvals);
159: PetscMemzero(nvals,nzmax*sizeof(PetscScalar));
160: for (i=0; i<m; i++) {
161: MatGetRow(B,i+rstart,&nzl,&cols,&vals);
162: for (j=0; j<nzl; j++) {
163: ncols[2*j] = 2*cols[j];
164: ncols[2*j+1] = 2*cols[j]+1;
165: }
166: nrow = 2*(i+rstart);
167: MatSetValues(J,1,&nrow,2*nzl,ncols,nvals,INSERT_VALUES);
168: nrow = 2*(i+rstart) + 1;
169: MatSetValues(J,1,&nrow,2*nzl,ncols,nvals,INSERT_VALUES);
170: MatRestoreRow(B,i+rstart,&nzl,&cols,&vals);
171: }
172: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
173: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
174: if (viewMats) {
175: if (!rank) printf("Jacobian matrix structure:\n");
176: MatView(J,PETSC_VIEWER_DRAW_WORLD);
177: }
178: MatDestroy(&J);
179: PetscFree2(ncols,nvals);
180: }
182: /*
183: Free work space. All PETSc objects should be destroyed when they
184: are no longer needed.
185: */
186: MatDestroy(&B);
187: VecDestroy(&xin);
188: VecDestroy(&xout);
189: PetscFinalize();
190: return 0;
191: }