1 static char help[] = "Tests MatIncreaseOverlap() and MatCreateSubmatrices() for the parallel case.\n\ 2 This example is similar to ex40.c; here the index sets used are random.\n\ 3 Input arguments are:\n\ 4 -f <input_file> : file to load. For example see $PETSC_DIR/share/petsc/datafiles/matrices\n\ 5 -nd <size> : > 0 no of domains per processor \n\ 6 -ov <overlap> : >=0 amount of overlap between domains\n\n"; 7 8 #include <petscmat.h> 9 10 int main(int argc, char **args) 11 { 12 PetscInt nd = 2, ov = 1, i, j, lsize, m, n, *idx, bs; 13 PetscMPIInt rank, size; 14 PetscBool flg; 15 Mat A, B, *submatA, *submatB; 16 char file[PETSC_MAX_PATH_LEN]; 17 PetscViewer fd; 18 IS *is1, *is2; 19 PetscRandom r; 20 PetscBool test_unsorted = PETSC_FALSE; 21 PetscScalar rand; 22 23 PetscFunctionBeginUser; 24 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 25 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 26 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 27 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), NULL)); 28 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nd", &nd, NULL)); 29 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ov", &ov, NULL)); 30 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_unsorted", &test_unsorted, NULL)); 31 32 /* Read matrix A and RHS */ 33 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd)); 34 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 35 PetscCall(MatSetType(A, MATAIJ)); 36 PetscCall(MatSetFromOptions(A)); 37 PetscCall(MatLoad(A, fd)); 38 PetscCall(PetscViewerDestroy(&fd)); 39 40 /* Read the same matrix as a seq matrix B */ 41 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, file, FILE_MODE_READ, &fd)); 42 PetscCall(MatCreate(PETSC_COMM_SELF, &B)); 43 PetscCall(MatSetType(B, MATSEQAIJ)); 44 PetscCall(MatSetFromOptions(B)); 45 PetscCall(MatLoad(B, fd)); 46 PetscCall(PetscViewerDestroy(&fd)); 47 48 PetscCall(MatGetBlockSize(A, &bs)); 49 50 /* Create the Random no generator */ 51 PetscCall(MatGetSize(A, &m, &n)); 52 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &r)); 53 PetscCall(PetscRandomSetFromOptions(r)); 54 55 /* Create the IS corresponding to subdomains */ 56 PetscCall(PetscMalloc1(nd, &is1)); 57 PetscCall(PetscMalloc1(nd, &is2)); 58 PetscCall(PetscMalloc1(m, &idx)); 59 for (i = 0; i < m; i++) idx[i] = i; 60 61 /* Create the random Index Sets */ 62 for (i = 0; i < nd; i++) { 63 /* Skip a few,so that the IS on different procs are different*/ 64 for (j = 0; j < rank; j++) PetscCall(PetscRandomGetValue(r, &rand)); 65 PetscCall(PetscRandomGetValue(r, &rand)); 66 lsize = (PetscInt)(rand * (m / bs)); 67 /* shuffle */ 68 for (j = 0; j < lsize; j++) { 69 PetscInt k, swap, l; 70 71 PetscCall(PetscRandomGetValue(r, &rand)); 72 k = j + (PetscInt)(rand * ((m / bs) - j)); 73 for (l = 0; l < bs; l++) { 74 swap = idx[bs * j + l]; 75 idx[bs * j + l] = idx[bs * k + l]; 76 idx[bs * k + l] = swap; 77 } 78 } 79 if (!test_unsorted) PetscCall(PetscSortInt(lsize * bs, idx)); 80 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, lsize * bs, idx, PETSC_COPY_VALUES, is1 + i)); 81 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, lsize * bs, idx, PETSC_COPY_VALUES, is2 + i)); 82 PetscCall(ISSetBlockSize(is1[i], bs)); 83 PetscCall(ISSetBlockSize(is2[i], bs)); 84 } 85 86 if (!test_unsorted) { 87 PetscCall(MatIncreaseOverlap(A, nd, is1, ov)); 88 PetscCall(MatIncreaseOverlap(B, nd, is2, ov)); 89 90 for (i = 0; i < nd; ++i) { 91 PetscCall(ISSort(is1[i])); 92 PetscCall(ISSort(is2[i])); 93 } 94 } 95 96 PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_INITIAL_MATRIX, &submatA)); 97 PetscCall(MatCreateSubMatrices(B, nd, is2, is2, MAT_INITIAL_MATRIX, &submatB)); 98 99 /* Now see if the serial and parallel case have the same answers */ 100 for (i = 0; i < nd; ++i) { 101 PetscCall(MatEqual(submatA[i], submatB[i], &flg)); 102 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%" PetscInt_FMT "-th parallel submatA != seq submatB", i); 103 } 104 105 /* Free Allocated Memory */ 106 for (i = 0; i < nd; ++i) { 107 PetscCall(ISDestroy(&is1[i])); 108 PetscCall(ISDestroy(&is2[i])); 109 } 110 PetscCall(MatDestroySubMatrices(nd, &submatA)); 111 PetscCall(MatDestroySubMatrices(nd, &submatB)); 112 113 PetscCall(PetscRandomDestroy(&r)); 114 PetscCall(PetscFree(is1)); 115 PetscCall(PetscFree(is2)); 116 PetscCall(MatDestroy(&A)); 117 PetscCall(MatDestroy(&B)); 118 PetscCall(PetscFree(idx)); 119 PetscCall(PetscFinalize()); 120 return 0; 121 } 122 123 /*TEST 124 125 build: 126 requires: !complex 127 128 test: 129 nsize: 3 130 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) !complex 131 args: -f ${DATAFILESPATH}/matrices/arco1 -nd 5 -ov 2 132 133 test: 134 suffix: 2 135 args: -f ${DATAFILESPATH}/matrices/arco1 -nd 8 -ov 2 136 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) !complex 137 138 test: 139 suffix: unsorted_baij_mpi 140 nsize: 3 141 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) !complex 142 args: -f ${DATAFILESPATH}/matrices/cfd.1.10 -nd 8 -mat_type baij -test_unsorted 143 144 test: 145 suffix: unsorted_baij_seq 146 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) !complex 147 args: -f ${DATAFILESPATH}/matrices/cfd.1.10 -nd 8 -mat_type baij -test_unsorted 148 149 test: 150 suffix: unsorted_mpi 151 nsize: 3 152 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) !complex 153 args: -f ${DATAFILESPATH}/matrices/arco1 -nd 8 -test_unsorted 154 155 test: 156 suffix: unsorted_seq 157 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) !complex 158 args: -f ${DATAFILESPATH}/matrices/arco1 -nd 8 -test_unsorted 159 160 TEST*/ 161