1 static char help[] = "Tests ILU and ICC factorization with and without matrix ordering on seqsbaij format. Modified from ex30.c\n\ 2 Input parameters are:\n\ 3 -lf <level> : level of fill for ILU (default is 0)\n\ 4 -lu : use full LU or Cholesky factorization\n\ 5 -m <value>,-n <value> : grid dimensions\n\ 6 Note that most users should employ the KSP interface to the\n\ 7 linear solvers instead of using the factorization routines\n\ 8 directly.\n\n"; 9 10 #include <petscmat.h> 11 12 int main(int argc, char **args) 13 { 14 Mat C, sC, sA; 15 PetscInt i, j, m = 5, n = 5, Ii, J, lf = 0; 16 PetscBool CHOLESKY = PETSC_FALSE, TRIANGULAR = PETSC_FALSE, flg; 17 PetscScalar v; 18 IS row, col; 19 MatFactorInfo info; 20 Vec x, y, b, ytmp; 21 PetscReal norm2, tol = 100 * PETSC_MACHINE_EPSILON; 22 PetscRandom rdm; 23 PetscMPIInt size; 24 25 PetscFunctionBeginUser; 26 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 27 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 28 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); 29 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL)); 30 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 31 PetscCall(PetscOptionsGetInt(NULL, NULL, "-lf", &lf, NULL)); 32 33 PetscCall(MatCreate(PETSC_COMM_SELF, &C)); 34 PetscCall(MatSetSizes(C, m * n, m * n, m * n, m * n)); 35 PetscCall(MatSetFromOptions(C)); 36 PetscCall(MatSetUp(C)); 37 38 /* Create matrix C in seqaij format and sC in seqsbaij. (This is five-point stencil with some extra elements) */ 39 for (i = 0; i < m; i++) { 40 for (j = 0; j < n; j++) { 41 v = -1.0; 42 Ii = j + n * i; 43 J = Ii - n; 44 if (J >= 0) PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 45 J = Ii + n; 46 if (J < m * n) PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 47 J = Ii - 1; 48 if (J >= 0) PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 49 J = Ii + 1; 50 if (J < m * n) PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES)); 51 v = 4.0; 52 PetscCall(MatSetValues(C, 1, &Ii, 1, &Ii, &v, INSERT_VALUES)); 53 } 54 } 55 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); 56 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); 57 58 PetscCall(MatIsSymmetric(C, 0.0, &flg)); 59 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "C is non-symmetric"); 60 PetscCall(MatConvert(C, MATSEQSBAIJ, MAT_INITIAL_MATRIX, &sC)); 61 62 /* Create vectors for error checking */ 63 PetscCall(MatCreateVecs(C, &x, &b)); 64 PetscCall(VecDuplicate(x, &y)); 65 PetscCall(VecDuplicate(x, &ytmp)); 66 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm)); 67 PetscCall(PetscRandomSetFromOptions(rdm)); 68 PetscCall(VecSetRandom(x, rdm)); 69 PetscCall(MatMult(C, x, b)); 70 71 PetscCall(MatGetOrdering(C, MATORDERINGNATURAL, &row, &col)); 72 73 /* Compute CHOLESKY or ICC factor sA */ 74 PetscCall(MatFactorInfoInitialize(&info)); 75 76 info.fill = 1.0; 77 info.diagonal_fill = 0; 78 info.zeropivot = 0.0; 79 80 PetscCall(PetscOptionsHasName(NULL, NULL, "-cholesky", &CHOLESKY)); 81 if (CHOLESKY) { 82 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Test CHOLESKY...\n")); 83 PetscCall(MatGetFactor(sC, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sA)); 84 PetscCall(MatCholeskyFactorSymbolic(sA, sC, row, &info)); 85 } else { 86 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Test ICC...\n")); 87 info.levels = lf; 88 89 PetscCall(MatGetFactor(sC, MATSOLVERPETSC, MAT_FACTOR_ICC, &sA)); 90 PetscCall(MatICCFactorSymbolic(sA, sC, row, &info)); 91 } 92 PetscCall(MatCholeskyFactorNumeric(sA, sC, &info)); 93 94 /* test MatForwardSolve() and MatBackwardSolve() with matrix reordering on aij matrix C */ 95 if (CHOLESKY) { 96 PetscCall(PetscOptionsHasName(NULL, NULL, "-triangular_solve", &TRIANGULAR)); 97 if (TRIANGULAR) { 98 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Test MatForwardSolve...\n")); 99 PetscCall(MatForwardSolve(sA, b, ytmp)); 100 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Test MatBackwardSolve...\n")); 101 PetscCall(MatBackwardSolve(sA, ytmp, y)); 102 PetscCall(VecAXPY(y, -1.0, x)); 103 PetscCall(VecNorm(y, NORM_2, &norm2)); 104 if (norm2 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatForwardSolve and BackwardSolve: Norm of error=%g\n", (double)norm2)); 105 } 106 } 107 108 PetscCall(MatSolve(sA, b, y)); 109 PetscCall(MatDestroy(&sC)); 110 PetscCall(MatDestroy(&sA)); 111 PetscCall(VecAXPY(y, -1.0, x)); 112 PetscCall(VecNorm(y, NORM_2, &norm2)); 113 if (lf == -1 && norm2 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, " reordered SEQAIJ: Cholesky/ICC levels %" PetscInt_FMT ", residual %g\n", lf, (double)norm2)); 114 115 /* Free data structures */ 116 PetscCall(MatDestroy(&C)); 117 PetscCall(ISDestroy(&row)); 118 PetscCall(ISDestroy(&col)); 119 PetscCall(PetscRandomDestroy(&rdm)); 120 PetscCall(VecDestroy(&x)); 121 PetscCall(VecDestroy(&y)); 122 PetscCall(VecDestroy(&ytmp)); 123 PetscCall(VecDestroy(&b)); 124 PetscCall(PetscFinalize()); 125 return 0; 126 } 127 128 /*TEST 129 130 test: 131 output_file: output/ex128.out 132 133 test: 134 suffix: 2 135 args: -cholesky -triangular_solve 136 137 TEST*/ 138