Actual source code: ex73.c
petsc-3.12.4 2020-02-04
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*/