Actual source code: ex73.c

petsc-3.12.5 2020-03-29
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:     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*/