1 2 static char help[] = "Tests various routines in MatMPISBAIJ format.\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc, char **args) { 7 Vec x, y, u, s1, s2; 8 Mat A, sA, sB; 9 PetscRandom rctx; 10 PetscReal r1, r2, rnorm, tol = PETSC_SQRT_MACHINE_EPSILON; 11 PetscScalar one = 1.0, neg_one = -1.0, value[3], four = 4.0, alpha = 0.1; 12 PetscInt n, col[3], n1, block, row, i, j, i2, j2, Ii, J, rstart, rend, bs = 1, mbs = 16, d_nz = 3, o_nz = 3, prob = 1; 13 PetscMPIInt size, rank; 14 PetscBool flg; 15 16 PetscFunctionBeginUser; 17 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 18 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mbs", &mbs, NULL)); 19 PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL)); 20 PetscCall(PetscOptionsGetInt(NULL, NULL, "-prob", &prob, NULL)); 21 22 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 23 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 24 25 /* Create a BAIJ matrix A */ 26 n = mbs * bs; 27 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 28 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n)); 29 PetscCall(MatSetType(A, MATBAIJ)); 30 PetscCall(MatSetFromOptions(A)); 31 PetscCall(MatMPIBAIJSetPreallocation(A, bs, d_nz, NULL, o_nz, NULL)); 32 PetscCall(MatSeqBAIJSetPreallocation(A, bs, d_nz, NULL)); 33 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 34 35 if (bs == 1) { 36 if (prob == 1) { /* tridiagonal matrix */ 37 value[0] = -1.0; 38 value[1] = 0.0; 39 value[2] = -1.0; 40 for (i = 1; i < n - 1; i++) { 41 col[0] = i - 1; 42 col[1] = i; 43 col[2] = i + 1; 44 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES)); 45 } 46 i = n - 1; 47 col[0] = 0; 48 col[1] = n - 2; 49 col[2] = n - 1; 50 value[0] = 0.1; 51 value[1] = -1.0; 52 value[2] = 0.0; 53 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES)); 54 55 i = 0; 56 col[0] = 0; 57 col[1] = 1; 58 col[2] = n - 1; 59 value[0] = 0.0; 60 value[1] = -1.0; 61 value[2] = 0.1; 62 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES)); 63 } else if (prob == 2) { /* matrix for the five point stencil */ 64 n1 = (int)PetscSqrtReal((PetscReal)n); 65 for (i = 0; i < n1; i++) { 66 for (j = 0; j < n1; j++) { 67 Ii = j + n1 * i; 68 if (i > 0) { 69 J = Ii - n1; 70 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 71 } 72 if (i < n1 - 1) { 73 J = Ii + n1; 74 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 75 } 76 if (j > 0) { 77 J = Ii - 1; 78 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 79 } 80 if (j < n1 - 1) { 81 J = Ii + 1; 82 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 83 } 84 PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &four, INSERT_VALUES)); 85 } 86 } 87 } 88 /* end of if (bs == 1) */ 89 } else { /* bs > 1 */ 90 for (block = 0; block < n / bs; block++) { 91 /* diagonal blocks */ 92 value[0] = -1.0; 93 value[1] = 4.0; 94 value[2] = -1.0; 95 for (i = 1 + block * bs; i < bs - 1 + block * bs; i++) { 96 col[0] = i - 1; 97 col[1] = i; 98 col[2] = i + 1; 99 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES)); 100 } 101 i = bs - 1 + block * bs; 102 col[0] = bs - 2 + block * bs; 103 col[1] = bs - 1 + block * bs; 104 value[0] = -1.0; 105 value[1] = 4.0; 106 PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES)); 107 108 i = 0 + block * bs; 109 col[0] = 0 + block * bs; 110 col[1] = 1 + block * bs; 111 value[0] = 4.0; 112 value[1] = -1.0; 113 PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES)); 114 } 115 /* off-diagonal blocks */ 116 value[0] = -1.0; 117 for (i = 0; i < (n / bs - 1) * bs; i++) { 118 col[0] = i + bs; 119 PetscCall(MatSetValues(A, 1, &i, 1, col, value, INSERT_VALUES)); 120 col[0] = i; 121 row = i + bs; 122 PetscCall(MatSetValues(A, 1, &row, 1, col, value, INSERT_VALUES)); 123 } 124 } 125 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 126 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 127 PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE)); 128 129 /* Get SBAIJ matrix sA from A */ 130 PetscCall(MatConvert(A, MATSBAIJ, MAT_INITIAL_MATRIX, &sA)); 131 132 /* Test MatGetSize(), MatGetLocalSize() */ 133 PetscCall(MatGetSize(sA, &i, &j)); 134 PetscCall(MatGetSize(A, &i2, &j2)); 135 i -= i2; 136 j -= j2; 137 if (i || j) { 138 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatGetSize()\n", rank)); 139 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 140 } 141 142 PetscCall(MatGetLocalSize(sA, &i, &j)); 143 PetscCall(MatGetLocalSize(A, &i2, &j2)); 144 i2 -= i; 145 j2 -= j; 146 if (i2 || j2) { 147 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatGetLocalSize()\n", rank)); 148 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 149 } 150 151 /* vectors */ 152 /*--------------------*/ 153 /* i is obtained from MatGetLocalSize() */ 154 PetscCall(VecCreate(PETSC_COMM_WORLD, &x)); 155 PetscCall(VecSetSizes(x, i, PETSC_DECIDE)); 156 PetscCall(VecSetFromOptions(x)); 157 PetscCall(VecDuplicate(x, &y)); 158 PetscCall(VecDuplicate(x, &u)); 159 PetscCall(VecDuplicate(x, &s1)); 160 PetscCall(VecDuplicate(x, &s2)); 161 162 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx)); 163 PetscCall(PetscRandomSetFromOptions(rctx)); 164 PetscCall(VecSetRandom(x, rctx)); 165 PetscCall(VecSet(u, one)); 166 167 /* Test MatNorm() */ 168 PetscCall(MatNorm(A, NORM_FROBENIUS, &r1)); 169 PetscCall(MatNorm(sA, NORM_FROBENIUS, &r2)); 170 rnorm = PetscAbsReal(r1 - r2) / r2; 171 if (rnorm > tol && rank == 0) { PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm_FROBENIUS(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n", (double)r1, (double)r2, bs)); } 172 PetscCall(MatNorm(A, NORM_INFINITY, &r1)); 173 PetscCall(MatNorm(sA, NORM_INFINITY, &r2)); 174 rnorm = PetscAbsReal(r1 - r2) / r2; 175 if (rnorm > tol && rank == 0) { PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: MatNorm_INFINITY(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n", (double)r1, (double)r2, bs)); } 176 PetscCall(MatNorm(A, NORM_1, &r1)); 177 PetscCall(MatNorm(sA, NORM_1, &r2)); 178 rnorm = PetscAbsReal(r1 - r2) / r2; 179 if (rnorm > tol && rank == 0) { PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: MatNorm_1(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n", (double)r1, (double)r2, bs)); } 180 181 /* Test MatGetOwnershipRange() */ 182 PetscCall(MatGetOwnershipRange(sA, &rstart, &rend)); 183 PetscCall(MatGetOwnershipRange(A, &i2, &j2)); 184 i2 -= rstart; 185 j2 -= rend; 186 if (i2 || j2) { 187 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MaGetOwnershipRange()\n", rank)); 188 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 189 } 190 191 /* Test MatDiagonalScale() */ 192 PetscCall(MatDiagonalScale(A, x, x)); 193 PetscCall(MatDiagonalScale(sA, x, x)); 194 PetscCall(MatMultEqual(A, sA, 10, &flg)); 195 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "Error in MatDiagonalScale"); 196 197 /* Test MatGetDiagonal(), MatScale() */ 198 PetscCall(MatGetDiagonal(A, s1)); 199 PetscCall(MatGetDiagonal(sA, s2)); 200 PetscCall(VecNorm(s1, NORM_1, &r1)); 201 PetscCall(VecNorm(s2, NORM_1, &r2)); 202 r1 -= r2; 203 if (r1 < -tol || r1 > tol) { 204 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatDiagonalScale() or MatGetDiagonal(), r1=%g \n", rank, (double)r1)); 205 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 206 } 207 208 PetscCall(MatScale(A, alpha)); 209 PetscCall(MatScale(sA, alpha)); 210 211 /* Test MatGetRowMaxAbs() */ 212 PetscCall(MatGetRowMaxAbs(A, s1, NULL)); 213 PetscCall(MatGetRowMaxAbs(sA, s2, NULL)); 214 215 PetscCall(VecNorm(s1, NORM_1, &r1)); 216 PetscCall(VecNorm(s2, NORM_1, &r2)); 217 r1 -= r2; 218 if (r1 < -tol || r1 > tol) { PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetRowMaxAbs() \n")); } 219 220 /* Test MatMult(), MatMultAdd() */ 221 PetscCall(MatMultEqual(A, sA, 10, &flg)); 222 if (!flg) { 223 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatMult() or MatScale()\n", rank)); 224 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 225 } 226 227 PetscCall(MatMultAddEqual(A, sA, 10, &flg)); 228 if (!flg) { 229 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatMultAdd()\n", rank)); 230 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 231 } 232 233 /* Test MatMultTranspose(), MatMultTransposeAdd() */ 234 for (i = 0; i < 10; i++) { 235 PetscCall(VecSetRandom(x, rctx)); 236 PetscCall(MatMultTranspose(A, x, s1)); 237 PetscCall(MatMultTranspose(sA, x, s2)); 238 PetscCall(VecNorm(s1, NORM_1, &r1)); 239 PetscCall(VecNorm(s2, NORM_1, &r2)); 240 r1 -= r2; 241 if (r1 < -tol || r1 > tol) { 242 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatMult() or MatScale(), err=%g\n", rank, (double)r1)); 243 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 244 } 245 } 246 for (i = 0; i < 10; i++) { 247 PetscCall(VecSetRandom(x, rctx)); 248 PetscCall(VecSetRandom(y, rctx)); 249 PetscCall(MatMultTransposeAdd(A, x, y, s1)); 250 PetscCall(MatMultTransposeAdd(sA, x, y, s2)); 251 PetscCall(VecNorm(s1, NORM_1, &r1)); 252 PetscCall(VecNorm(s2, NORM_1, &r2)); 253 r1 -= r2; 254 if (r1 < -tol || r1 > tol) { 255 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatMultAdd(), err=%g \n", rank, (double)r1)); 256 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 257 } 258 } 259 260 /* Test MatDuplicate() */ 261 PetscCall(MatDuplicate(sA, MAT_COPY_VALUES, &sB)); 262 PetscCall(MatEqual(sA, sB, &flg)); 263 if (!flg) { PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Error in MatDuplicate(), sA != sB \n")); } 264 PetscCall(MatMultEqual(sA, sB, 5, &flg)); 265 if (!flg) { 266 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatDuplicate() or MatMult()\n", rank)); 267 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 268 } 269 PetscCall(MatMultAddEqual(sA, sB, 5, &flg)); 270 if (!flg) { 271 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatDuplicate() or MatMultAdd(()\n", rank)); 272 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT)); 273 } 274 PetscCall(MatDestroy(&sB)); 275 PetscCall(VecDestroy(&u)); 276 PetscCall(VecDestroy(&x)); 277 PetscCall(VecDestroy(&y)); 278 PetscCall(VecDestroy(&s1)); 279 PetscCall(VecDestroy(&s2)); 280 PetscCall(MatDestroy(&sA)); 281 PetscCall(MatDestroy(&A)); 282 PetscCall(PetscRandomDestroy(&rctx)); 283 PetscCall(PetscFinalize()); 284 return 0; 285 } 286 287 /*TEST 288 289 test: 290 nsize: {{1 3}} 291 args: -bs {{1 2 3 5 7 8}} -mat_ignore_lower_triangular -prob {{1 2}} 292 293 TEST*/ 294