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