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, Bt, *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_MPI(S)BAIJ handles incompletely filled blocks */ 89 PetscCall(MatConvert(A, MATBAIJ, MAT_INITIAL_MATRIX, &E)); 90 PetscCall(MatDestroy(&E)); 91 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Bt)); 92 PetscCall(MatAXPY(Bt, 1.0, B, DIFFERENT_NONZERO_PATTERN)); 93 PetscCall(MatSetOption(Bt, MAT_SYMMETRIC, PETSC_TRUE)); 94 PetscCall(MatConvert(Bt, MATSBAIJ, MAT_INITIAL_MATRIX, &E)); 95 PetscCall(MatDestroy(&E)); 96 PetscCall(MatDestroy(&Bt)); 97 98 /* Test MatIncreaseOverlap() */ 99 PetscCall(PetscMalloc1(nd, &is1)); 100 PetscCall(PetscMalloc1(nd, &is2)); 101 102 emptynd = PETSC_FALSE; 103 if (rank == 0 && test_nd0) emptynd = PETSC_TRUE; /* test case */ 104 105 for (i = 0; i < nd; i++) { 106 PetscCall(PetscRandomGetValue(rdm, &rval)); 107 sz = (int)(PetscRealPart(rval) * m); 108 for (j = 0; j < sz; j++) { 109 PetscCall(PetscRandomGetValue(rdm, &rval)); 110 idx[j * bs] = bs * (int)(PetscRealPart(rval) * Mbs); 111 for (k = 1; k < bs; k++) idx[j * bs + k] = idx[j * bs] + k; 112 } 113 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, emptynd ? 0 : sz * bs, idx, PETSC_COPY_VALUES, is1 + i)); 114 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, emptynd ? 0 : sz * bs, idx, PETSC_COPY_VALUES, is2 + i)); 115 } 116 PetscCall(MatIncreaseOverlap(A, nd, is1, ov)); 117 PetscCall(MatIncreaseOverlap(B, nd, is2, ov)); 118 119 for (i = 0; i < nd; ++i) { 120 PetscCall(ISEqual(is1[i], is2[i], &flg)); 121 122 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)); 123 } 124 125 for (i = 0; i < nd; ++i) { 126 PetscCall(ISSort(is1[i])); 127 PetscCall(ISSort(is2[i])); 128 } 129 130 PetscCall(MatCreateSubMatrices(B, nd, is2, is2, MAT_INITIAL_MATRIX, &submatB)); 131 PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_INITIAL_MATRIX, &submatA)); 132 133 /* Test MatMult() */ 134 for (i = 0; i < nd; i++) { 135 PetscCall(MatGetSize(submatA[i], &mm, &nn)); 136 PetscCall(VecCreateSeq(PETSC_COMM_SELF, mm, &xx)); 137 PetscCall(VecDuplicate(xx, &s1)); 138 PetscCall(VecDuplicate(xx, &s2)); 139 for (j = 0; j < 3; j++) { 140 PetscCall(VecSetRandom(xx, rdm)); 141 PetscCall(MatMult(submatA[i], xx, s1)); 142 PetscCall(MatMult(submatB[i], xx, s2)); 143 PetscCall(VecNorm(s1, NORM_2, &s1norm)); 144 PetscCall(VecNorm(s2, NORM_2, &s2norm)); 145 rnorm = s2norm - s1norm; 146 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)); 147 } 148 PetscCall(VecDestroy(&xx)); 149 PetscCall(VecDestroy(&s1)); 150 PetscCall(VecDestroy(&s2)); 151 } 152 153 /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */ 154 PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_REUSE_MATRIX, &submatA)); 155 PetscCall(MatCreateSubMatrices(B, nd, is2, is2, MAT_REUSE_MATRIX, &submatB)); 156 157 /* Test MatMult() */ 158 for (i = 0; i < nd; i++) { 159 PetscCall(MatGetSize(submatA[i], &mm, &nn)); 160 PetscCall(VecCreateSeq(PETSC_COMM_SELF, mm, &xx)); 161 PetscCall(VecDuplicate(xx, &s1)); 162 PetscCall(VecDuplicate(xx, &s2)); 163 for (j = 0; j < 3; j++) { 164 PetscCall(VecSetRandom(xx, rdm)); 165 PetscCall(MatMult(submatA[i], xx, s1)); 166 PetscCall(MatMult(submatB[i], xx, s2)); 167 PetscCall(VecNorm(s1, NORM_2, &s1norm)); 168 PetscCall(VecNorm(s2, NORM_2, &s2norm)); 169 rnorm = s2norm - s1norm; 170 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)); 171 } 172 PetscCall(VecDestroy(&xx)); 173 PetscCall(VecDestroy(&s1)); 174 PetscCall(VecDestroy(&s2)); 175 } 176 177 /* Free allocated memory */ 178 for (i = 0; i < nd; ++i) { 179 PetscCall(ISDestroy(&is1[i])); 180 PetscCall(ISDestroy(&is2[i])); 181 } 182 PetscCall(MatDestroySubMatrices(nd, &submatA)); 183 PetscCall(MatDestroySubMatrices(nd, &submatB)); 184 185 PetscCall(PetscFree(is1)); 186 PetscCall(PetscFree(is2)); 187 PetscCall(PetscFree(idx)); 188 PetscCall(PetscFree(rows)); 189 PetscCall(PetscFree(cols)); 190 PetscCall(PetscFree(vals)); 191 PetscCall(MatDestroy(&A)); 192 PetscCall(MatDestroy(&B)); 193 PetscCall(PetscRandomDestroy(&rdm)); 194 PetscCall(PetscFinalize()); 195 return 0; 196 } 197 198 /*TEST 199 200 test: 201 nsize: {{1 3}} 202 args: -mat_block_size {{1 3 4 6 8}} -ov {{1 3}} -mat_size {{11 13}} -nd 7 203 output_file: output/ex54.out 204 205 test: 206 suffix: 2 207 args: -nd 2 -test_nd0 208 output_file: output/ex54.out 209 210 test: 211 suffix: 3 212 nsize: 3 213 args: -nd 2 -test_nd0 214 output_file: output/ex54.out 215 216 TEST*/ 217