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