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