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