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