1 2 static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for parallel AIJ and BAIJ formats.\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc, char **args) 7 { 8 Mat E, A, B, *submatA, *submatB; 9 PetscInt bs = 1, m = 11, ov = 1, i, j, k, *rows, *cols, nd = 5, *idx, rstart, rend, sz, mm, nn, M, N, Mbs; 10 PetscMPIInt size, rank; 11 PetscScalar *vals, rval; 12 IS *is1, *is2; 13 PetscRandom rdm; 14 Vec xx, s1, s2; 15 PetscReal s1norm, s2norm, rnorm, tol = 100 * PETSC_SMALL; 16 PetscBool flg, test_nd0 = PETSC_FALSE, emptynd; 17 18 PetscFunctionBeginUser; 19 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 20 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 21 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 22 23 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL)); 24 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_size", &m, NULL)); 25 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ov", &ov, NULL)); 26 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nd", &nd, NULL)); 27 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_nd0", &test_nd0, NULL)); 28 29 /* Create a AIJ matrix A */ 30 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 31 PetscCall(MatSetSizes(A, m * bs, m * bs, PETSC_DECIDE, PETSC_DECIDE)); 32 PetscCall(MatSetType(A, MATAIJ)); 33 PetscCall(MatSetBlockSize(A, bs)); 34 PetscCall(MatSeqAIJSetPreallocation(A, PETSC_DEFAULT, NULL)); 35 PetscCall(MatMPIAIJSetPreallocation(A, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL)); 36 PetscCall(MatSetFromOptions(A)); 37 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 38 39 /* Create a BAIJ matrix B */ 40 PetscCall(MatCreate(PETSC_COMM_WORLD, &B)); 41 PetscCall(MatSetSizes(B, m * bs, m * bs, PETSC_DECIDE, PETSC_DECIDE)); 42 PetscCall(MatSetType(B, MATBAIJ)); 43 PetscCall(MatSeqBAIJSetPreallocation(B, bs, PETSC_DEFAULT, NULL)); 44 PetscCall(MatMPIBAIJSetPreallocation(B, bs, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL)); 45 PetscCall(MatSetFromOptions(B)); 46 PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 47 48 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm)); 49 PetscCall(PetscRandomSetFromOptions(rdm)); 50 51 PetscCall(MatGetOwnershipRange(A, &rstart, &rend)); 52 PetscCall(MatGetSize(A, &M, &N)); 53 Mbs = M / bs; 54 55 PetscCall(PetscMalloc1(bs, &rows)); 56 PetscCall(PetscMalloc1(bs, &cols)); 57 PetscCall(PetscMalloc1(bs * bs, &vals)); 58 PetscCall(PetscMalloc1(M, &idx)); 59 60 /* Now set blocks of values */ 61 for (i = 0; i < 40 * bs; i++) { 62 PetscInt nr = 1, nc = 1; 63 PetscCall(PetscRandomGetValue(rdm, &rval)); 64 cols[0] = bs * (int)(PetscRealPart(rval) * Mbs); 65 PetscCall(PetscRandomGetValue(rdm, &rval)); 66 rows[0] = rstart + bs * (int)(PetscRealPart(rval) * m); 67 for (j = 1; j < bs; j++) { 68 PetscCall(PetscRandomGetValue(rdm, &rval)); 69 if (PetscRealPart(rval) > .5) rows[nr++] = rows[0] + j - 1; 70 } 71 for (j = 1; j < bs; j++) { 72 PetscCall(PetscRandomGetValue(rdm, &rval)); 73 if (PetscRealPart(rval) > .5) cols[nc++] = cols[0] + j - 1; 74 } 75 76 for (j = 0; j < nr * nc; j++) { 77 PetscCall(PetscRandomGetValue(rdm, &rval)); 78 vals[j] = rval; 79 } 80 PetscCall(MatSetValues(A, nr, rows, nc, cols, vals, ADD_VALUES)); 81 PetscCall(MatSetValues(B, nr, rows, nc, cols, vals, ADD_VALUES)); 82 } 83 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 84 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 85 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 86 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 87 88 /* Test MatConvert_MPIAIJ_MPIBAIJ handles incompletely filled blocks */ 89 PetscCall(MatConvert(A, MATBAIJ, MAT_INITIAL_MATRIX, &E)); 90 PetscCall(MatDestroy(&E)); 91 92 /* Test MatIncreaseOverlap() */ 93 PetscCall(PetscMalloc1(nd, &is1)); 94 PetscCall(PetscMalloc1(nd, &is2)); 95 96 emptynd = PETSC_FALSE; 97 if (rank == 0 && test_nd0) emptynd = PETSC_TRUE; /* test case */ 98 99 for (i = 0; i < nd; i++) { 100 PetscCall(PetscRandomGetValue(rdm, &rval)); 101 sz = (int)(PetscRealPart(rval) * m); 102 for (j = 0; j < sz; j++) { 103 PetscCall(PetscRandomGetValue(rdm, &rval)); 104 idx[j * bs] = bs * (int)(PetscRealPart(rval) * Mbs); 105 for (k = 1; k < bs; k++) idx[j * bs + k] = idx[j * bs] + k; 106 } 107 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, emptynd ? 0 : sz * bs, idx, PETSC_COPY_VALUES, is1 + i)); 108 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, emptynd ? 0 : sz * bs, idx, PETSC_COPY_VALUES, is2 + i)); 109 } 110 PetscCall(MatIncreaseOverlap(A, nd, is1, ov)); 111 PetscCall(MatIncreaseOverlap(B, nd, is2, ov)); 112 113 for (i = 0; i < nd; ++i) { 114 PetscCall(ISEqual(is1[i], is2[i], &flg)); 115 116 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "i=%" PetscInt_FMT ", flg=%d :bs=%" PetscInt_FMT " m=%" PetscInt_FMT " ov=%" PetscInt_FMT " nd=%" PetscInt_FMT " np=%d\n", i, flg, bs, m, ov, nd, size)); 117 } 118 119 for (i = 0; i < nd; ++i) { 120 PetscCall(ISSort(is1[i])); 121 PetscCall(ISSort(is2[i])); 122 } 123 124 PetscCall(MatCreateSubMatrices(B, nd, is2, is2, MAT_INITIAL_MATRIX, &submatB)); 125 PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_INITIAL_MATRIX, &submatA)); 126 127 /* Test MatMult() */ 128 for (i = 0; i < nd; i++) { 129 PetscCall(MatGetSize(submatA[i], &mm, &nn)); 130 PetscCall(VecCreateSeq(PETSC_COMM_SELF, mm, &xx)); 131 PetscCall(VecDuplicate(xx, &s1)); 132 PetscCall(VecDuplicate(xx, &s2)); 133 for (j = 0; j < 3; j++) { 134 PetscCall(VecSetRandom(xx, rdm)); 135 PetscCall(MatMult(submatA[i], xx, s1)); 136 PetscCall(MatMult(submatB[i], xx, s2)); 137 PetscCall(VecNorm(s1, NORM_2, &s1norm)); 138 PetscCall(VecNorm(s2, NORM_2, &s2norm)); 139 rnorm = s2norm - s1norm; 140 if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", rank, (double)s1norm, (double)s2norm)); 141 } 142 PetscCall(VecDestroy(&xx)); 143 PetscCall(VecDestroy(&s1)); 144 PetscCall(VecDestroy(&s2)); 145 } 146 147 /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */ 148 PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_REUSE_MATRIX, &submatA)); 149 PetscCall(MatCreateSubMatrices(B, nd, is2, is2, MAT_REUSE_MATRIX, &submatB)); 150 151 /* Test MatMult() */ 152 for (i = 0; i < nd; i++) { 153 PetscCall(MatGetSize(submatA[i], &mm, &nn)); 154 PetscCall(VecCreateSeq(PETSC_COMM_SELF, mm, &xx)); 155 PetscCall(VecDuplicate(xx, &s1)); 156 PetscCall(VecDuplicate(xx, &s2)); 157 for (j = 0; j < 3; j++) { 158 PetscCall(VecSetRandom(xx, rdm)); 159 PetscCall(MatMult(submatA[i], xx, s1)); 160 PetscCall(MatMult(submatB[i], xx, s2)); 161 PetscCall(VecNorm(s1, NORM_2, &s1norm)); 162 PetscCall(VecNorm(s2, NORM_2, &s2norm)); 163 rnorm = s2norm - s1norm; 164 if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", rank, (double)s1norm, (double)s2norm)); 165 } 166 PetscCall(VecDestroy(&xx)); 167 PetscCall(VecDestroy(&s1)); 168 PetscCall(VecDestroy(&s2)); 169 } 170 171 /* Free allocated memory */ 172 for (i = 0; i < nd; ++i) { 173 PetscCall(ISDestroy(&is1[i])); 174 PetscCall(ISDestroy(&is2[i])); 175 } 176 PetscCall(MatDestroySubMatrices(nd, &submatA)); 177 PetscCall(MatDestroySubMatrices(nd, &submatB)); 178 179 PetscCall(PetscFree(is1)); 180 PetscCall(PetscFree(is2)); 181 PetscCall(PetscFree(idx)); 182 PetscCall(PetscFree(rows)); 183 PetscCall(PetscFree(cols)); 184 PetscCall(PetscFree(vals)); 185 PetscCall(MatDestroy(&A)); 186 PetscCall(MatDestroy(&B)); 187 PetscCall(PetscRandomDestroy(&rdm)); 188 PetscCall(PetscFinalize()); 189 return 0; 190 } 191 192 /*TEST 193 194 test: 195 nsize: {{1 3}} 196 args: -mat_block_size {{1 3 4 6 8}} -ov {{1 3}} -mat_size {{11 13}} -nd 7 197 output_file: output/ex54.out 198 199 test: 200 suffix: 2 201 args: -nd 2 -test_nd0 202 output_file: output/ex54.out 203 204 test: 205 suffix: 3 206 nsize: 3 207 args: -nd 2 -test_nd0 208 output_file: output/ex54.out 209 210 TEST*/ 211