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