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 PetscErrorCode ierr; 35 36 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 37 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 38 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 39 PetscCheckFalse(size>2,PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG, "A uniprocessor or two-processor example only."); 40 41 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A)); 42 if (size > 1) { 43 n = 8; N = 16; 44 } else { 45 n = 16; N = 16; 46 } 47 CHKERRQ(MatSetSizes(A,n,n,N,N)); 48 CHKERRQ(MatSetFromOptions(A)); 49 CHKERRQ(MatSetUp(A)); 50 51 /* Don't care if the entries are set multiple times by different procs. */ 52 for (i=0; i<4; ++i) { 53 for (j = 0; j<4; ++j) { 54 row = j*4+i; 55 v = -1.0; 56 if (i>0) { 57 col = row-1; CHKERRQ(MatSetValues(A,1,&row,1,&col,&v,INSERT_VALUES)); 58 } 59 if (i<3) { 60 col = row+1; CHKERRQ(MatSetValues(A,1,&row,1,&col,&v,INSERT_VALUES)); 61 } 62 if (j>0) { 63 col = row-4; CHKERRQ(MatSetValues(A,1,&row,1,&col,&v,INSERT_VALUES)); 64 } 65 if (j<3) { 66 col = row+4; CHKERRQ(MatSetValues(A,1,&row,1,&col,&v,INSERT_VALUES)); 67 } 68 v = 4.0; 69 CHKERRQ(MatSetValues(A,1,&row,1,&row,&v,INSERT_VALUES)); 70 } 71 } 72 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 73 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 74 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "Original matrix\n")); 75 CHKERRQ(MatView(A,PETSC_VIEWER_STDOUT_WORLD)); 76 77 if (size > 1) { 78 nsub = 1; /* one subdomain per rank */ 79 } 80 else { 81 nsub = 2; /* both subdomains on rank 0 */ 82 } 83 if (rank) { 84 jlow = Jlow+1; jhigh = Jhigh+1; 85 } 86 else { 87 jlow = Jlow; jhigh = Jhigh; 88 } 89 sort_rows = PETSC_FALSE; 90 CHKERRQ(PetscOptionsGetBool(NULL,NULL, "-sort_rows", &sort_rows, NULL)); 91 sort_cols = PETSC_FALSE; 92 CHKERRQ(PetscOptionsGetBool(NULL,NULL, "-sort_cols", &sort_cols, NULL)); 93 for (l = 0; l < nsub; ++l) { 94 CHKERRQ(PetscMalloc1(12, &subindices)); 95 k = 0; 96 for (i = 0; i < 4; ++i) { 97 for (j = jlow[l]; j < jhigh[l]; ++j) { 98 subindices[k] = j*4+i; 99 k++; 100 } 101 } 102 CHKERRQ(ISCreateGeneral(PETSC_COMM_SELF, 12, subindices, PETSC_OWN_POINTER, rowis+l)); 103 if ((sort_rows && !sort_cols) || (!sort_rows && sort_cols)) { 104 CHKERRQ(ISDuplicate(rowis[l],colis+l)); 105 } else { 106 CHKERRQ(PetscObjectReference((PetscObject)rowis[l])); 107 colis[l] = rowis[l]; 108 } 109 if (sort_rows) { 110 CHKERRQ(ISSort(rowis[l])); 111 } 112 if (sort_cols) { 113 CHKERRQ(ISSort(colis[l])); 114 } 115 } 116 117 CHKERRQ(MatCreateSubMatrices(A,nsub,rowis,colis,MAT_INITIAL_MATRIX, &S)); 118 119 show_inversions = PETSC_FALSE; 120 121 CHKERRQ(PetscOptionsGetBool(NULL,NULL, "-show_inversions", &show_inversions, NULL)); 122 123 inversions = 0; 124 for (p = 0; p < size; ++p) { 125 if (p == rank) { 126 CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "[%" PetscInt_FMT ":%" PetscInt_FMT "]: Number of subdomains: %" PetscInt_FMT ":\n", rank, size, nsub)); 127 for (l = 0; l < nsub; ++l) { 128 PetscInt i0, i1; 129 CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "[%" PetscInt_FMT ":%" PetscInt_FMT "]: Subdomain row IS %" PetscInt_FMT ":\n", rank, size, l)); 130 CHKERRQ(ISView(rowis[l],PETSC_VIEWER_STDOUT_SELF)); 131 CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "[%" PetscInt_FMT ":%" PetscInt_FMT "]: Subdomain col IS %" PetscInt_FMT ":\n", rank, size, l)); 132 CHKERRQ(ISView(colis[l],PETSC_VIEWER_STDOUT_SELF)); 133 CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "[%" PetscInt_FMT ":%" PetscInt_FMT "]: Submatrix %" PetscInt_FMT ":\n", rank, size, l)); 134 CHKERRQ(MatView(S[l],PETSC_VIEWER_STDOUT_SELF)); 135 if (show_inversions) { 136 CHKERRQ(MatGetOwnershipRange(S[l], &i0,&i1)); 137 for (i = i0; i < i1; ++i) { 138 CHKERRQ(MatGetRow(S[l], i, &ncols, &cols, NULL)); 139 for (j = 1; j < ncols; ++j) { 140 if (cols[j] < cols[j-1]) { 141 CHKERRQ(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)); 142 inversions++; 143 } 144 } 145 CHKERRQ(MatRestoreRow(S[l], i, &ncols, &cols, NULL)); 146 } 147 } 148 } 149 } 150 CHKERRMPI(MPI_Barrier(PETSC_COMM_WORLD)); 151 } 152 if (show_inversions) { 153 CHKERRMPI(MPI_Reduce(&inversions,&total_inversions,1,MPIU_INT, MPI_SUM,0,PETSC_COMM_WORLD)); 154 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "*Total inversions: %" PetscInt_FMT "\n", total_inversions)); 155 } 156 CHKERRQ(MatDestroy(&A)); 157 158 for (l = 0; l < nsub; ++l) { 159 CHKERRQ(ISDestroy(&(rowis[l]))); 160 CHKERRQ(ISDestroy(&(colis[l]))); 161 } 162 CHKERRQ(MatDestroySubMatrices(nsub,&S)); 163 ierr = PetscFinalize(); 164 return ierr; 165 } 166