1 2 static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for sequential MatSBAIJ format. Derived from ex51.c\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc, char **args) 7 { 8 Mat A, Atrans, sA, *submatA, *submatsA; 9 PetscInt bs = 1, m = 43, ov = 1, i, j, k, *rows, *cols, M, nd = 5, *idx, mm, nn; 10 PetscMPIInt size; 11 PetscScalar *vals, rval, one = 1.0; 12 IS *is1, *is2; 13 PetscRandom rand; 14 Vec xx, s1, s2; 15 PetscReal s1norm, s2norm, rnorm, tol = 10 * PETSC_SMALL; 16 PetscBool flg; 17 18 PetscFunctionBeginUser; 19 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 20 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL)); 21 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_size", &m, NULL)); 22 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ov", &ov, NULL)); 23 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nd", &nd, NULL)); 24 25 /* create a SeqBAIJ matrix A */ 26 M = m * bs; 27 PetscCall(MatCreateSeqBAIJ(PETSC_COMM_SELF, bs, M, M, 1, NULL, &A)); 28 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 29 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rand)); 30 PetscCall(PetscRandomSetFromOptions(rand)); 31 32 PetscCall(PetscMalloc1(bs, &rows)); 33 PetscCall(PetscMalloc1(bs, &cols)); 34 PetscCall(PetscMalloc1(bs * bs, &vals)); 35 PetscCall(PetscMalloc1(M, &idx)); 36 37 /* Now set blocks of random values */ 38 /* first, set diagonal blocks as zero */ 39 for (j = 0; j < bs * bs; j++) vals[j] = 0.0; 40 for (i = 0; i < m; i++) { 41 cols[0] = i * bs; 42 rows[0] = i * bs; 43 for (j = 1; j < bs; j++) { 44 rows[j] = rows[j - 1] + 1; 45 cols[j] = cols[j - 1] + 1; 46 } 47 PetscCall(MatSetValues(A, bs, rows, bs, cols, vals, ADD_VALUES)); 48 } 49 /* second, add random blocks */ 50 for (i = 0; i < 20 * bs; i++) { 51 PetscCall(PetscRandomGetValue(rand, &rval)); 52 cols[0] = bs * (int)(PetscRealPart(rval) * m); 53 PetscCall(PetscRandomGetValue(rand, &rval)); 54 rows[0] = bs * (int)(PetscRealPart(rval) * m); 55 for (j = 1; j < bs; j++) { 56 rows[j] = rows[j - 1] + 1; 57 cols[j] = cols[j - 1] + 1; 58 } 59 60 for (j = 0; j < bs * bs; j++) { 61 PetscCall(PetscRandomGetValue(rand, &rval)); 62 vals[j] = rval; 63 } 64 PetscCall(MatSetValues(A, bs, rows, bs, cols, vals, ADD_VALUES)); 65 } 66 67 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 68 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 69 70 /* make A a symmetric matrix: A <- A^T + A */ 71 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Atrans)); 72 PetscCall(MatAXPY(A, one, Atrans, DIFFERENT_NONZERO_PATTERN)); 73 PetscCall(MatDestroy(&Atrans)); 74 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Atrans)); 75 PetscCall(MatEqual(A, Atrans, &flg)); 76 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "A+A^T is non-symmetric"); 77 PetscCall(MatDestroy(&Atrans)); 78 79 /* create a SeqSBAIJ matrix sA (= A) */ 80 PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE)); 81 PetscCall(MatConvert(A, MATSEQSBAIJ, MAT_INITIAL_MATRIX, &sA)); 82 83 /* Test sA==A through MatMult() */ 84 for (i = 0; i < nd; i++) { 85 PetscCall(MatGetSize(A, &mm, &nn)); 86 PetscCall(VecCreateSeq(PETSC_COMM_SELF, mm, &xx)); 87 PetscCall(VecDuplicate(xx, &s1)); 88 PetscCall(VecDuplicate(xx, &s2)); 89 for (j = 0; j < 3; j++) { 90 PetscCall(VecSetRandom(xx, rand)); 91 PetscCall(MatMult(A, xx, s1)); 92 PetscCall(MatMult(sA, xx, s2)); 93 PetscCall(VecNorm(s1, NORM_2, &s1norm)); 94 PetscCall(VecNorm(s2, NORM_2, &s2norm)); 95 rnorm = s2norm - s1norm; 96 if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", (double)s1norm, (double)s2norm)); 97 } 98 PetscCall(VecDestroy(&xx)); 99 PetscCall(VecDestroy(&s1)); 100 PetscCall(VecDestroy(&s2)); 101 } 102 103 /* Test MatIncreaseOverlap() */ 104 PetscCall(PetscMalloc1(nd, &is1)); 105 PetscCall(PetscMalloc1(nd, &is2)); 106 107 for (i = 0; i < nd; i++) { 108 PetscCall(PetscRandomGetValue(rand, &rval)); 109 size = (int)(PetscRealPart(rval) * m); 110 for (j = 0; j < size; j++) { 111 PetscCall(PetscRandomGetValue(rand, &rval)); 112 idx[j * bs] = bs * (int)(PetscRealPart(rval) * m); 113 for (k = 1; k < bs; k++) idx[j * bs + k] = idx[j * bs] + k; 114 } 115 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size * bs, idx, PETSC_COPY_VALUES, is1 + i)); 116 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size * bs, idx, PETSC_COPY_VALUES, is2 + i)); 117 } 118 /* for debugging */ 119 /* 120 PetscCall(MatView(A,PETSC_VIEWER_STDOUT_SELF)); 121 PetscCall(MatView(sA,PETSC_VIEWER_STDOUT_SELF)); 122 */ 123 124 PetscCall(MatIncreaseOverlap(A, nd, is1, ov)); 125 PetscCall(MatIncreaseOverlap(sA, nd, is2, ov)); 126 127 for (i = 0; i < nd; ++i) { 128 PetscCall(ISSort(is1[i])); 129 PetscCall(ISSort(is2[i])); 130 } 131 132 for (i = 0; i < nd; ++i) { 133 PetscCall(ISEqual(is1[i], is2[i], &flg)); 134 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "i=%" PetscInt_FMT ", is1 != is2", i); 135 } 136 137 PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_INITIAL_MATRIX, &submatA)); 138 PetscCall(MatCreateSubMatrices(sA, nd, is2, is2, MAT_INITIAL_MATRIX, &submatsA)); 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, rand)); 148 PetscCall(MatMult(submatA[i], xx, s1)); 149 PetscCall(MatMult(submatsA[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, "Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", (double)s1norm, (double)s2norm)); 154 } 155 PetscCall(VecDestroy(&xx)); 156 PetscCall(VecDestroy(&s1)); 157 PetscCall(VecDestroy(&s2)); 158 } 159 160 /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */ 161 PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_REUSE_MATRIX, &submatA)); 162 PetscCall(MatCreateSubMatrices(sA, nd, is2, is2, MAT_REUSE_MATRIX, &submatsA)); 163 164 /* Test MatMult() */ 165 for (i = 0; i < nd; i++) { 166 PetscCall(MatGetSize(submatA[i], &mm, &nn)); 167 PetscCall(VecCreateSeq(PETSC_COMM_SELF, mm, &xx)); 168 PetscCall(VecDuplicate(xx, &s1)); 169 PetscCall(VecDuplicate(xx, &s2)); 170 for (j = 0; j < 3; j++) { 171 PetscCall(VecSetRandom(xx, rand)); 172 PetscCall(MatMult(submatA[i], xx, s1)); 173 PetscCall(MatMult(submatsA[i], xx, s2)); 174 PetscCall(VecNorm(s1, NORM_2, &s1norm)); 175 PetscCall(VecNorm(s2, NORM_2, &s2norm)); 176 rnorm = s2norm - s1norm; 177 if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", (double)s1norm, (double)s2norm)); 178 } 179 PetscCall(VecDestroy(&xx)); 180 PetscCall(VecDestroy(&s1)); 181 PetscCall(VecDestroy(&s2)); 182 } 183 184 /* Free allocated memory */ 185 for (i = 0; i < nd; ++i) { 186 PetscCall(ISDestroy(&is1[i])); 187 PetscCall(ISDestroy(&is2[i])); 188 } 189 PetscCall(MatDestroySubMatrices(nd, &submatA)); 190 PetscCall(MatDestroySubMatrices(nd, &submatsA)); 191 192 PetscCall(PetscFree(is1)); 193 PetscCall(PetscFree(is2)); 194 PetscCall(PetscFree(idx)); 195 PetscCall(PetscFree(rows)); 196 PetscCall(PetscFree(cols)); 197 PetscCall(PetscFree(vals)); 198 PetscCall(MatDestroy(&A)); 199 PetscCall(MatDestroy(&sA)); 200 PetscCall(PetscRandomDestroy(&rand)); 201 PetscCall(PetscFinalize()); 202 return 0; 203 } 204 205 /*TEST 206 207 test: 208 args: -ov 2 209 210 TEST*/ 211