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