Actual source code: ex73.c
petsc-3.10.5 2019-03-28
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*/
11: /*
12: Include "petscmat.h" so that we can use matrices. Note that this file
13: automatically includes:
14: petscsys.h - base PETSc routines petscvec.h - vectors
15: petscmat.h - matrices
16: petscis.h - index sets
17: petscviewer.h - viewers
19: Example of usage:
20: mpiexec -n 3 ex73 -f <matfile> -mat_partitioning_type parmetis/scotch -viewer_binary_skip_info -nox
21: */
22: #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,useND,noVecLoad = PETSC_FALSE;
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);if (ierr) return ierr;
39: MPI_Comm_size(PETSC_COMM_WORLD,&size);
40: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
41: PetscOptionsHasName(NULL,NULL, "-view_mats", &viewMats);
42: PetscOptionsHasName(NULL,NULL, "-view_is", &viewIS);
43: PetscOptionsHasName(NULL,NULL, "-view_vecs", &viewVecs);
44: PetscOptionsHasName(NULL,NULL, "-use_nd", &useND);
45: PetscOptionsHasName(NULL,NULL, "-novec_load", &noVecLoad);
47: /*
48: Determine file from which we read the matrix
49: */
50: PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);
52: /*
53: Open binary file. Note that we use FILE_MODE_READ to indicate
54: reading from this file.
55: */
56: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
58: /*
59: Load the matrix and vector; then destroy the viewer.
60: */
61: MatCreate(PETSC_COMM_WORLD,&A);
62: MatSetType(A,mtype);
63: MatLoad(A,fd);
64: if (!noVecLoad) {
65: VecCreate(PETSC_COMM_WORLD,&xin);
66: VecLoad(xin,fd);
67: } else {
68: MatCreateVecs(A,&xin,NULL);
69: VecSetRandom(xin,NULL);
70: }
71: PetscViewerDestroy(&fd);
72: if (viewMats) {
73: PetscPrintf(PETSC_COMM_WORLD,"Original matrix:\n");
74: MatView(A,PETSC_VIEWER_DRAW_WORLD);
75: }
76: if (viewVecs) {
77: PetscPrintf(PETSC_COMM_WORLD,"Original vector:\n");
78: VecView(xin,PETSC_VIEWER_STDOUT_WORLD);
79: }
81: /* Partition the graph of the matrix */
82: MatPartitioningCreate(PETSC_COMM_WORLD,&part);
83: MatPartitioningSetAdjacency(part,A);
84: MatPartitioningSetFromOptions(part);
86: /* get new processor owner number of each vertex */
87: if (useND) {
88: MatPartitioningApplyND(part,&is);
89: } else {
90: MatPartitioningApply(part,&is);
91: }
92: if (viewIS) {
93: PetscPrintf(PETSC_COMM_WORLD,"IS1 - new processor ownership:\n");
94: ISView(is,PETSC_VIEWER_STDOUT_WORLD);
95: }
97: /* get new global number of each old global number */
98: ISPartitioningToNumbering(is,&isn);
99: if (viewIS) {
100: PetscPrintf(PETSC_COMM_WORLD,"IS2 - new global numbering:\n");
101: ISView(isn,PETSC_VIEWER_STDOUT_WORLD);
102: }
104: /* get number of new vertices for each processor */
105: PetscMalloc1(size,&nlocal);
106: ISPartitioningCount(is,size,nlocal);
107: ISDestroy(&is);
109: /* get old global number of each new global number */
110: ISInvertPermutation(isn,useND ? PETSC_DECIDE : nlocal[rank],&is);
111: if (viewIS) {
112: PetscPrintf(PETSC_COMM_WORLD,"IS3=inv(IS2) - old global number of each new global number:\n");
113: ISView(is,PETSC_VIEWER_STDOUT_WORLD);
114: }
116: /* move the matrix rows to the new processes they have been assigned to by the permutation */
117: MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&B);
118: PetscFree(nlocal);
119: ISDestroy(&isn);
120: MatDestroy(&A);
121: MatPartitioningDestroy(&part);
122: if (viewMats) {
123: PetscPrintf(PETSC_COMM_WORLD,"Partitioned matrix:\n");
124: MatView(B,PETSC_VIEWER_DRAW_WORLD);
125: }
127: /* move the vector rows to the new processes they have been assigned to */
128: MatGetLocalSize(B,&m,&n);
129: VecCreateMPI(PETSC_COMM_WORLD,m,PETSC_DECIDE,&xout);
130: VecScatterCreate(xin,is,xout,NULL,&scat);
131: VecScatterBegin(scat,xin,xout,INSERT_VALUES,SCATTER_FORWARD);
132: VecScatterEnd(scat,xin,xout,INSERT_VALUES,SCATTER_FORWARD);
133: VecScatterDestroy(&scat);
134: if (viewVecs) {
135: PetscPrintf(PETSC_COMM_WORLD,"Mapped vector:\n");
136: VecView(xout,PETSC_VIEWER_STDOUT_WORLD);
137: }
138: VecDestroy(&xout);
139: ISDestroy(&is);
141: {
142: PetscInt rstart,i,*nzd,*nzo,nzl,nzmax = 0,*ncols,nrow,j;
143: Mat J;
144: const PetscInt *cols;
145: const PetscScalar *vals;
146: PetscScalar *nvals;
148: MatGetOwnershipRange(B,&rstart,NULL);
149: PetscCalloc1(2*m,&nzd);
150: PetscCalloc1(2*m,&nzo);
151: for (i=0; i<m; i++) {
152: MatGetRow(B,i+rstart,&nzl,&cols,NULL);
153: for (j=0; j<nzl; j++) {
154: if (cols[j] >= rstart && cols[j] < rstart+n) {
155: nzd[2*i] += 2;
156: nzd[2*i+1] += 2;
157: } else {
158: nzo[2*i] += 2;
159: nzo[2*i+1] += 2;
160: }
161: }
162: nzmax = PetscMax(nzmax,nzd[2*i]+nzo[2*i]);
163: MatRestoreRow(B,i+rstart,&nzl,&cols,NULL);
164: }
165: MatCreateAIJ(PETSC_COMM_WORLD,2*m,2*m,PETSC_DECIDE,PETSC_DECIDE,0,nzd,0,nzo,&J);
166: PetscInfo(0,"Created empty Jacobian matrix\n");
167: PetscFree(nzd);
168: PetscFree(nzo);
169: PetscMalloc2(nzmax,&ncols,nzmax,&nvals);
170: PetscMemzero(nvals,nzmax*sizeof(PetscScalar));
171: for (i=0; i<m; i++) {
172: MatGetRow(B,i+rstart,&nzl,&cols,&vals);
173: for (j=0; j<nzl; j++) {
174: ncols[2*j] = 2*cols[j];
175: ncols[2*j+1] = 2*cols[j]+1;
176: }
177: nrow = 2*(i+rstart);
178: MatSetValues(J,1,&nrow,2*nzl,ncols,nvals,INSERT_VALUES);
179: nrow = 2*(i+rstart) + 1;
180: MatSetValues(J,1,&nrow,2*nzl,ncols,nvals,INSERT_VALUES);
181: MatRestoreRow(B,i+rstart,&nzl,&cols,&vals);
182: }
183: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
184: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
185: if (viewMats) {
186: PetscPrintf(PETSC_COMM_WORLD,"Jacobian matrix structure:\n");
187: MatView(J,PETSC_VIEWER_DRAW_WORLD);
188: }
189: MatDestroy(&J);
190: PetscFree2(ncols,nvals);
191: }
193: /*
194: Free work space. All PETSc objects should be destroyed when they
195: are no longer needed.
196: */
197: MatDestroy(&B);
198: VecDestroy(&xin);
199: PetscFinalize();
200: return ierr;
201: }
203: /*TEST
205: test:
206: nsize: 3
207: requires: parmetis datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
208: args: -nox -f ${DATAFILESPATH}/matrices/arco1 -mat_partitioning_type parmetis -viewer_binary_skip_info -novec_load
210: test:
211: requires: parmetis !complex double !define(PETSC_USE_64BIT_INDICES)
212: output_file: output/ex73_1.out
213: suffix: parmetis_nd_32
214: nsize: 3
215: args: -nox -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/spd-real-int32-float64 -mat_partitioning_type parmetis -viewer_binary_skip_info -use_nd -novec_load
217: test:
218: requires: parmetis !complex double define(PETSC_USE_64BIT_INDICES)
219: output_file: output/ex73_1.out
220: suffix: parmetis_nd_64
221: nsize: 3
222: args: -nox -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/spd-real-int64-float64 -mat_partitioning_type parmetis -viewer_binary_skip_info -use_nd -novec_load
224: test:
225: requires: ptscotch !complex double !define(PETSC_USE_64BIT_INDICES) define(PETSC_HAVE_SCOTCH_PARMETIS_V3_NODEND)
226: output_file: output/ex73_1.out
227: suffix: ptscotch_nd_32
228: nsize: 4
229: args: -nox -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/spd-real-int32-float64 -mat_partitioning_type ptscotch -viewer_binary_skip_info -use_nd -novec_load
231: test:
232: requires: ptscotch !complex double define(PETSC_USE_64BIT_INDICES) define(PETSC_HAVE_SCOTCH_PARMETIS_V3_NODEND)
233: output_file: output/ex73_1.out
234: suffix: ptscotch_nd_64
235: nsize: 4
236: args: -nox -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/spd-real-int64-float64 -mat_partitioning_type ptscotch -viewer_binary_skip_info -use_nd -novec_load
238: TEST*/