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