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