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