1 2 static char help[] = "Reads a PETSc matrix from a file partitions it\n\n"; 3 4 /*T 5 Concepts: partitioning 6 Processors: n 7 T*/ 8 9 10 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 18 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> 23 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; 37 38 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 39 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 40 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); 41 ierr = PetscOptionsHasName(NULL,NULL, "-view_mats", &viewMats);CHKERRQ(ierr); 42 ierr = PetscOptionsHasName(NULL,NULL, "-view_is", &viewIS);CHKERRQ(ierr); 43 ierr = PetscOptionsHasName(NULL,NULL, "-view_vecs", &viewVecs);CHKERRQ(ierr); 44 ierr = PetscOptionsHasName(NULL,NULL, "-use_nd", &useND);CHKERRQ(ierr); 45 ierr = PetscOptionsHasName(NULL,NULL, "-novec_load", &noVecLoad);CHKERRQ(ierr); 46 47 /* 48 Determine file from which we read the matrix 49 */ 50 ierr = PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 51 52 /* 53 Open binary file. Note that we use FILE_MODE_READ to indicate 54 reading from this file. 55 */ 56 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);CHKERRQ(ierr); 57 58 /* 59 Load the matrix and vector; then destroy the viewer. 60 */ 61 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 62 ierr = MatSetType(A,mtype);CHKERRQ(ierr); 63 ierr = MatLoad(A,fd);CHKERRQ(ierr); 64 if (!noVecLoad) { 65 ierr = VecCreate(PETSC_COMM_WORLD,&xin);CHKERRQ(ierr); 66 ierr = VecLoad(xin,fd);CHKERRQ(ierr); 67 } else { 68 ierr = MatCreateVecs(A,&xin,NULL);CHKERRQ(ierr); 69 ierr = VecSetRandom(xin,NULL);CHKERRQ(ierr); 70 } 71 ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr); 72 if (viewMats) { 73 ierr = PetscPrintf(PETSC_COMM_WORLD,"Original matrix:\n");CHKERRQ(ierr); 74 ierr = MatView(A,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 75 } 76 if (viewVecs) { 77 ierr = PetscPrintf(PETSC_COMM_WORLD,"Original vector:\n");CHKERRQ(ierr); 78 ierr = VecView(xin,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 79 } 80 81 /* Partition the graph of the matrix */ 82 ierr = MatPartitioningCreate(PETSC_COMM_WORLD,&part);CHKERRQ(ierr); 83 ierr = MatPartitioningSetAdjacency(part,A);CHKERRQ(ierr); 84 ierr = MatPartitioningSetFromOptions(part);CHKERRQ(ierr); 85 86 /* get new processor owner number of each vertex */ 87 if (useND) { 88 ierr = MatPartitioningApplyND(part,&is);CHKERRQ(ierr); 89 } else { 90 ierr = MatPartitioningApply(part,&is);CHKERRQ(ierr); 91 } 92 if (viewIS) { 93 ierr = PetscPrintf(PETSC_COMM_WORLD,"IS1 - new processor ownership:\n");CHKERRQ(ierr); 94 ierr = ISView(is,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 95 } 96 97 /* get new global number of each old global number */ 98 ierr = ISPartitioningToNumbering(is,&isn);CHKERRQ(ierr); 99 if (viewIS) { 100 ierr = PetscPrintf(PETSC_COMM_WORLD,"IS2 - new global numbering:\n");CHKERRQ(ierr); 101 ierr = ISView(isn,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 102 } 103 104 /* get number of new vertices for each processor */ 105 ierr = PetscMalloc1(size,&nlocal);CHKERRQ(ierr); 106 ierr = ISPartitioningCount(is,size,nlocal);CHKERRQ(ierr); 107 ierr = ISDestroy(&is);CHKERRQ(ierr); 108 109 /* get old global number of each new global number */ 110 ierr = ISInvertPermutation(isn,useND ? PETSC_DECIDE : nlocal[rank],&is);CHKERRQ(ierr); 111 if (viewIS) { 112 ierr = PetscPrintf(PETSC_COMM_WORLD,"IS3=inv(IS2) - old global number of each new global number:\n");CHKERRQ(ierr); 113 ierr = ISView(is,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 114 } 115 116 /* move the matrix rows to the new processes they have been assigned to by the permutation */ 117 ierr = MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 118 ierr = PetscFree(nlocal);CHKERRQ(ierr); 119 ierr = ISDestroy(&isn);CHKERRQ(ierr); 120 ierr = MatDestroy(&A);CHKERRQ(ierr); 121 ierr = MatPartitioningDestroy(&part);CHKERRQ(ierr); 122 if (viewMats) { 123 ierr = PetscPrintf(PETSC_COMM_WORLD,"Partitioned matrix:\n");CHKERRQ(ierr); 124 ierr = MatView(B,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 125 } 126 127 /* move the vector rows to the new processes they have been assigned to */ 128 ierr = MatGetLocalSize(B,&m,&n);CHKERRQ(ierr); 129 ierr = VecCreateMPI(PETSC_COMM_WORLD,m,PETSC_DECIDE,&xout);CHKERRQ(ierr); 130 ierr = VecScatterCreate(xin,is,xout,NULL,&scat);CHKERRQ(ierr); 131 ierr = VecScatterBegin(scat,xin,xout,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 132 ierr = VecScatterEnd(scat,xin,xout,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 133 ierr = VecScatterDestroy(&scat);CHKERRQ(ierr); 134 if (viewVecs) { 135 ierr = PetscPrintf(PETSC_COMM_WORLD,"Mapped vector:\n");CHKERRQ(ierr); 136 ierr = VecView(xout,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 137 } 138 ierr = VecDestroy(&xout);CHKERRQ(ierr); 139 ierr = ISDestroy(&is);CHKERRQ(ierr); 140 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; 147 148 ierr = MatGetOwnershipRange(B,&rstart,NULL);CHKERRQ(ierr); 149 ierr = PetscCalloc2(2*m,&nzd,2*m,&nzo);CHKERRQ(ierr); 150 for (i=0; i<m; i++) { 151 ierr = MatGetRow(B,i+rstart,&nzl,&cols,NULL);CHKERRQ(ierr); 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 ierr = MatRestoreRow(B,i+rstart,&nzl,&cols,NULL);CHKERRQ(ierr); 163 } 164 ierr = MatCreateAIJ(PETSC_COMM_WORLD,2*m,2*m,PETSC_DECIDE,PETSC_DECIDE,0,nzd,0,nzo,&J);CHKERRQ(ierr); 165 ierr = PetscInfo(0,"Created empty Jacobian matrix\n");CHKERRQ(ierr); 166 ierr = PetscFree2(nzd,nzo);CHKERRQ(ierr); 167 ierr = PetscMalloc2(nzmax,&ncols,nzmax,&nvals);CHKERRQ(ierr); 168 ierr = PetscArrayzero(nvals,nzmax);CHKERRQ(ierr); 169 for (i=0; i<m; i++) { 170 ierr = MatGetRow(B,i+rstart,&nzl,&cols,&vals);CHKERRQ(ierr); 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 ierr = MatSetValues(J,1,&nrow,2*nzl,ncols,nvals,INSERT_VALUES);CHKERRQ(ierr); 177 nrow = 2*(i+rstart) + 1; 178 ierr = MatSetValues(J,1,&nrow,2*nzl,ncols,nvals,INSERT_VALUES);CHKERRQ(ierr); 179 ierr = MatRestoreRow(B,i+rstart,&nzl,&cols,&vals);CHKERRQ(ierr); 180 } 181 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 182 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 183 if (viewMats) { 184 ierr = PetscPrintf(PETSC_COMM_WORLD,"Jacobian matrix structure:\n");CHKERRQ(ierr); 185 ierr = MatView(J,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 186 } 187 ierr = MatDestroy(&J);CHKERRQ(ierr); 188 ierr = PetscFree2(ncols,nvals);CHKERRQ(ierr); 189 } 190 191 /* 192 Free work space. All PETSc objects should be destroyed when they 193 are no longer needed. 194 */ 195 ierr = MatDestroy(&B);CHKERRQ(ierr); 196 ierr = VecDestroy(&xin);CHKERRQ(ierr); 197 ierr = PetscFinalize(); 198 return ierr; 199 } 200 201 /*TEST 202 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 207 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 214 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 221 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 228 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 235 236 TEST*/ 237