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