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