1 2 static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for parallel MatSBAIJ format.\n"; 3 /* Example of usage: 4 mpiexec -n 2 ./ex92 -nd 2 -ov 3 -mat_block_size 2 -view_id 0 -test_overlap -test_submat 5 */ 6 #include <petscmat.h> 7 8 int main(int argc, char **args) 9 { 10 Mat A, Atrans, sA, *submatA, *submatsA; 11 PetscMPIInt size, rank; 12 PetscInt bs = 1, mbs = 10, ov = 1, i, j, k, *rows, *cols, nd = 2, *idx, rstart, rend, sz, M, N, Mbs; 13 PetscScalar *vals, rval, one = 1.0; 14 IS *is1, *is2; 15 PetscRandom rand; 16 PetscBool flg, TestOverlap, TestSubMat, TestAllcols, test_sorted = PETSC_FALSE; 17 PetscInt vid = -1; 18 PetscLogStage stages[2]; 19 20 PetscFunctionBeginUser; 21 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 22 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 23 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 24 25 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL)); 26 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_mbs", &mbs, NULL)); 27 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ov", &ov, NULL)); 28 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nd", &nd, NULL)); 29 PetscCall(PetscOptionsGetInt(NULL, NULL, "-view_id", &vid, NULL)); 30 PetscCall(PetscOptionsHasName(NULL, NULL, "-test_overlap", &TestOverlap)); 31 PetscCall(PetscOptionsHasName(NULL, NULL, "-test_submat", &TestSubMat)); 32 PetscCall(PetscOptionsHasName(NULL, NULL, "-test_allcols", &TestAllcols)); 33 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_sorted", &test_sorted, NULL)); 34 35 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 36 PetscCall(MatSetSizes(A, mbs * bs, mbs * bs, PETSC_DECIDE, PETSC_DECIDE)); 37 PetscCall(MatSetType(A, MATBAIJ)); 38 PetscCall(MatSeqBAIJSetPreallocation(A, bs, PETSC_DEFAULT, NULL)); 39 PetscCall(MatMPIBAIJSetPreallocation(A, bs, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL)); 40 41 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand)); 42 PetscCall(PetscRandomSetFromOptions(rand)); 43 44 PetscCall(MatGetOwnershipRange(A, &rstart, &rend)); 45 PetscCall(MatGetSize(A, &M, &N)); 46 Mbs = M / bs; 47 48 PetscCall(PetscMalloc1(bs, &rows)); 49 PetscCall(PetscMalloc1(bs, &cols)); 50 PetscCall(PetscMalloc1(bs * bs, &vals)); 51 PetscCall(PetscMalloc1(M, &idx)); 52 53 /* Now set blocks of values */ 54 for (j = 0; j < bs * bs; j++) vals[j] = 0.0; 55 for (i = 0; i < Mbs; i++) { 56 cols[0] = i * bs; 57 rows[0] = i * bs; 58 for (j = 1; j < bs; j++) { 59 rows[j] = rows[j - 1] + 1; 60 cols[j] = cols[j - 1] + 1; 61 } 62 PetscCall(MatSetValues(A, bs, rows, bs, cols, vals, ADD_VALUES)); 63 } 64 /* second, add random blocks */ 65 for (i = 0; i < 20 * bs; i++) { 66 PetscCall(PetscRandomGetValue(rand, &rval)); 67 cols[0] = bs * (PetscInt)(PetscRealPart(rval) * Mbs); 68 PetscCall(PetscRandomGetValue(rand, &rval)); 69 rows[0] = rstart + bs * (PetscInt)(PetscRealPart(rval) * mbs); 70 for (j = 1; j < bs; j++) { 71 rows[j] = rows[j - 1] + 1; 72 cols[j] = cols[j - 1] + 1; 73 } 74 75 for (j = 0; j < bs * bs; j++) { 76 PetscCall(PetscRandomGetValue(rand, &rval)); 77 vals[j] = rval; 78 } 79 PetscCall(MatSetValues(A, bs, rows, bs, cols, vals, ADD_VALUES)); 80 } 81 82 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 83 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 84 85 /* make A a symmetric matrix: A <- A^T + A */ 86 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Atrans)); 87 PetscCall(MatAXPY(A, one, Atrans, DIFFERENT_NONZERO_PATTERN)); 88 PetscCall(MatDestroy(&Atrans)); 89 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Atrans)); 90 PetscCall(MatEqual(A, Atrans, &flg)); 91 if (flg) { 92 PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE)); 93 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "A+A^T is non-symmetric"); 94 PetscCall(MatDestroy(&Atrans)); 95 96 /* create a SeqSBAIJ matrix sA (= A) */ 97 PetscCall(MatConvert(A, MATSBAIJ, MAT_INITIAL_MATRIX, &sA)); 98 if (vid >= 0 && vid < size) { 99 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "A:\n")); 100 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 101 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "sA:\n")); 102 PetscCall(MatView(sA, PETSC_VIEWER_STDOUT_WORLD)); 103 } 104 105 /* Test sA==A through MatMult() */ 106 PetscCall(MatMultEqual(A, sA, 10, &flg)); 107 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Error in MatConvert(): A != sA"); 108 109 /* Test MatIncreaseOverlap() */ 110 PetscCall(PetscMalloc1(nd, &is1)); 111 PetscCall(PetscMalloc1(nd, &is2)); 112 113 for (i = 0; i < nd; i++) { 114 if (!TestAllcols) { 115 PetscCall(PetscRandomGetValue(rand, &rval)); 116 sz = (PetscInt)((0.5 + 0.2 * PetscRealPart(rval)) * mbs); /* 0.5*mbs < sz < 0.7*mbs */ 117 118 for (j = 0; j < sz; j++) { 119 PetscCall(PetscRandomGetValue(rand, &rval)); 120 idx[j * bs] = bs * (PetscInt)(PetscRealPart(rval) * Mbs); 121 for (k = 1; k < bs; k++) idx[j * bs + k] = idx[j * bs] + k; 122 } 123 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, sz * bs, idx, PETSC_COPY_VALUES, is1 + i)); 124 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, sz * bs, idx, PETSC_COPY_VALUES, is2 + i)); 125 if (rank == vid) { 126 PetscCall(PetscPrintf(PETSC_COMM_SELF, " [%d] IS sz[%" PetscInt_FMT "]: %" PetscInt_FMT "\n", rank, i, sz)); 127 PetscCall(ISView(is2[i], PETSC_VIEWER_STDOUT_SELF)); 128 } 129 } else { /* Test all rows and columns */ 130 sz = M; 131 PetscCall(ISCreateStride(PETSC_COMM_SELF, sz, 0, 1, is1 + i)); 132 PetscCall(ISCreateStride(PETSC_COMM_SELF, sz, 0, 1, is2 + i)); 133 134 if (rank == vid) { 135 PetscBool colflag; 136 PetscCall(ISIdentity(is2[i], &colflag)); 137 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] is2[%" PetscInt_FMT "], colflag %d\n", rank, i, colflag)); 138 PetscCall(ISView(is2[i], PETSC_VIEWER_STDOUT_SELF)); 139 } 140 } 141 } 142 143 PetscCall(PetscLogStageRegister("MatOv_SBAIJ", &stages[0])); 144 PetscCall(PetscLogStageRegister("MatOv_BAIJ", &stages[1])); 145 146 /* Test MatIncreaseOverlap */ 147 if (TestOverlap) { 148 PetscCall(PetscLogStagePush(stages[0])); 149 PetscCall(MatIncreaseOverlap(sA, nd, is2, ov)); 150 PetscCall(PetscLogStagePop()); 151 152 PetscCall(PetscLogStagePush(stages[1])); 153 PetscCall(MatIncreaseOverlap(A, nd, is1, ov)); 154 PetscCall(PetscLogStagePop()); 155 156 if (rank == vid) { 157 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n[%d] IS from BAIJ:\n", rank)); 158 PetscCall(ISView(is1[0], PETSC_VIEWER_STDOUT_SELF)); 159 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n[%d] IS from SBAIJ:\n", rank)); 160 PetscCall(ISView(is2[0], PETSC_VIEWER_STDOUT_SELF)); 161 } 162 163 for (i = 0; i < nd; ++i) { 164 PetscCall(ISEqual(is1[i], is2[i], &flg)); 165 if (!flg) { 166 if (rank == 0) { 167 PetscCall(ISSort(is1[i])); 168 PetscCall(ISSort(is2[i])); 169 } 170 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "i=%" PetscInt_FMT ", is1 != is2", i); 171 } 172 } 173 } 174 175 /* Test MatCreateSubmatrices */ 176 if (TestSubMat) { 177 if (test_sorted) { 178 for (i = 0; i < nd; ++i) PetscCall(ISSort(is1[i])); 179 } 180 PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_INITIAL_MATRIX, &submatA)); 181 PetscCall(MatCreateSubMatrices(sA, nd, is1, is1, MAT_INITIAL_MATRIX, &submatsA)); 182 183 PetscCall(MatMultEqual(A, sA, 10, &flg)); 184 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "A != sA"); 185 186 /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */ 187 PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_REUSE_MATRIX, &submatA)); 188 PetscCall(MatCreateSubMatrices(sA, nd, is1, is1, MAT_REUSE_MATRIX, &submatsA)); 189 PetscCall(MatMultEqual(A, sA, 10, &flg)); 190 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "MatCreateSubmatrices(): A != sA"); 191 192 PetscCall(MatDestroySubMatrices(nd, &submatA)); 193 PetscCall(MatDestroySubMatrices(nd, &submatsA)); 194 } 195 196 /* Free allocated memory */ 197 for (i = 0; i < nd; ++i) { 198 PetscCall(ISDestroy(&is1[i])); 199 PetscCall(ISDestroy(&is2[i])); 200 } 201 PetscCall(PetscFree(is1)); 202 PetscCall(PetscFree(is2)); 203 PetscCall(PetscFree(idx)); 204 PetscCall(PetscFree(rows)); 205 PetscCall(PetscFree(cols)); 206 PetscCall(PetscFree(vals)); 207 PetscCall(MatDestroy(&A)); 208 PetscCall(MatDestroy(&sA)); 209 PetscCall(PetscRandomDestroy(&rand)); 210 PetscCall(PetscFinalize()); 211 return 0; 212 } 213 214 /*TEST 215 216 test: 217 args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_submat 218 output_file: output/ex92_1.out 219 220 test: 221 suffix: 2 222 nsize: {{3 4}} 223 args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_submat 224 output_file: output/ex92_1.out 225 226 test: 227 suffix: 3 228 nsize: {{3 4}} 229 args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_allcols 230 output_file: output/ex92_1.out 231 232 test: 233 suffix: 3_sorted 234 nsize: {{3 4}} 235 args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_allcols -test_sorted 236 output_file: output/ex92_1.out 237 238 test: 239 suffix: 4 240 nsize: {{3 4}} 241 args: -ov {{1 3}} -mat_block_size {{2 8}} -test_submat -test_allcols 242 output_file: output/ex92_1.out 243 244 TEST*/ 245