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 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr); 38 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 39 if (size>2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG, "A uniprocessor or two-processor example only.\n"); 40 41 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 42 if (size > 1) { 43 n = 8; N = 16; 44 } else { 45 n = 16; N = 16; 46 } 47 ierr = MatSetSizes(A,n,n,N,N);CHKERRQ(ierr); 48 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 49 ierr = MatSetUp(A);CHKERRQ(ierr); 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; ierr = MatSetValues(A,1,&row,1,&col,&v,INSERT_VALUES);CHKERRQ(ierr); 58 } 59 if (i<3) { 60 col = row+1; ierr = MatSetValues(A,1,&row,1,&col,&v,INSERT_VALUES);CHKERRQ(ierr); 61 } 62 if (j>0) { 63 col = row-4; ierr = MatSetValues(A,1,&row,1,&col,&v,INSERT_VALUES);CHKERRQ(ierr); 64 } 65 if (j<3) { 66 col = row+4; ierr = MatSetValues(A,1,&row,1,&col,&v,INSERT_VALUES);CHKERRQ(ierr); 67 } 68 v = 4.0; 69 ierr = MatSetValues(A,1,&row,1,&row,&v,INSERT_VALUES);CHKERRQ(ierr); 70 } 71 } 72 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 73 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 74 ierr = PetscPrintf(PETSC_COMM_WORLD, "Original matrix\n");CHKERRQ(ierr); 75 ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 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 ierr = PetscOptionsGetBool(NULL,NULL, "-sort_rows", &sort_rows, NULL);CHKERRQ(ierr); 91 sort_cols = PETSC_FALSE; 92 ierr = PetscOptionsGetBool(NULL,NULL, "-sort_cols", &sort_cols, NULL);CHKERRQ(ierr); 93 for (l = 0; l < nsub; ++l) { 94 ierr = PetscMalloc1(12, &subindices);CHKERRQ(ierr); 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 ierr = ISCreateGeneral(PETSC_COMM_SELF, 12, subindices, PETSC_OWN_POINTER, rowis+l);CHKERRQ(ierr); 103 if ((sort_rows && !sort_cols) || (!sort_rows && sort_cols)) { 104 ierr = ISDuplicate(rowis[l],colis+l);CHKERRQ(ierr); 105 } else { 106 ierr = PetscObjectReference((PetscObject)rowis[l]);CHKERRQ(ierr); 107 colis[l] = rowis[l]; 108 } 109 if (sort_rows) { 110 ierr = ISSort(rowis[l]);CHKERRQ(ierr); 111 } 112 if (sort_cols) { 113 ierr = ISSort(colis[l]);CHKERRQ(ierr); 114 } 115 } 116 117 ierr = MatCreateSubMatrices(A,nsub,rowis,colis,MAT_INITIAL_MATRIX, &S);CHKERRQ(ierr); 118 119 show_inversions = PETSC_FALSE; 120 121 ierr = PetscOptionsGetBool(NULL,NULL, "-show_inversions", &show_inversions, NULL);CHKERRQ(ierr); 122 123 inversions = 0; 124 for (p = 0; p < size; ++p) { 125 if (p == rank) { 126 ierr = PetscPrintf(PETSC_COMM_SELF, "[%D:%D]: Number of subdomains: %D:\n", rank, size, nsub);CHKERRQ(ierr); 127 for (l = 0; l < nsub; ++l) { 128 PetscInt i0, i1; 129 ierr = PetscPrintf(PETSC_COMM_SELF, "[%D:%D]: Subdomain row IS %D:\n", rank, size, l);CHKERRQ(ierr); 130 ierr = ISView(rowis[l],PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 131 ierr = PetscPrintf(PETSC_COMM_SELF, "[%D:%D]: Subdomain col IS %D:\n", rank, size, l);CHKERRQ(ierr); 132 ierr = ISView(colis[l],PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 133 ierr = PetscPrintf(PETSC_COMM_SELF, "[%D:%D]: Submatrix %D:\n", rank, size, l);CHKERRQ(ierr); 134 ierr = MatView(S[l],PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 135 if (show_inversions) { 136 ierr = MatGetOwnershipRange(S[l], &i0,&i1);CHKERRQ(ierr); 137 for (i = i0; i < i1; ++i) { 138 ierr = MatGetRow(S[l], i, &ncols, &cols, NULL);CHKERRQ(ierr); 139 for (j = 1; j < ncols; ++j) { 140 if (cols[j] < cols[j-1]) { 141 ierr = PetscPrintf(PETSC_COMM_SELF, "***Inversion in row %D: col[%D] = %D < %D = col[%D]\n", i, j, cols[j], cols[j-1], j-1);CHKERRQ(ierr); 142 inversions++; 143 } 144 } 145 ierr = MatRestoreRow(S[l], i, &ncols, &cols, NULL);CHKERRQ(ierr); 146 } 147 } 148 } 149 } 150 ierr = MPI_Barrier(PETSC_COMM_WORLD);CHKERRMPI(ierr); 151 } 152 if (show_inversions) { 153 ierr = MPI_Reduce(&inversions,&total_inversions,1,MPIU_INT, MPI_SUM,0,PETSC_COMM_WORLD);CHKERRMPI(ierr); 154 ierr = PetscPrintf(PETSC_COMM_WORLD, "*Total inversions: %D\n", total_inversions);CHKERRQ(ierr); 155 } 156 ierr = MatDestroy(&A);CHKERRQ(ierr); 157 158 for (l = 0; l < nsub; ++l) { 159 ierr = ISDestroy(&(rowis[l]));CHKERRQ(ierr); 160 ierr = ISDestroy(&(colis[l]));CHKERRQ(ierr); 161 } 162 ierr = MatDestroySubMatrices(nsub,&S);CHKERRQ(ierr); 163 ierr = PetscFinalize(); 164 return ierr; 165 } 166