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