1 static char help[] = "Illustrate how to do one symbolic factorization and multiple numeric factorizations using same matrix nonzero structure. \n\n"; 2 3 #include <petscmat.h> 4 5 int main(int argc, char **args) 6 { 7 PetscInt i, rstart, rend, N = 10, num_numfac = 5, col[3], k; 8 Mat A[5], F; 9 Vec u, x, b; 10 PetscMPIInt rank; 11 PetscScalar value[3]; 12 PetscReal norm, tol = 100 * PETSC_MACHINE_EPSILON; 13 IS perm, iperm; 14 MatFactorInfo info; 15 MatFactorType facttype = MAT_FACTOR_LU; 16 char solvertype[64]; 17 char factortype[64]; 18 19 PetscFunctionBeginUser; 20 PetscCall(PetscInitialize(&argc, &args, NULL, help)); 21 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 22 PetscCall(PetscOptionsGetReal(NULL, NULL, "-tol", &tol, NULL)); 23 24 /* Create and assemble matrices, all have same data structure */ 25 for (k = 0; k < num_numfac; k++) { 26 PetscCall(MatCreate(PETSC_COMM_WORLD, &A[k])); 27 PetscCall(MatSetSizes(A[k], PETSC_DECIDE, PETSC_DECIDE, N, N)); 28 PetscCall(MatSetFromOptions(A[k])); 29 PetscCall(MatSetUp(A[k])); 30 PetscCall(MatGetOwnershipRange(A[k], &rstart, &rend)); 31 32 value[0] = -1.0 * (k + 1); 33 value[1] = 2.0 * (k + 1); 34 value[2] = -1.0 * (k + 1); 35 for (i = rstart; i < rend; i++) { 36 col[0] = i - 1; 37 col[1] = i; 38 col[2] = i + 1; 39 if (i == 0) { 40 PetscCall(MatSetValues(A[k], 1, &i, 2, col + 1, value + 1, INSERT_VALUES)); 41 } else if (i == N - 1) { 42 PetscCall(MatSetValues(A[k], 1, &i, 2, col, value, INSERT_VALUES)); 43 } else { 44 PetscCall(MatSetValues(A[k], 1, &i, 3, col, value, INSERT_VALUES)); 45 } 46 } 47 PetscCall(MatAssemblyBegin(A[k], MAT_FINAL_ASSEMBLY)); 48 PetscCall(MatAssemblyEnd(A[k], MAT_FINAL_ASSEMBLY)); 49 PetscCall(MatSetOption(A[k], MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 50 } 51 52 /* Create vectors */ 53 PetscCall(MatCreateVecs(A[0], &x, &b)); 54 PetscCall(VecDuplicate(x, &u)); 55 56 /* Set rhs vector b */ 57 PetscCall(VecSet(b, 1.0)); 58 59 /* Get a symbolic factor F from A[0] */ 60 PetscCall(PetscStrncpy(solvertype, "petsc", sizeof(solvertype))); 61 PetscCall(PetscOptionsGetString(NULL, NULL, "-mat_solver_type", solvertype, sizeof(solvertype), NULL)); 62 PetscCall(PetscOptionsGetEnum(NULL, NULL, "-mat_factor_type", MatFactorTypes, (PetscEnum *)&facttype, NULL)); 63 64 PetscCall(MatGetFactor(A[0], solvertype, facttype, &F)); 65 /* test mumps options */ 66 PetscCall(MatMumpsSetIcntl(F, 7, 5)); 67 PetscCall(PetscStrncpy(factortype, MatFactorTypes[facttype], sizeof(factortype))); 68 PetscCall(PetscStrtoupper(solvertype)); 69 PetscCall(PetscStrtoupper(factortype)); 70 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %s %s:\n", solvertype, factortype)); 71 72 PetscCall(MatFactorInfoInitialize(&info)); 73 info.fill = 5.0; 74 PetscCall(MatGetOrdering(A[0], MATORDERINGNATURAL, &perm, &iperm)); 75 switch (facttype) { 76 case MAT_FACTOR_LU: 77 PetscCall(MatLUFactorSymbolic(F, A[0], perm, iperm, &info)); 78 break; 79 case MAT_FACTOR_ILU: 80 PetscCall(MatILUFactorSymbolic(F, A[0], perm, iperm, &info)); 81 break; 82 case MAT_FACTOR_ICC: 83 PetscCall(MatICCFactorSymbolic(F, A[0], perm, &info)); 84 break; 85 case MAT_FACTOR_CHOLESKY: 86 PetscCall(MatCholeskyFactorSymbolic(F, A[0], perm, &info)); 87 break; 88 default: 89 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not for factor type %s", factortype); 90 } 91 92 /* Compute numeric factors using same F, then solve */ 93 for (k = 0; k < num_numfac; k++) { 94 switch (facttype) { 95 case MAT_FACTOR_LU: 96 case MAT_FACTOR_ILU: 97 PetscCall(MatLUFactorNumeric(F, A[k], &info)); 98 break; 99 case MAT_FACTOR_ICC: 100 case MAT_FACTOR_CHOLESKY: 101 PetscCall(MatCholeskyFactorNumeric(F, A[k], &info)); 102 break; 103 default: 104 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not for factor type %s", factortype); 105 } 106 107 /* Solve A[k] * x = b */ 108 PetscCall(MatSolve(F, b, x)); 109 110 /* Check the residual */ 111 PetscCall(MatMult(A[k], x, u)); 112 PetscCall(VecAXPY(u, -1.0, b)); 113 PetscCall(VecNorm(u, NORM_INFINITY, &norm)); 114 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the %s numfact and solve: residual %g\n", k, factortype, (double)norm)); 115 } 116 117 /* Free data structures */ 118 for (k = 0; k < num_numfac; k++) PetscCall(MatDestroy(&A[k])); 119 PetscCall(MatDestroy(&F)); 120 PetscCall(ISDestroy(&perm)); 121 PetscCall(ISDestroy(&iperm)); 122 PetscCall(VecDestroy(&x)); 123 PetscCall(VecDestroy(&b)); 124 PetscCall(VecDestroy(&u)); 125 PetscCall(PetscFinalize()); 126 return 0; 127 } 128 129 /*TEST 130 131 test: 132 133 test: 134 suffix: 2 135 args: -mat_solver_type superlu 136 requires: superlu 137 138 testset: 139 nsize: 2 140 requires: mumps 141 args: -mat_solver_type mumps 142 output_file: output/ex28_3.out 143 144 test: 145 suffix: 3 146 requires: !__float128 147 test: 148 suffix: 3_fp128 149 requires: __float128 150 args: -tol 1e-14 151 152 test: 153 suffix: 4 154 args: -mat_solver_type cusparse -mat_type aijcusparse -mat_factor_type {{lu cholesky ilu icc}separate output} 155 requires: cuda 156 157 TEST*/ 158