Actual source code: ex73.c
petsc-3.13.6 2020-09-29
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: PetscCalloc2(2*m,&nzd,2*m,&nzo);
150: for (i=0; i<m; i++) {
151: MatGetRow(B,i+rstart,&nzl,&cols,NULL);
152: for (j=0; j<nzl; j++) {
153: if (cols[j] >= rstart && cols[j] < rstart+n) {
154: nzd[2*i] += 2;
155: nzd[2*i+1] += 2;
156: } else {
157: nzo[2*i] += 2;
158: nzo[2*i+1] += 2;
159: }
160: }
161: nzmax = PetscMax(nzmax,nzd[2*i]+nzo[2*i]);
162: MatRestoreRow(B,i+rstart,&nzl,&cols,NULL);
163: }
164: MatCreateAIJ(PETSC_COMM_WORLD,2*m,2*m,PETSC_DECIDE,PETSC_DECIDE,0,nzd,0,nzo,&J);
165: PetscInfo(0,"Created empty Jacobian matrix\n");
166: PetscFree2(nzd,nzo);
167: PetscMalloc2(nzmax,&ncols,nzmax,&nvals);
168: PetscArrayzero(nvals,nzmax);
169: for (i=0; i<m; i++) {
170: MatGetRow(B,i+rstart,&nzl,&cols,&vals);
171: for (j=0; j<nzl; j++) {
172: ncols[2*j] = 2*cols[j];
173: ncols[2*j+1] = 2*cols[j]+1;
174: }
175: nrow = 2*(i+rstart);
176: MatSetValues(J,1,&nrow,2*nzl,ncols,nvals,INSERT_VALUES);
177: nrow = 2*(i+rstart) + 1;
178: MatSetValues(J,1,&nrow,2*nzl,ncols,nvals,INSERT_VALUES);
179: MatRestoreRow(B,i+rstart,&nzl,&cols,&vals);
180: }
181: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
182: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
183: if (viewMats) {
184: PetscPrintf(PETSC_COMM_WORLD,"Jacobian matrix structure:\n");
185: MatView(J,PETSC_VIEWER_DRAW_WORLD);
186: }
187: MatDestroy(&J);
188: PetscFree2(ncols,nvals);
189: }
191: /*
192: Free work space. All PETSc objects should be destroyed when they
193: are no longer needed.
194: */
195: MatDestroy(&B);
196: VecDestroy(&xin);
197: PetscFinalize();
198: return ierr;
199: }
201: /*TEST
203: test:
204: nsize: 3
205: requires: parmetis datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
206: args: -nox -f ${DATAFILESPATH}/matrices/arco1 -mat_partitioning_type parmetis -viewer_binary_skip_info -novec_load
208: test:
209: requires: parmetis !complex double !define(PETSC_USE_64BIT_INDICES)
210: output_file: output/ex73_1.out
211: suffix: parmetis_nd_32
212: nsize: 3
213: 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
215: test:
216: requires: parmetis !complex double define(PETSC_USE_64BIT_INDICES)
217: output_file: output/ex73_1.out
218: suffix: parmetis_nd_64
219: nsize: 3
220: 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
222: test:
223: requires: ptscotch !complex double !define(PETSC_USE_64BIT_INDICES) define(PETSC_HAVE_SCOTCH_PARMETIS_V3_NODEND)
224: output_file: output/ex73_1.out
225: suffix: ptscotch_nd_32
226: nsize: 4
227: 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
229: test:
230: requires: ptscotch !complex double define(PETSC_USE_64BIT_INDICES) define(PETSC_HAVE_SCOTCH_PARMETIS_V3_NODEND)
231: output_file: output/ex73_1.out
232: suffix: ptscotch_nd_64
233: nsize: 4
234: 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
236: TEST*/