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 { 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 #if defined(PETSC_HAVE_MUMPS) 66 PetscCall(MatMumpsSetIcntl(F, 7, 5)); 67 #endif 68 PetscCall(PetscStrncpy(factortype, MatFactorTypes[facttype], sizeof(factortype))); 69 PetscCall(PetscStrtoupper(solvertype)); 70 PetscCall(PetscStrtoupper(factortype)); 71 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %s %s:\n", solvertype, factortype)); 72 73 PetscCall(MatFactorInfoInitialize(&info)); 74 info.fill = 5.0; 75 PetscCall(MatGetOrdering(A[0], MATORDERINGNATURAL, &perm, &iperm)); 76 switch (facttype) { 77 case MAT_FACTOR_LU: 78 PetscCall(MatLUFactorSymbolic(F, A[0], perm, iperm, &info)); 79 break; 80 case MAT_FACTOR_ILU: 81 PetscCall(MatILUFactorSymbolic(F, A[0], perm, iperm, &info)); 82 break; 83 case MAT_FACTOR_ICC: 84 PetscCall(MatICCFactorSymbolic(F, A[0], perm, &info)); 85 break; 86 case MAT_FACTOR_CHOLESKY: 87 PetscCall(MatCholeskyFactorSymbolic(F, A[0], perm, &info)); 88 break; 89 default: 90 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not for factor type %s", factortype); 91 } 92 93 /* Compute numeric factors using same F, then solve */ 94 for (k = 0; k < num_numfac; k++) { 95 switch (facttype) { 96 case MAT_FACTOR_LU: 97 case MAT_FACTOR_ILU: 98 PetscCall(MatLUFactorNumeric(F, A[k], &info)); 99 break; 100 case MAT_FACTOR_ICC: 101 case MAT_FACTOR_CHOLESKY: 102 PetscCall(MatCholeskyFactorNumeric(F, A[k], &info)); 103 break; 104 default: 105 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not for factor type %s", factortype); 106 } 107 108 /* Solve A[k] * x = b */ 109 PetscCall(MatSolve(F, b, x)); 110 111 /* Check the residual */ 112 PetscCall(MatMult(A[k], x, u)); 113 PetscCall(VecAXPY(u, -1.0, b)); 114 PetscCall(VecNorm(u, NORM_INFINITY, &norm)); 115 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the %s numfact and solve: residual %g\n", k, factortype, (double)norm)); 116 } 117 118 /* Free data structures */ 119 for (k = 0; k < num_numfac; k++) PetscCall(MatDestroy(&A[k])); 120 PetscCall(MatDestroy(&F)); 121 PetscCall(ISDestroy(&perm)); 122 PetscCall(ISDestroy(&iperm)); 123 PetscCall(VecDestroy(&x)); 124 PetscCall(VecDestroy(&b)); 125 PetscCall(VecDestroy(&u)); 126 PetscCall(PetscFinalize()); 127 return 0; 128 } 129 130 /*TEST 131 132 test: 133 134 test: 135 suffix: 2 136 args: -mat_solver_type superlu 137 requires: superlu 138 139 test: 140 suffix: 3 141 nsize: 2 142 requires: mumps 143 args: -mat_solver_type mumps 144 145 test: 146 suffix: 4 147 args: -mat_solver_type cusparse -mat_type aijcusparse -mat_factor_type {{lu cholesky ilu icc}separate output} 148 requires: cuda 149 150 TEST*/ 151