1 static char help[] = "Tests the various sequential routines in MATSEQSBAIJ format.\n"; 2 3 #include <petscmat.h> 4 5 int main(int argc, char **args) 6 { 7 PetscMPIInt size; 8 Vec x, y, b, s1, s2; 9 Mat A; /* linear system matrix */ 10 Mat sA, sB, sFactor, B, C; /* symmetric matrices */ 11 PetscInt n, mbs = 16, bs = 1, nz = 3, prob = 1, i, j, k1, k2, col[3], lf, block, row, Ii, J, n1, inc; 12 PetscReal norm1, norm2, rnorm, tol = 10 * PETSC_SMALL; 13 PetscScalar neg_one = -1.0, four = 4.0, value[3]; 14 IS perm, iscol; 15 PetscRandom rdm; 16 PetscBool doIcc = PETSC_TRUE, equal; 17 MatInfo minfo1, minfo2; 18 MatFactorInfo factinfo; 19 MatType type; 20 21 PetscFunctionBeginUser; 22 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 23 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 24 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); 25 PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL)); 26 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mbs", &mbs, NULL)); 27 28 n = mbs * bs; 29 PetscCall(MatCreate(PETSC_COMM_SELF, &A)); 30 PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE)); 31 PetscCall(MatSetType(A, MATSEQBAIJ)); 32 PetscCall(MatSetFromOptions(A)); 33 PetscCall(MatSeqBAIJSetPreallocation(A, bs, nz, NULL)); 34 35 PetscCall(MatCreate(PETSC_COMM_SELF, &sA)); 36 PetscCall(MatSetSizes(sA, n, n, PETSC_DETERMINE, PETSC_DETERMINE)); 37 PetscCall(MatSetType(sA, MATSEQSBAIJ)); 38 PetscCall(MatSetFromOptions(sA)); 39 PetscCall(MatGetType(sA, &type)); 40 PetscCall(PetscObjectTypeCompare((PetscObject)sA, MATSEQSBAIJ, &doIcc)); 41 PetscCall(MatSeqSBAIJSetPreallocation(sA, bs, nz, NULL)); 42 PetscCall(MatSetOption(sA, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE)); 43 44 /* Test MatGetOwnershipRange() */ 45 PetscCall(MatGetOwnershipRange(A, &Ii, &J)); 46 PetscCall(MatGetOwnershipRange(sA, &i, &j)); 47 if (i - Ii || j - J) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetOwnershipRange() in MatSBAIJ format\n")); 48 49 /* Assemble matrix */ 50 if (bs == 1) { 51 PetscCall(PetscOptionsGetInt(NULL, NULL, "-test_problem", &prob, NULL)); 52 if (prob == 1) { /* tridiagonal matrix */ 53 value[0] = -1.0; 54 value[1] = 2.0; 55 value[2] = -1.0; 56 for (i = 1; i < n - 1; i++) { 57 col[0] = i - 1; 58 col[1] = i; 59 col[2] = i + 1; 60 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES)); 61 PetscCall(MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES)); 62 } 63 i = n - 1; 64 col[0] = 0; 65 col[1] = n - 2; 66 col[2] = n - 1; 67 68 value[0] = 0.1; 69 value[1] = -1; 70 value[2] = 2; 71 72 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES)); 73 PetscCall(MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES)); 74 75 i = 0; 76 col[0] = n - 1; 77 col[1] = 1; 78 col[2] = 0; 79 value[0] = 0.1; 80 value[1] = -1.0; 81 value[2] = 2; 82 83 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES)); 84 PetscCall(MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES)); 85 86 } else if (prob == 2) { /* matrix for the five point stencil */ 87 n1 = (PetscInt)(PetscSqrtReal((PetscReal)n) + 0.001); 88 PetscCheck(n1 * n1 == n, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "sqrt(n) must be a positive integer!"); 89 for (i = 0; i < n1; i++) { 90 for (j = 0; j < n1; j++) { 91 Ii = j + n1 * i; 92 if (i > 0) { 93 J = Ii - n1; 94 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 95 PetscCall(MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 96 } 97 if (i < n1 - 1) { 98 J = Ii + n1; 99 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 100 PetscCall(MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 101 } 102 if (j > 0) { 103 J = Ii - 1; 104 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 105 PetscCall(MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 106 } 107 if (j < n1 - 1) { 108 J = Ii + 1; 109 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 110 PetscCall(MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 111 } 112 PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &four, INSERT_VALUES)); 113 PetscCall(MatSetValues(sA, 1, &Ii, 1, &Ii, &four, INSERT_VALUES)); 114 } 115 } 116 } 117 118 } else { /* bs > 1 */ 119 for (block = 0; block < n / bs; block++) { 120 /* diagonal blocks */ 121 value[0] = -1.0; 122 value[1] = 4.0; 123 value[2] = -1.0; 124 for (i = 1 + block * bs; i < bs - 1 + block * bs; i++) { 125 col[0] = i - 1; 126 col[1] = i; 127 col[2] = i + 1; 128 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES)); 129 PetscCall(MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES)); 130 } 131 i = bs - 1 + block * bs; 132 col[0] = bs - 2 + block * bs; 133 col[1] = bs - 1 + block * bs; 134 135 value[0] = -1.0; 136 value[1] = 4.0; 137 138 PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES)); 139 PetscCall(MatSetValues(sA, 1, &i, 2, col, value, INSERT_VALUES)); 140 141 i = 0 + block * bs; 142 col[0] = 0 + block * bs; 143 col[1] = 1 + block * bs; 144 145 value[0] = 4.0; 146 value[1] = -1.0; 147 148 PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES)); 149 PetscCall(MatSetValues(sA, 1, &i, 2, col, value, INSERT_VALUES)); 150 } 151 /* off-diagonal blocks */ 152 value[0] = -1.0; 153 for (i = 0; i < (n / bs - 1) * bs; i++) { 154 col[0] = i + bs; 155 156 PetscCall(MatSetValues(A, 1, &i, 1, col, value, INSERT_VALUES)); 157 PetscCall(MatSetValues(sA, 1, &i, 1, col, value, INSERT_VALUES)); 158 159 col[0] = i; 160 row = i + bs; 161 162 PetscCall(MatSetValues(A, 1, &row, 1, col, value, INSERT_VALUES)); 163 PetscCall(MatSetValues(sA, 1, &row, 1, col, value, INSERT_VALUES)); 164 } 165 } 166 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 167 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 168 169 PetscCall(MatAssemblyBegin(sA, MAT_FINAL_ASSEMBLY)); 170 PetscCall(MatAssemblyEnd(sA, MAT_FINAL_ASSEMBLY)); 171 172 /* Test MatGetInfo() of A and sA */ 173 PetscCall(MatGetInfo(A, MAT_LOCAL, &minfo1)); 174 PetscCall(MatGetInfo(sA, MAT_LOCAL, &minfo2)); 175 i = (int)(minfo1.nz_used - minfo2.nz_used); 176 j = (int)(minfo1.nz_allocated - minfo2.nz_allocated); 177 k1 = (int)(minfo1.nz_allocated - minfo1.nz_used); 178 k2 = (int)(minfo2.nz_allocated - minfo2.nz_used); 179 if (i < 0 || j < 0 || k1 < 0 || k2 < 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error (compare A and sA): MatGetInfo()\n")); 180 181 /* Test MatDuplicate() */ 182 PetscCall(MatNorm(A, NORM_FROBENIUS, &norm1)); 183 PetscCall(MatDuplicate(sA, MAT_COPY_VALUES, &sB)); 184 PetscCall(MatEqual(sA, sB, &equal)); 185 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "Error in MatDuplicate()"); 186 187 /* Test MatNorm() */ 188 PetscCall(MatNorm(A, NORM_FROBENIUS, &norm1)); 189 PetscCall(MatNorm(sB, NORM_FROBENIUS, &norm2)); 190 rnorm = PetscAbsReal(norm1 - norm2) / norm2; 191 if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm_FROBENIUS, NormA=%16.14e NormsB=%16.14e\n", (double)norm1, (double)norm2)); 192 PetscCall(MatNorm(A, NORM_INFINITY, &norm1)); 193 PetscCall(MatNorm(sB, NORM_INFINITY, &norm2)); 194 rnorm = PetscAbsReal(norm1 - norm2) / norm2; 195 if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm_INFINITY(), NormA=%16.14e NormsB=%16.14e\n", (double)norm1, (double)norm2)); 196 PetscCall(MatNorm(A, NORM_1, &norm1)); 197 PetscCall(MatNorm(sB, NORM_1, &norm2)); 198 rnorm = PetscAbsReal(norm1 - norm2) / norm2; 199 if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm_INFINITY(), NormA=%16.14e NormsB=%16.14e\n", (double)norm1, (double)norm2)); 200 201 /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */ 202 PetscCall(MatGetInfo(A, MAT_LOCAL, &minfo1)); 203 PetscCall(MatGetInfo(sB, MAT_LOCAL, &minfo2)); 204 i = (int)(minfo1.nz_used - minfo2.nz_used); 205 j = (int)(minfo1.nz_allocated - minfo2.nz_allocated); 206 k1 = (int)(minfo1.nz_allocated - minfo1.nz_used); 207 k2 = (int)(minfo2.nz_allocated - minfo2.nz_used); 208 if (i < 0 || j < 0 || k1 < 0 || k2 < 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error(compare A and sB): MatGetInfo()\n")); 209 210 PetscCall(MatGetSize(A, &Ii, &J)); 211 PetscCall(MatGetSize(sB, &i, &j)); 212 if (i - Ii || j - J) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetSize()\n")); 213 214 PetscCall(MatGetBlockSize(A, &Ii)); 215 PetscCall(MatGetBlockSize(sB, &i)); 216 if (i - Ii) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetBlockSize()\n")); 217 218 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm)); 219 PetscCall(PetscRandomSetFromOptions(rdm)); 220 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &x)); 221 PetscCall(VecDuplicate(x, &s1)); 222 PetscCall(VecDuplicate(x, &s2)); 223 PetscCall(VecDuplicate(x, &y)); 224 PetscCall(VecDuplicate(x, &b)); 225 PetscCall(VecSetRandom(x, rdm)); 226 227 /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */ 228 #if !defined(PETSC_USE_COMPLEX) 229 /* Scaling matrix with complex numbers results non-spd matrix, 230 causing crash of MatForwardSolve() and MatBackwardSolve() */ 231 PetscCall(MatDiagonalScale(A, x, x)); 232 PetscCall(MatDiagonalScale(sB, x, x)); 233 PetscCall(MatMultEqual(A, sB, 10, &equal)); 234 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "Error in MatDiagonalScale"); 235 236 PetscCall(MatGetDiagonal(A, s1)); 237 PetscCall(MatGetDiagonal(sB, s2)); 238 PetscCall(VecAXPY(s2, neg_one, s1)); 239 PetscCall(VecNorm(s2, NORM_1, &norm1)); 240 if (norm1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatGetDiagonal(), ||s1-s2||=%g\n", (double)norm1)); 241 242 { 243 PetscScalar alpha = 0.1; 244 PetscCall(MatScale(A, alpha)); 245 PetscCall(MatScale(sB, alpha)); 246 } 247 #endif 248 249 /* Test MatGetRowMaxAbs() */ 250 PetscCall(MatGetRowMaxAbs(A, s1, NULL)); 251 PetscCall(MatGetRowMaxAbs(sB, s2, NULL)); 252 PetscCall(VecNorm(s1, NORM_1, &norm1)); 253 PetscCall(VecNorm(s2, NORM_1, &norm2)); 254 norm1 -= norm2; 255 if (norm1 < -tol || norm1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatGetRowMaxAbs() \n")); 256 257 /* Test MatMult() */ 258 for (i = 0; i < 40; i++) { 259 PetscCall(VecSetRandom(x, rdm)); 260 PetscCall(MatMult(A, x, s1)); 261 PetscCall(MatMult(sB, x, s2)); 262 PetscCall(VecNorm(s1, NORM_1, &norm1)); 263 PetscCall(VecNorm(s2, NORM_1, &norm2)); 264 norm1 -= norm2; 265 if (norm1 < -tol || norm1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatMult(), norm1-norm2: %g\n", (double)norm1)); 266 } 267 268 /* MatMultAdd() */ 269 for (i = 0; i < 40; i++) { 270 PetscCall(VecSetRandom(x, rdm)); 271 PetscCall(VecSetRandom(y, rdm)); 272 PetscCall(MatMultAdd(A, x, y, s1)); 273 PetscCall(MatMultAdd(sB, x, y, s2)); 274 PetscCall(VecNorm(s1, NORM_1, &norm1)); 275 PetscCall(VecNorm(s2, NORM_1, &norm2)); 276 norm1 -= norm2; 277 if (norm1 < -tol || norm1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatMultAdd(), norm1-norm2: %g\n", (double)norm1)); 278 } 279 280 /* Test MatMatMult() for sbaij and dense matrices */ 281 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, n, 5 * n, NULL, &B)); 282 PetscCall(MatSetRandom(B, rdm)); 283 PetscCall(MatMatMult(sA, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C)); 284 PetscCall(MatMatMultEqual(sA, B, C, 5 * n, &equal)); 285 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error: MatMatMult()"); 286 PetscCall(MatDestroy(&C)); 287 PetscCall(MatDestroy(&B)); 288 289 /* Test MatCholeskyFactor(), MatICCFactor() with natural ordering */ 290 PetscCall(MatGetOrdering(A, MATORDERINGNATURAL, &perm, &iscol)); 291 PetscCall(ISDestroy(&iscol)); 292 norm1 = tol; 293 inc = bs; 294 295 /* initialize factinfo */ 296 PetscCall(PetscMemzero(&factinfo, sizeof(MatFactorInfo))); 297 298 for (lf = -1; lf < 10; lf += inc) { 299 if (lf == -1) { /* Cholesky factor of sB (duplicate sA) */ 300 factinfo.fill = 5.0; 301 302 PetscCall(MatGetFactor(sB, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sFactor)); 303 PetscCall(MatCholeskyFactorSymbolic(sFactor, sB, perm, &factinfo)); 304 } else if (!doIcc) break; 305 else { /* incomplete Cholesky factor */ factinfo.fill = 5.0; 306 factinfo.levels = lf; 307 308 PetscCall(MatGetFactor(sB, MATSOLVERPETSC, MAT_FACTOR_ICC, &sFactor)); 309 PetscCall(MatICCFactorSymbolic(sFactor, sB, perm, &factinfo)); 310 } 311 PetscCall(MatCholeskyFactorNumeric(sFactor, sB, &factinfo)); 312 /* MatView(sFactor, PETSC_VIEWER_DRAW_WORLD); */ 313 314 /* test MatGetDiagonal on numeric factor */ 315 /* 316 if (lf == -1) { 317 PetscCall(MatGetDiagonal(sFactor,s1)); 318 printf(" in ex74.c, diag: \n"); 319 PetscCall(VecView(s1,PETSC_VIEWER_STDOUT_SELF)); 320 } 321 */ 322 323 PetscCall(MatMult(sB, x, b)); 324 325 /* test MatForwardSolve() and MatBackwardSolve() */ 326 if (lf == -1) { 327 PetscCall(MatForwardSolve(sFactor, b, s1)); 328 PetscCall(MatBackwardSolve(sFactor, s1, s2)); 329 PetscCall(VecAXPY(s2, neg_one, x)); 330 PetscCall(VecNorm(s2, NORM_2, &norm2)); 331 if (10 * norm1 < norm2) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatForwardSolve and BackwardSolve: Norm of error=%g, bs=%" PetscInt_FMT "\n", (double)norm2, bs)); 332 } 333 334 /* test MatSolve() */ 335 PetscCall(MatSolve(sFactor, b, y)); 336 PetscCall(MatDestroy(&sFactor)); 337 /* Check the error */ 338 PetscCall(VecAXPY(y, neg_one, x)); 339 PetscCall(VecNorm(y, NORM_2, &norm2)); 340 if (10 * norm1 < norm2 && lf - inc != -1) PetscCall(PetscPrintf(PETSC_COMM_SELF, "lf=%" PetscInt_FMT ", %" PetscInt_FMT ", Norm of error=%g, %g\n", lf - inc, lf, (double)norm1, (double)norm2)); 341 norm1 = norm2; 342 if (norm2 < tol && lf != -1) break; 343 } 344 345 #if defined(PETSC_HAVE_MUMPS) 346 PetscCall(MatGetFactor(sA, MATSOLVERMUMPS, MAT_FACTOR_CHOLESKY, &sFactor)); 347 PetscCall(MatCholeskyFactorSymbolic(sFactor, sA, NULL, NULL)); 348 PetscCall(MatCholeskyFactorNumeric(sFactor, sA, NULL)); 349 for (i = 0; i < 10; i++) { 350 PetscCall(VecSetRandom(b, rdm)); 351 PetscCall(MatSolve(sFactor, b, y)); 352 /* Check the error */ 353 PetscCall(MatMult(sA, y, x)); 354 PetscCall(VecAXPY(x, neg_one, b)); 355 PetscCall(VecNorm(x, NORM_2, &norm2)); 356 if (norm2 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatSolve(), norm2: %g\n", (double)norm2)); 357 } 358 PetscCall(MatDestroy(&sFactor)); 359 #endif 360 361 PetscCall(ISDestroy(&perm)); 362 363 PetscCall(MatDestroy(&A)); 364 PetscCall(MatDestroy(&sB)); 365 PetscCall(MatDestroy(&sA)); 366 PetscCall(VecDestroy(&x)); 367 PetscCall(VecDestroy(&y)); 368 PetscCall(VecDestroy(&s1)); 369 PetscCall(VecDestroy(&s2)); 370 PetscCall(VecDestroy(&b)); 371 PetscCall(PetscRandomDestroy(&rdm)); 372 373 PetscCall(PetscFinalize()); 374 return 0; 375 } 376 377 /*TEST 378 379 test: 380 args: -bs {{1 2 3 4 5 6 7 8}} 381 382 TEST*/ 383