1 2 static char help[] = "Tests cholesky, icc factorization and solve on sequential aij, baij and sbaij matrices. \n"; 3 4 #include <petscmat.h> 5 6 int main(int argc, char **args) 7 { 8 Vec x, y, b; 9 Mat A; /* linear system matrix */ 10 Mat sA, sC; /* symmetric part of the matrices */ 11 PetscInt n, mbs = 16, bs = 1, nz = 3, prob = 1, i, j, col[3], block, row, Ii, J, n1, lvl; 12 PetscMPIInt size; 13 PetscReal norm2; 14 PetscScalar neg_one = -1.0, four = 4.0, value[3]; 15 IS perm, cperm; 16 PetscRandom rdm; 17 PetscBool reorder = PETSC_FALSE, displ = PETSC_FALSE; 18 MatFactorInfo factinfo; 19 PetscBool equal; 20 PetscBool TestAIJ = PETSC_FALSE, TestBAIJ = PETSC_TRUE; 21 PetscInt TestShift = 0; 22 23 PetscFunctionBeginUser; 24 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 25 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 26 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); 27 PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL)); 28 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mbs", &mbs, NULL)); 29 PetscCall(PetscOptionsGetBool(NULL, NULL, "-reorder", &reorder, NULL)); 30 PetscCall(PetscOptionsGetBool(NULL, NULL, "-testaij", &TestAIJ, NULL)); 31 PetscCall(PetscOptionsGetInt(NULL, NULL, "-testShift", &TestShift, NULL)); 32 PetscCall(PetscOptionsGetBool(NULL, NULL, "-displ", &displ, NULL)); 33 34 n = mbs * bs; 35 if (TestAIJ) { /* A is in aij format */ 36 PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, n, n, nz, NULL, &A)); 37 TestBAIJ = PETSC_FALSE; 38 } else { /* A is in baij format */ 39 PetscCall(MatCreateSeqBAIJ(PETSC_COMM_WORLD, bs, n, n, nz, NULL, &A)); 40 TestAIJ = PETSC_FALSE; 41 } 42 43 /* Assemble matrix */ 44 if (bs == 1) { 45 PetscCall(PetscOptionsGetInt(NULL, NULL, "-test_problem", &prob, NULL)); 46 if (prob == 1) { /* tridiagonal matrix */ 47 value[0] = -1.0; 48 value[1] = 2.0; 49 value[2] = -1.0; 50 for (i = 1; i < n - 1; i++) { 51 col[0] = i - 1; 52 col[1] = i; 53 col[2] = i + 1; 54 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES)); 55 } 56 i = n - 1; 57 col[0] = 0; 58 col[1] = n - 2; 59 col[2] = n - 1; 60 61 value[0] = 0.1; 62 value[1] = -1; 63 value[2] = 2; 64 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES)); 65 66 i = 0; 67 col[0] = 0; 68 col[1] = 1; 69 col[2] = n - 1; 70 71 value[0] = 2.0; 72 value[1] = -1.0; 73 value[2] = 0.1; 74 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES)); 75 } else if (prob == 2) { /* matrix for the five point stencil */ 76 n1 = (PetscInt)(PetscSqrtReal((PetscReal)n) + 0.001); 77 PetscCheck(n1 * n1 == n, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "sqrt(n) must be a positive integer!"); 78 for (i = 0; i < n1; i++) { 79 for (j = 0; j < n1; j++) { 80 Ii = j + n1 * i; 81 if (i > 0) { 82 J = Ii - n1; 83 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 84 } 85 if (i < n1 - 1) { 86 J = Ii + n1; 87 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 88 } 89 if (j > 0) { 90 J = Ii - 1; 91 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 92 } 93 if (j < n1 - 1) { 94 J = Ii + 1; 95 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES)); 96 } 97 PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &four, INSERT_VALUES)); 98 } 99 } 100 } 101 } else { /* bs > 1 */ 102 for (block = 0; block < n / bs; block++) { 103 /* diagonal blocks */ 104 value[0] = -1.0; 105 value[1] = 4.0; 106 value[2] = -1.0; 107 for (i = 1 + block * bs; i < bs - 1 + block * bs; i++) { 108 col[0] = i - 1; 109 col[1] = i; 110 col[2] = i + 1; 111 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES)); 112 } 113 i = bs - 1 + block * bs; 114 col[0] = bs - 2 + block * bs; 115 col[1] = bs - 1 + block * bs; 116 117 value[0] = -1.0; 118 value[1] = 4.0; 119 PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES)); 120 121 i = 0 + block * bs; 122 col[0] = 0 + block * bs; 123 col[1] = 1 + block * bs; 124 125 value[0] = 4.0; 126 value[1] = -1.0; 127 PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES)); 128 } 129 /* off-diagonal blocks */ 130 value[0] = -1.0; 131 for (i = 0; i < (n / bs - 1) * bs; i++) { 132 col[0] = i + bs; 133 PetscCall(MatSetValues(A, 1, &i, 1, col, value, INSERT_VALUES)); 134 col[0] = i; 135 row = i + bs; 136 PetscCall(MatSetValues(A, 1, &row, 1, col, value, INSERT_VALUES)); 137 } 138 } 139 140 if (TestShift) { 141 /* set diagonals in the 0-th block as 0 for testing shift numerical factor */ 142 for (i = 0; i < bs; i++) { 143 row = i; 144 col[0] = i; 145 value[0] = 0.0; 146 PetscCall(MatSetValues(A, 1, &row, 1, col, value, INSERT_VALUES)); 147 } 148 } 149 150 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 151 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 152 153 /* Test MatConvert */ 154 PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE)); 155 PetscCall(MatConvert(A, MATSEQSBAIJ, MAT_INITIAL_MATRIX, &sA)); 156 PetscCall(MatMultEqual(A, sA, 20, &equal)); 157 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_USER, "A != sA"); 158 159 /* Test MatGetOwnershipRange() */ 160 PetscCall(MatGetOwnershipRange(A, &Ii, &J)); 161 PetscCall(MatGetOwnershipRange(sA, &i, &j)); 162 PetscCheck(i == Ii && j == J, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatGetOwnershipRange() in MatSBAIJ format"); 163 164 /* Vectors */ 165 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm)); 166 PetscCall(PetscRandomSetFromOptions(rdm)); 167 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &x)); 168 PetscCall(VecDuplicate(x, &b)); 169 PetscCall(VecDuplicate(x, &y)); 170 PetscCall(VecSetRandom(x, rdm)); 171 172 /* Test MatReordering() - not work on sbaij matrix */ 173 if (reorder) { 174 PetscCall(MatGetOrdering(A, MATORDERINGRCM, &perm, &cperm)); 175 } else { 176 PetscCall(MatGetOrdering(A, MATORDERINGNATURAL, &perm, &cperm)); 177 } 178 PetscCall(ISDestroy(&cperm)); 179 180 /* initialize factinfo */ 181 PetscCall(MatFactorInfoInitialize(&factinfo)); 182 if (TestShift == 1) { 183 factinfo.shifttype = (PetscReal)MAT_SHIFT_NONZERO; 184 factinfo.shiftamount = 0.1; 185 } else if (TestShift == 2) { 186 factinfo.shifttype = (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE; 187 } 188 189 /* Test MatCholeskyFactor(), MatICCFactor() */ 190 /*------------------------------------------*/ 191 /* Test aij matrix A */ 192 if (TestAIJ) { 193 if (displ) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "AIJ: \n")); 194 i = 0; 195 for (lvl = -1; lvl < 10; lvl++) { 196 if (lvl == -1) { /* Cholesky factor */ 197 factinfo.fill = 5.0; 198 199 PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sC)); 200 PetscCall(MatCholeskyFactorSymbolic(sC, A, perm, &factinfo)); 201 } else { /* incomplete Cholesky factor */ 202 factinfo.fill = 5.0; 203 factinfo.levels = lvl; 204 205 PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_ICC, &sC)); 206 PetscCall(MatICCFactorSymbolic(sC, A, perm, &factinfo)); 207 } 208 PetscCall(MatCholeskyFactorNumeric(sC, A, &factinfo)); 209 210 PetscCall(MatMult(A, x, b)); 211 PetscCall(MatSolve(sC, b, y)); 212 PetscCall(MatDestroy(&sC)); 213 214 /* Check the residual */ 215 PetscCall(VecAXPY(y, neg_one, x)); 216 PetscCall(VecNorm(y, NORM_2, &norm2)); 217 218 if (displ) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " lvl: %" PetscInt_FMT ", residual: %g\n", lvl, (double)norm2)); 219 } 220 } 221 222 /* Test baij matrix A */ 223 if (TestBAIJ) { 224 if (displ) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "BAIJ: \n")); 225 i = 0; 226 for (lvl = -1; lvl < 10; lvl++) { 227 if (lvl == -1) { /* Cholesky factor */ 228 factinfo.fill = 5.0; 229 230 PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sC)); 231 PetscCall(MatCholeskyFactorSymbolic(sC, A, perm, &factinfo)); 232 } else { /* incomplete Cholesky factor */ 233 factinfo.fill = 5.0; 234 factinfo.levels = lvl; 235 236 PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_ICC, &sC)); 237 PetscCall(MatICCFactorSymbolic(sC, A, perm, &factinfo)); 238 } 239 PetscCall(MatCholeskyFactorNumeric(sC, A, &factinfo)); 240 241 PetscCall(MatMult(A, x, b)); 242 PetscCall(MatSolve(sC, b, y)); 243 PetscCall(MatDestroy(&sC)); 244 245 /* Check the residual */ 246 PetscCall(VecAXPY(y, neg_one, x)); 247 PetscCall(VecNorm(y, NORM_2, &norm2)); 248 if (displ) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " lvl: %" PetscInt_FMT ", residual: %g\n", lvl, (double)norm2)); 249 } 250 } 251 252 /* Test sbaij matrix sA */ 253 if (displ) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "SBAIJ: \n")); 254 i = 0; 255 for (lvl = -1; lvl < 10; lvl++) { 256 if (lvl == -1) { /* Cholesky factor */ 257 factinfo.fill = 5.0; 258 259 PetscCall(MatGetFactor(sA, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sC)); 260 PetscCall(MatCholeskyFactorSymbolic(sC, sA, perm, &factinfo)); 261 } else { /* incomplete Cholesky factor */ 262 factinfo.fill = 5.0; 263 factinfo.levels = lvl; 264 265 PetscCall(MatGetFactor(sA, MATSOLVERPETSC, MAT_FACTOR_ICC, &sC)); 266 PetscCall(MatICCFactorSymbolic(sC, sA, perm, &factinfo)); 267 } 268 PetscCall(MatCholeskyFactorNumeric(sC, sA, &factinfo)); 269 270 if (lvl == 0 && bs == 1) { /* Test inplace ICC(0) for sbaij sA - does not work for new datastructure */ 271 /* 272 Mat B; 273 PetscCall(MatDuplicate(sA,MAT_COPY_VALUES,&B)); 274 PetscCall(MatICCFactor(B,perm,&factinfo)); 275 PetscCall(MatEqual(sC,B,&equal)); 276 if (!equal) { 277 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"in-place Cholesky factor != out-place Cholesky factor"); 278 } 279 PetscCall(MatDestroy(&B)); 280 */ 281 } 282 283 PetscCall(MatMult(sA, x, b)); 284 PetscCall(MatSolve(sC, b, y)); 285 286 /* Test MatSolves() */ 287 if (bs == 1) { 288 Vecs xx, bb; 289 PetscCall(VecsCreateSeq(PETSC_COMM_SELF, n, 4, &xx)); 290 PetscCall(VecsDuplicate(xx, &bb)); 291 PetscCall(MatSolves(sC, bb, xx)); 292 PetscCall(VecsDestroy(xx)); 293 PetscCall(VecsDestroy(bb)); 294 } 295 PetscCall(MatDestroy(&sC)); 296 297 /* Check the residual */ 298 PetscCall(VecAXPY(y, neg_one, x)); 299 PetscCall(VecNorm(y, NORM_2, &norm2)); 300 if (displ) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " lvl: %" PetscInt_FMT ", residual: %g\n", lvl, (double)norm2)); 301 } 302 303 PetscCall(ISDestroy(&perm)); 304 PetscCall(MatDestroy(&A)); 305 PetscCall(MatDestroy(&sA)); 306 PetscCall(VecDestroy(&x)); 307 PetscCall(VecDestroy(&y)); 308 PetscCall(VecDestroy(&b)); 309 PetscCall(PetscRandomDestroy(&rdm)); 310 311 PetscCall(PetscFinalize()); 312 return 0; 313 } 314 315 /*TEST 316 317 test: 318 args: -bs {{1 2 3 4 5 6 7 8}} 319 320 test: 321 suffix: 3 322 args: -testaij 323 output_file: output/ex76_1.out 324 325 TEST*/ 326