1 2 static char help[] = "Extract submatrices using unsorted indices. For SEQSBAIJ either sort both rows and columns, or sort none.\n\n"; 3 /* 4 Take a 4x4 grid and form a 5-point stencil graph Laplacian over it. 5 Partition the grid into two subdomains by splitting into two in the j-direction (slowest varying). 6 Impose an overlap of 1 and order the subdomains with the j-direction varying fastest. 7 Extract the subdomain submatrices, one per rank. 8 */ 9 /* Results: 10 Sequential: 11 - seqaij: will error out, if rows or columns are unsorted 12 - seqbaij: will extract submatrices correctly even for unsorted row or column indices 13 - seqsbaij: will extract submatrices correctly even for unsorted row and column indices (both must be sorted or not); 14 CANNOT automatically report inversions, because MatGetRow is not available. 15 MPI: 16 - mpiaij: will error out, if columns are unsorted 17 - mpibaij: will error out, if columns are unsorted. 18 - mpisbaij: will error out, if columns are unsorted; even with unsorted rows will produce correct submatrices; 19 CANNOT automatically report inversions, because MatGetRow is not available. 20 */ 21 22 #include <petscmat.h> 23 #include <petscis.h> 24 25 int main(int argc,char **args) 26 { 27 Mat A, *S; 28 IS rowis[2], colis[2]; 29 PetscInt n,N,i,j,k,l,nsub,Jlow[2] = {0,1}, *jlow, Jhigh[2] = {3,4}, *jhigh, row, col, *subindices, ncols; 30 const PetscInt *cols; 31 PetscScalar v; 32 PetscMPIInt rank, size, p, inversions, total_inversions; 33 PetscBool sort_rows, sort_cols, show_inversions; 34 35 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 36 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 37 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 38 PetscCheck(size < 3,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE, "A uniprocessor or two-processor example only."); 39 40 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 41 if (size > 1) { 42 n = 8; N = 16; 43 } else { 44 n = 16; N = 16; 45 } 46 PetscCall(MatSetSizes(A,n,n,N,N)); 47 PetscCall(MatSetFromOptions(A)); 48 PetscCall(MatSetUp(A)); 49 50 /* Don't care if the entries are set multiple times by different procs. */ 51 for (i=0; i<4; ++i) { 52 for (j = 0; j<4; ++j) { 53 row = j*4+i; 54 v = -1.0; 55 if (i>0) { 56 col = row-1; PetscCall(MatSetValues(A,1,&row,1,&col,&v,INSERT_VALUES)); 57 } 58 if (i<3) { 59 col = row+1; PetscCall(MatSetValues(A,1,&row,1,&col,&v,INSERT_VALUES)); 60 } 61 if (j>0) { 62 col = row-4; PetscCall(MatSetValues(A,1,&row,1,&col,&v,INSERT_VALUES)); 63 } 64 if (j<3) { 65 col = row+4; PetscCall(MatSetValues(A,1,&row,1,&col,&v,INSERT_VALUES)); 66 } 67 v = 4.0; 68 PetscCall(MatSetValues(A,1,&row,1,&row,&v,INSERT_VALUES)); 69 } 70 } 71 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 72 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 73 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Original matrix\n")); 74 PetscCall(MatView(A,PETSC_VIEWER_STDOUT_WORLD)); 75 76 if (size > 1) { 77 nsub = 1; /* one subdomain per rank */ 78 } 79 else { 80 nsub = 2; /* both subdomains on rank 0 */ 81 } 82 if (rank) { 83 jlow = Jlow+1; jhigh = Jhigh+1; 84 } 85 else { 86 jlow = Jlow; jhigh = Jhigh; 87 } 88 sort_rows = PETSC_FALSE; 89 PetscCall(PetscOptionsGetBool(NULL,NULL, "-sort_rows", &sort_rows, NULL)); 90 sort_cols = PETSC_FALSE; 91 PetscCall(PetscOptionsGetBool(NULL,NULL, "-sort_cols", &sort_cols, NULL)); 92 for (l = 0; l < nsub; ++l) { 93 PetscCall(PetscMalloc1(12, &subindices)); 94 k = 0; 95 for (i = 0; i < 4; ++i) { 96 for (j = jlow[l]; j < jhigh[l]; ++j) { 97 subindices[k] = j*4+i; 98 k++; 99 } 100 } 101 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 12, subindices, PETSC_OWN_POINTER, rowis+l)); 102 if ((sort_rows && !sort_cols) || (!sort_rows && sort_cols)) { 103 PetscCall(ISDuplicate(rowis[l],colis+l)); 104 } else { 105 PetscCall(PetscObjectReference((PetscObject)rowis[l])); 106 colis[l] = rowis[l]; 107 } 108 if (sort_rows) { 109 PetscCall(ISSort(rowis[l])); 110 } 111 if (sort_cols) { 112 PetscCall(ISSort(colis[l])); 113 } 114 } 115 116 PetscCall(MatCreateSubMatrices(A,nsub,rowis,colis,MAT_INITIAL_MATRIX, &S)); 117 118 show_inversions = PETSC_FALSE; 119 120 PetscCall(PetscOptionsGetBool(NULL,NULL, "-show_inversions", &show_inversions, NULL)); 121 122 inversions = 0; 123 for (p = 0; p < size; ++p) { 124 if (p == rank) { 125 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%" PetscInt_FMT ":%" PetscInt_FMT "]: Number of subdomains: %" PetscInt_FMT ":\n", rank, size, nsub)); 126 for (l = 0; l < nsub; ++l) { 127 PetscInt i0, i1; 128 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%" PetscInt_FMT ":%" PetscInt_FMT "]: Subdomain row IS %" PetscInt_FMT ":\n", rank, size, l)); 129 PetscCall(ISView(rowis[l],PETSC_VIEWER_STDOUT_SELF)); 130 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%" PetscInt_FMT ":%" PetscInt_FMT "]: Subdomain col IS %" PetscInt_FMT ":\n", rank, size, l)); 131 PetscCall(ISView(colis[l],PETSC_VIEWER_STDOUT_SELF)); 132 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%" PetscInt_FMT ":%" PetscInt_FMT "]: Submatrix %" PetscInt_FMT ":\n", rank, size, l)); 133 PetscCall(MatView(S[l],PETSC_VIEWER_STDOUT_SELF)); 134 if (show_inversions) { 135 PetscCall(MatGetOwnershipRange(S[l], &i0,&i1)); 136 for (i = i0; i < i1; ++i) { 137 PetscCall(MatGetRow(S[l], i, &ncols, &cols, NULL)); 138 for (j = 1; j < ncols; ++j) { 139 if (cols[j] < cols[j-1]) { 140 PetscCall(PetscPrintf(PETSC_COMM_SELF, "***Inversion in row %" PetscInt_FMT ": col[%" PetscInt_FMT "] = %" PetscInt_FMT " < %" PetscInt_FMT " = col[%" PetscInt_FMT "]\n", i, j, cols[j], cols[j-1], j-1)); 141 inversions++; 142 } 143 } 144 PetscCall(MatRestoreRow(S[l], i, &ncols, &cols, NULL)); 145 } 146 } 147 } 148 } 149 PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD)); 150 } 151 if (show_inversions) { 152 PetscCallMPI(MPI_Reduce(&inversions,&total_inversions,1,MPIU_INT, MPI_SUM,0,PETSC_COMM_WORLD)); 153 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "*Total inversions: %" PetscInt_FMT "\n", total_inversions)); 154 } 155 PetscCall(MatDestroy(&A)); 156 157 for (l = 0; l < nsub; ++l) { 158 PetscCall(ISDestroy(&(rowis[l]))); 159 PetscCall(ISDestroy(&(colis[l]))); 160 } 161 PetscCall(MatDestroySubMatrices(nsub,&S)); 162 PetscCall(PetscFinalize()); 163 return 0; 164 } 165