1 2 static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for MatBAIJ format.\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 = 43, ov = 1, i, j, k, *rows, *cols, M, nd = 5, *idx, mm, nn, lsize; 9 PetscScalar *vals, rval; 10 IS *is1, *is2; 11 PetscRandom rdm; 12 Vec xx, s1, s2; 13 PetscReal s1norm, s2norm, rnorm, tol = PETSC_SQRT_MACHINE_EPSILON; 14 PetscBool flg; 15 16 PetscFunctionBeginUser; 17 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 18 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL)); 19 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_size", &m, NULL)); 20 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ov", &ov, NULL)); 21 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nd", &nd, NULL)); 22 M = m * bs; 23 24 PetscCall(MatCreateSeqBAIJ(PETSC_COMM_SELF, bs, M, M, 1, NULL, &A)); 25 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 26 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, M, M, 15, NULL, &B)); 27 PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 28 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm)); 29 PetscCall(PetscRandomSetFromOptions(rdm)); 30 31 PetscCall(PetscMalloc1(bs, &rows)); 32 PetscCall(PetscMalloc1(bs, &cols)); 33 PetscCall(PetscMalloc1(bs * bs, &vals)); 34 PetscCall(PetscMalloc1(M, &idx)); 35 36 /* Now set blocks of values */ 37 for (i = 0; i < 20 * bs; i++) { 38 PetscCall(PetscRandomGetValue(rdm, &rval)); 39 cols[0] = bs * (int)(PetscRealPart(rval) * m); 40 PetscCall(PetscRandomGetValue(rdm, &rval)); 41 rows[0] = bs * (int)(PetscRealPart(rval) * m); 42 for (j = 1; j < bs; j++) { 43 rows[j] = rows[j - 1] + 1; 44 cols[j] = cols[j - 1] + 1; 45 } 46 47 for (j = 0; j < bs * bs; j++) { 48 PetscCall(PetscRandomGetValue(rdm, &rval)); 49 vals[j] = rval; 50 } 51 PetscCall(MatSetValues(A, bs, rows, bs, cols, vals, ADD_VALUES)); 52 PetscCall(MatSetValues(B, bs, rows, bs, cols, vals, ADD_VALUES)); 53 } 54 55 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 56 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 57 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 58 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 59 60 /* Test MatIncreaseOverlap() */ 61 PetscCall(PetscMalloc1(nd, &is1)); 62 PetscCall(PetscMalloc1(nd, &is2)); 63 64 for (i = 0; i < nd; i++) { 65 PetscCall(PetscRandomGetValue(rdm, &rval)); 66 lsize = (int)(PetscRealPart(rval) * m); 67 for (j = 0; j < lsize; j++) { 68 PetscCall(PetscRandomGetValue(rdm, &rval)); 69 idx[j * bs] = bs * (int)(PetscRealPart(rval) * m); 70 for (k = 1; k < bs; k++) idx[j * bs + k] = idx[j * bs] + k; 71 } 72 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, lsize * bs, idx, PETSC_COPY_VALUES, is1 + i)); 73 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, lsize * bs, idx, PETSC_COPY_VALUES, is2 + i)); 74 } 75 PetscCall(MatIncreaseOverlap(A, nd, is1, ov)); 76 PetscCall(MatIncreaseOverlap(B, nd, is2, ov)); 77 78 for (i = 0; i < nd; ++i) { 79 PetscCall(ISEqual(is1[i], is2[i], &flg)); 80 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "i=%" PetscInt_FMT ", flg =%d", i, (int)flg); 81 } 82 83 for (i = 0; i < nd; ++i) { 84 PetscCall(ISSort(is1[i])); 85 PetscCall(ISSort(is2[i])); 86 } 87 88 PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_INITIAL_MATRIX, &submatA)); 89 PetscCall(MatCreateSubMatrices(B, nd, is2, is2, MAT_INITIAL_MATRIX, &submatB)); 90 91 /* Test MatMult() */ 92 for (i = 0; i < nd; i++) { 93 PetscCall(MatGetSize(submatA[i], &mm, &nn)); 94 PetscCall(VecCreateSeq(PETSC_COMM_SELF, mm, &xx)); 95 PetscCall(VecDuplicate(xx, &s1)); 96 PetscCall(VecDuplicate(xx, &s2)); 97 for (j = 0; j < 3; j++) { 98 PetscCall(VecSetRandom(xx, rdm)); 99 PetscCall(MatMult(submatA[i], xx, s1)); 100 PetscCall(MatMult(submatB[i], xx, s2)); 101 PetscCall(VecNorm(s1, NORM_2, &s1norm)); 102 PetscCall(VecNorm(s2, NORM_2, &s2norm)); 103 rnorm = s2norm - s1norm; 104 if (rnorm < -tol || rnorm > tol) { PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", (double)s1norm, (double)s2norm)); } 105 } 106 PetscCall(VecDestroy(&xx)); 107 PetscCall(VecDestroy(&s1)); 108 PetscCall(VecDestroy(&s2)); 109 } 110 /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */ 111 PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_REUSE_MATRIX, &submatA)); 112 PetscCall(MatCreateSubMatrices(B, nd, is2, is2, MAT_REUSE_MATRIX, &submatB)); 113 114 /* Test MatMult() */ 115 for (i = 0; i < nd; i++) { 116 PetscCall(MatGetSize(submatA[i], &mm, &nn)); 117 PetscCall(VecCreateSeq(PETSC_COMM_SELF, mm, &xx)); 118 PetscCall(VecDuplicate(xx, &s1)); 119 PetscCall(VecDuplicate(xx, &s2)); 120 for (j = 0; j < 3; j++) { 121 PetscCall(VecSetRandom(xx, rdm)); 122 PetscCall(MatMult(submatA[i], xx, s1)); 123 PetscCall(MatMult(submatB[i], xx, s2)); 124 PetscCall(VecNorm(s1, NORM_2, &s1norm)); 125 PetscCall(VecNorm(s2, NORM_2, &s2norm)); 126 rnorm = s2norm - s1norm; 127 if (rnorm < -tol || rnorm > tol) { PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", (double)s1norm, (double)s2norm)); } 128 } 129 PetscCall(VecDestroy(&xx)); 130 PetscCall(VecDestroy(&s1)); 131 PetscCall(VecDestroy(&s2)); 132 } 133 134 /* Free allocated memory */ 135 for (i = 0; i < nd; ++i) { 136 PetscCall(ISDestroy(&is1[i])); 137 PetscCall(ISDestroy(&is2[i])); 138 } 139 PetscCall(MatDestroySubMatrices(nd, &submatA)); 140 PetscCall(MatDestroySubMatrices(nd, &submatB)); 141 PetscCall(PetscFree(is1)); 142 PetscCall(PetscFree(is2)); 143 PetscCall(PetscFree(idx)); 144 PetscCall(PetscFree(rows)); 145 PetscCall(PetscFree(cols)); 146 PetscCall(PetscFree(vals)); 147 PetscCall(MatDestroy(&A)); 148 PetscCall(MatDestroy(&B)); 149 PetscCall(PetscRandomDestroy(&rdm)); 150 PetscCall(PetscFinalize()); 151 return 0; 152 } 153 154 /*TEST 155 156 test: 157 args: -mat_block_size {{1 2 5 7 8}} -ov {{1 3}} -mat_size {{11 13}} -nd {{7}} 158 159 TEST*/ 160