Actual source code: ex73.c

petsc-3.11.4 2019-09-28
Report Typos and Errors

  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*/