1 static char help[] = "Reads a PETSc matrix from a file partitions it\n\n"; 2 3 /* 4 Include "petscmat.h" so that we can use matrices. Note that this file 5 automatically includes: 6 petscsys.h - base PETSc routines petscvec.h - vectors 7 petscmat.h - matrices 8 petscis.h - index sets 9 petscviewer.h - viewers 10 11 Example of usage: 12 mpiexec -n 3 ex73 -f <matfile> -mat_partitioning_type parmetis/scotch -viewer_binary_skip_info -nox 13 */ 14 #include <petscmat.h> 15 16 int main(int argc, char **args) 17 { 18 MatType mtype = MATMPIAIJ; /* matrix format */ 19 Mat A, B; /* matrix */ 20 PetscViewer fd; /* viewer */ 21 char file[PETSC_MAX_PATH_LEN]; /* input file name */ 22 PetscBool flg, viewMats, viewIS, viewVecs, useND, noVecLoad = PETSC_FALSE; 23 PetscInt *nlocal, m, n; 24 PetscMPIInt rank, size; 25 MatPartitioning part; 26 IS is, isn; 27 Vec xin, xout; 28 VecScatter scat; 29 30 PetscFunctionBeginUser; 31 PetscCall(PetscInitialize(&argc, &args, NULL, help)); 32 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 33 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 34 PetscCall(PetscOptionsHasName(NULL, NULL, "-view_mats", &viewMats)); 35 PetscCall(PetscOptionsHasName(NULL, NULL, "-view_is", &viewIS)); 36 PetscCall(PetscOptionsHasName(NULL, NULL, "-view_vecs", &viewVecs)); 37 PetscCall(PetscOptionsHasName(NULL, NULL, "-use_nd", &useND)); 38 PetscCall(PetscOptionsHasName(NULL, NULL, "-novec_load", &noVecLoad)); 39 40 /* 41 Determine file from which we read the matrix 42 */ 43 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg)); 44 45 /* 46 Open binary file. Note that we use FILE_MODE_READ to indicate 47 reading from this file. 48 */ 49 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd)); 50 51 /* 52 Load the matrix and vector; then destroy the viewer. 53 */ 54 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 55 PetscCall(MatSetType(A, mtype)); 56 PetscCall(MatLoad(A, fd)); 57 if (!noVecLoad) { 58 PetscCall(VecCreate(PETSC_COMM_WORLD, &xin)); 59 PetscCall(VecLoad(xin, fd)); 60 } else { 61 PetscCall(MatCreateVecs(A, &xin, NULL)); 62 PetscCall(VecSetRandom(xin, NULL)); 63 } 64 PetscCall(PetscViewerDestroy(&fd)); 65 if (viewMats) { 66 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Original matrix:\n")); 67 PetscCall(MatView(A, PETSC_VIEWER_DRAW_WORLD)); 68 } 69 if (viewVecs) { 70 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Original vector:\n")); 71 PetscCall(VecView(xin, PETSC_VIEWER_STDOUT_WORLD)); 72 } 73 74 /* Partition the graph of the matrix */ 75 PetscCall(MatPartitioningCreate(PETSC_COMM_WORLD, &part)); 76 PetscCall(MatPartitioningSetAdjacency(part, A)); 77 PetscCall(MatPartitioningSetFromOptions(part)); 78 79 /* get new processor owner number of each vertex */ 80 if (useND) { 81 PetscCall(MatPartitioningApplyND(part, &is)); 82 } else { 83 PetscCall(MatPartitioningApply(part, &is)); 84 } 85 if (viewIS) { 86 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "IS1 - new processor ownership:\n")); 87 PetscCall(ISView(is, PETSC_VIEWER_STDOUT_WORLD)); 88 } 89 90 /* get new global number of each old global number */ 91 PetscCall(ISPartitioningToNumbering(is, &isn)); 92 if (viewIS) { 93 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "IS2 - new global numbering:\n")); 94 PetscCall(ISView(isn, PETSC_VIEWER_STDOUT_WORLD)); 95 } 96 97 /* get number of new vertices for each processor */ 98 PetscCall(PetscMalloc1(size, &nlocal)); 99 PetscCall(ISPartitioningCount(is, size, nlocal)); 100 PetscCall(ISDestroy(&is)); 101 102 /* get old global number of each new global number */ 103 PetscCall(ISInvertPermutation(isn, useND ? PETSC_DECIDE : nlocal[rank], &is)); 104 if (viewIS) { 105 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "IS3=inv(IS2) - old global number of each new global number:\n")); 106 PetscCall(ISView(is, PETSC_VIEWER_STDOUT_WORLD)); 107 } 108 109 /* move the matrix rows to the new processes they have been assigned to by the permutation */ 110 PetscCall(MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, &B)); 111 PetscCall(PetscFree(nlocal)); 112 PetscCall(ISDestroy(&isn)); 113 PetscCall(MatDestroy(&A)); 114 PetscCall(MatPartitioningDestroy(&part)); 115 if (viewMats) { 116 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Partitioned matrix:\n")); 117 PetscCall(MatView(B, PETSC_VIEWER_DRAW_WORLD)); 118 } 119 120 /* move the vector rows to the new processes they have been assigned to */ 121 PetscCall(MatGetLocalSize(B, &m, &n)); 122 PetscCall(VecCreateFromOptions(PETSC_COMM_WORLD, NULL, 1, m, PETSC_DECIDE, &xout)); 123 PetscCall(VecScatterCreate(xin, is, xout, NULL, &scat)); 124 PetscCall(VecScatterBegin(scat, xin, xout, INSERT_VALUES, SCATTER_FORWARD)); 125 PetscCall(VecScatterEnd(scat, xin, xout, INSERT_VALUES, SCATTER_FORWARD)); 126 PetscCall(VecScatterDestroy(&scat)); 127 if (viewVecs) { 128 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Mapped vector:\n")); 129 PetscCall(VecView(xout, PETSC_VIEWER_STDOUT_WORLD)); 130 } 131 PetscCall(VecDestroy(&xout)); 132 PetscCall(ISDestroy(&is)); 133 134 { 135 PetscInt rstart, i, *nzd, *nzo, nzl, nzmax = 0, *ncols, nrow, j; 136 Mat J; 137 const PetscInt *cols; 138 const PetscScalar *vals; 139 PetscScalar *nvals; 140 141 PetscCall(MatGetOwnershipRange(B, &rstart, NULL)); 142 PetscCall(PetscCalloc2(2 * m, &nzd, 2 * m, &nzo)); 143 for (i = 0; i < m; i++) { 144 PetscCall(MatGetRow(B, i + rstart, &nzl, &cols, NULL)); 145 for (j = 0; j < nzl; j++) { 146 if (cols[j] >= rstart && cols[j] < rstart + n) { 147 nzd[2 * i] += 2; 148 nzd[2 * i + 1] += 2; 149 } else { 150 nzo[2 * i] += 2; 151 nzo[2 * i + 1] += 2; 152 } 153 } 154 nzmax = PetscMax(nzmax, nzd[2 * i] + nzo[2 * i]); 155 PetscCall(MatRestoreRow(B, i + rstart, &nzl, &cols, NULL)); 156 } 157 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, 2 * m, 2 * m, PETSC_DECIDE, PETSC_DECIDE, 0, nzd, 0, nzo, &J)); 158 PetscCall(PetscInfo(0, "Created empty Jacobian matrix\n")); 159 PetscCall(PetscFree2(nzd, nzo)); 160 PetscCall(PetscMalloc2(nzmax, &ncols, nzmax, &nvals)); 161 PetscCall(PetscArrayzero(nvals, nzmax)); 162 for (i = 0; i < m; i++) { 163 PetscCall(MatGetRow(B, i + rstart, &nzl, &cols, &vals)); 164 for (j = 0; j < nzl; j++) { 165 ncols[2 * j] = 2 * cols[j]; 166 ncols[2 * j + 1] = 2 * cols[j] + 1; 167 } 168 nrow = 2 * (i + rstart); 169 PetscCall(MatSetValues(J, 1, &nrow, 2 * nzl, ncols, nvals, INSERT_VALUES)); 170 nrow = 2 * (i + rstart) + 1; 171 PetscCall(MatSetValues(J, 1, &nrow, 2 * nzl, ncols, nvals, INSERT_VALUES)); 172 PetscCall(MatRestoreRow(B, i + rstart, &nzl, &cols, &vals)); 173 } 174 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 175 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 176 if (viewMats) { 177 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Jacobian matrix nonzero structure:\n")); 178 PetscCall(MatView(J, PETSC_VIEWER_DRAW_WORLD)); 179 } 180 PetscCall(MatDestroy(&J)); 181 PetscCall(PetscFree2(ncols, nvals)); 182 } 183 184 /* 185 Free work space. All PETSc objects should be destroyed when they 186 are no longer needed. 187 */ 188 PetscCall(MatDestroy(&B)); 189 PetscCall(VecDestroy(&xin)); 190 PetscCall(PetscFinalize()); 191 return 0; 192 } 193 194 /*TEST 195 196 test: 197 nsize: 3 198 requires: parmetis datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 199 args: -nox -f ${DATAFILESPATH}/matrices/arco1 -mat_partitioning_type parmetis -viewer_binary_skip_info -novec_load 200 201 test: 202 requires: parmetis !complex double !defined(PETSC_USE_64BIT_INDICES) 203 output_file: output/ex73_1.out 204 suffix: parmetis_nd_32 205 nsize: 3 206 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 207 208 test: 209 requires: parmetis !complex double defined(PETSC_USE_64BIT_INDICES) 210 output_file: output/ex73_1.out 211 suffix: parmetis_nd_64 212 nsize: 3 213 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 214 215 test: 216 requires: ptscotch !complex double !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_SCOTCH_PARMETIS_V3_NODEND) 217 output_file: output/ex73_1.out 218 suffix: ptscotch_nd_32 219 nsize: 4 220 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 221 222 test: 223 requires: ptscotch !complex double defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_SCOTCH_PARMETIS_V3_NODEND) 224 output_file: output/ex73_1.out 225 suffix: ptscotch_nd_64 226 nsize: 4 227 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 228 229 TEST*/ 230