Actual source code: ex73.c

petsc-3.7.3 2016-08-01
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*/

  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>

 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;
 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);
 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);

 45:   /*
 46:      Determine file from which we read the matrix
 47:   */
 48:   PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);

 50:   /*
 51:        Open binary file.  Note that we use FILE_MODE_READ to indicate
 52:        reading from this file.
 53:   */
 54:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);

 56:   /*
 57:       Load the matrix and vector; then destroy the viewer.
 58:   */
 59:   MatCreate(PETSC_COMM_WORLD,&A);
 60:   MatSetType(A,mtype);
 61:   MatLoad(A,fd);
 62:   VecCreate(PETSC_COMM_WORLD,&xin);
 63:   VecLoad(xin,fd);
 64:   PetscViewerDestroy(&fd);
 65:   if (viewMats) {
 66:     if (!rank) printf("Original matrix:\n");
 67:     MatView(A,PETSC_VIEWER_DRAW_WORLD);
 68:   }
 69:   if (viewVecs) {
 70:     if (!rank) printf("Original vector:\n");
 71:     VecView(xin,PETSC_VIEWER_STDOUT_WORLD);
 72:   }

 74:   /* Partition the graph of the matrix */
 75:   MatPartitioningCreate(PETSC_COMM_WORLD,&part);
 76:   MatPartitioningSetAdjacency(part,A);
 77:   MatPartitioningSetFromOptions(part);

 79:   /* get new processor owner number of each vertex */
 80:   MatPartitioningApply(part,&is);
 81:   if (viewIS) {
 82:     if (!rank) printf("IS1 - new processor ownership:\n");
 83:     ISView(is,PETSC_VIEWER_STDOUT_WORLD);
 84:   }

 86:   /* get new global number of each old global number */
 87:   ISPartitioningToNumbering(is,&isn);
 88:   if (viewIS) {
 89:     if (!rank) printf("IS2 - new global numbering:\n");
 90:     ISView(isn,PETSC_VIEWER_STDOUT_WORLD);
 91:   }

 93:   /* get number of new vertices for each processor */
 94:   PetscMalloc1(size,&nlocal);
 95:   ISPartitioningCount(is,size,nlocal);
 96:   ISDestroy(&is);

 98:   /* get old global number of each new global number */
 99:   ISInvertPermutation(isn,nlocal[rank],&is);
100:   PetscFree(nlocal);
101:   ISDestroy(&isn);
102:   MatPartitioningDestroy(&part);
103:   if (viewIS) {
104:     if (!rank) printf("IS3=inv(IS2) - old global number of each new global number:\n");
105:     ISView(is,PETSC_VIEWER_STDOUT_WORLD);
106:   }

108:   /* move the matrix rows to the new processes they have been assigned to by the permutation */
109:   ISSort(is);
110:   MatGetSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&B);
111:   MatDestroy(&A);

113:   /* move the vector rows to the new processes they have been assigned to */
114:   MatGetLocalSize(B,&m,&n);
115:   VecCreateMPI(PETSC_COMM_WORLD,m,PETSC_DECIDE,&xout);
116:   VecScatterCreate(xin,is,xout,NULL,&scat);
117:   VecScatterBegin(scat,xin,xout,INSERT_VALUES,SCATTER_FORWARD);
118:   VecScatterEnd(scat,xin,xout,INSERT_VALUES,SCATTER_FORWARD);
119:   VecScatterDestroy(&scat);
120:   ISDestroy(&is);
121:   if (viewMats) {
122:     if (!rank) printf("Partitioned matrix:\n");
123:     MatView(B,PETSC_VIEWER_DRAW_WORLD);
124:   }
125:   if (viewVecs) {
126:     if (!rank) printf("Mapped vector:\n");
127:     VecView(xout,PETSC_VIEWER_STDOUT_WORLD);
128:   }

130:   {
131:     PetscInt          rstart,i,*nzd,*nzo,nzl,nzmax = 0,*ncols,nrow,j;
132:     Mat               J;
133:     const PetscInt    *cols;
134:     const PetscScalar *vals;
135:     PetscScalar       *nvals;

137:     MatGetOwnershipRange(B,&rstart,NULL);
138:     PetscCalloc1(2*m,&nzd);
139:     PetscCalloc1(2*m,&nzo);
140:     for (i=0; i<m; i++) {
141:       MatGetRow(B,i+rstart,&nzl,&cols,NULL);
142:       for (j=0; j<nzl; j++) {
143:         if (cols[j] >= rstart && cols[j] < rstart+n) {
144:           nzd[2*i] += 2;
145:           nzd[2*i+1] += 2;
146:         } else {
147:           nzo[2*i] += 2;
148:           nzo[2*i+1] += 2;
149:         }
150:       }
151:       nzmax = PetscMax(nzmax,nzd[2*i]+nzo[2*i]);
152:       MatRestoreRow(B,i+rstart,&nzl,&cols,NULL);
153:     }
154:     MatCreateAIJ(PETSC_COMM_WORLD,2*m,2*m,PETSC_DECIDE,PETSC_DECIDE,0,nzd,0,nzo,&J);
155:     PetscInfo(0,"Created empty Jacobian matrix\n");
156:     PetscFree(nzd);
157:     PetscFree(nzo);
158:     PetscMalloc2(nzmax,&ncols,nzmax,&nvals);
159:     PetscMemzero(nvals,nzmax*sizeof(PetscScalar));
160:     for (i=0; i<m; i++) {
161:       MatGetRow(B,i+rstart,&nzl,&cols,&vals);
162:       for (j=0; j<nzl; j++) {
163:         ncols[2*j]   = 2*cols[j];
164:         ncols[2*j+1] = 2*cols[j]+1;
165:       }
166:       nrow = 2*(i+rstart);
167:       MatSetValues(J,1,&nrow,2*nzl,ncols,nvals,INSERT_VALUES);
168:       nrow = 2*(i+rstart) + 1;
169:       MatSetValues(J,1,&nrow,2*nzl,ncols,nvals,INSERT_VALUES);
170:       MatRestoreRow(B,i+rstart,&nzl,&cols,&vals);
171:     }
172:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
173:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
174:     if (viewMats) {
175:       if (!rank) printf("Jacobian matrix structure:\n");
176:       MatView(J,PETSC_VIEWER_DRAW_WORLD);
177:     }
178:     MatDestroy(&J);
179:     PetscFree2(ncols,nvals);
180:   }

182:   /*
183:        Free work space.  All PETSc objects should be destroyed when they
184:        are no longer needed.
185:   */
186:   MatDestroy(&B);
187:   VecDestroy(&xin);
188:   VecDestroy(&xout);
189:   PetscFinalize();
190:   return 0;
191: }