1: static char help[] = "Illustrate how to do one symbolic factorization and multiple numeric factorizations using same matrix nonzero structure. \n\n";

  3: #include <petscmat.h>

  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];

 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));

 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));

 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) PetscCall(MatSetValues(A[k], 1, &i, 2, col + 1, value + 1, INSERT_VALUES));
 40:       else if (i == N - 1) PetscCall(MatSetValues(A[k], 1, &i, 2, col, value, INSERT_VALUES));
 41:       else PetscCall(MatSetValues(A[k], 1, &i, 3, col, value, INSERT_VALUES));
 42:     }
 43:     PetscCall(MatAssemblyBegin(A[k], MAT_FINAL_ASSEMBLY));
 44:     PetscCall(MatAssemblyEnd(A[k], MAT_FINAL_ASSEMBLY));
 45:     PetscCall(MatSetOption(A[k], MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
 46:   }

 48:   /* Create vectors */
 49:   PetscCall(MatCreateVecs(A[0], &x, &b));
 50:   PetscCall(VecDuplicate(x, &u));

 52:   /* Set rhs vector b */
 53:   PetscCall(VecSet(b, 1.0));

 55:   /* Get a symbolic factor F from A[0] */
 56:   PetscCall(PetscStrncpy(solvertype, "petsc", sizeof(solvertype)));
 57:   PetscCall(PetscOptionsGetString(NULL, NULL, "-mat_solver_type", solvertype, sizeof(solvertype), NULL));
 58:   PetscCall(PetscOptionsGetEnum(NULL, NULL, "-mat_factor_type", MatFactorTypes, (PetscEnum *)&facttype, NULL));

 60:   PetscCall(MatGetFactor(A[0], solvertype, facttype, &F));
 61:   /* test mumps options */
 62:   PetscCall(MatMumpsSetIcntl(F, 7, 5));
 63:   PetscCall(PetscStrncpy(factortype, MatFactorTypes[facttype], sizeof(factortype)));
 64:   PetscCall(PetscStrtoupper(solvertype));
 65:   PetscCall(PetscStrtoupper(factortype));
 66:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %s %s:\n", solvertype, factortype));

 68:   PetscCall(MatFactorInfoInitialize(&info));
 69:   info.fill = 5.0;
 70:   PetscCall(MatGetOrdering(A[0], MATORDERINGNATURAL, &perm, &iperm));
 71:   switch (facttype) {
 72:   case MAT_FACTOR_LU:
 73:     PetscCall(MatLUFactorSymbolic(F, A[0], perm, iperm, &info));
 74:     break;
 75:   case MAT_FACTOR_ILU:
 76:     PetscCall(MatILUFactorSymbolic(F, A[0], perm, iperm, &info));
 77:     break;
 78:   case MAT_FACTOR_ICC:
 79:     PetscCall(MatICCFactorSymbolic(F, A[0], perm, &info));
 80:     break;
 81:   case MAT_FACTOR_CHOLESKY:
 82:     PetscCall(MatCholeskyFactorSymbolic(F, A[0], perm, &info));
 83:     break;
 84:   default:
 85:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not for factor type %s", factortype);
 86:   }

 88:   /* Compute numeric factors using same F, then solve */
 89:   for (k = 0; k < num_numfac; k++) {
 90:     switch (facttype) {
 91:     case MAT_FACTOR_LU:
 92:     case MAT_FACTOR_ILU:
 93:       PetscCall(MatLUFactorNumeric(F, A[k], &info));
 94:       break;
 95:     case MAT_FACTOR_ICC:
 96:     case MAT_FACTOR_CHOLESKY:
 97:       PetscCall(MatCholeskyFactorNumeric(F, A[k], &info));
 98:       break;
 99:     default:
100:       SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not for factor type %s", factortype);
101:     }

103:     /* Solve A[k] * x = b */
104:     PetscCall(MatSolve(F, b, x));

106:     /* Check the residual */
107:     PetscCall(MatMult(A[k], x, u));
108:     PetscCall(VecAXPY(u, -1.0, b));
109:     PetscCall(VecNorm(u, NORM_INFINITY, &norm));
110:     if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the %s numfact and solve: residual %g\n", k, factortype, (double)norm));
111:   }

113:   /* Free data structures */
114:   for (k = 0; k < num_numfac; k++) PetscCall(MatDestroy(&A[k]));
115:   PetscCall(MatDestroy(&F));
116:   PetscCall(ISDestroy(&perm));
117:   PetscCall(ISDestroy(&iperm));
118:   PetscCall(VecDestroy(&x));
119:   PetscCall(VecDestroy(&b));
120:   PetscCall(VecDestroy(&u));
121:   PetscCall(PetscFinalize());
122:   return 0;
123: }

125: /*TEST

127:    test:

129:    test:
130:       suffix: 2
131:       args: -mat_solver_type superlu
132:       requires: superlu

134:    testset:
135:       nsize: 2
136:       requires: mumps
137:       args: -mat_solver_type mumps
138:       output_file: output/ex28_3.out

140:       test:
141:         suffix: 3
142:         requires: !__float128
143:       test:
144:         suffix: 3_fp128
145:         requires: __float128
146:         args: -tol 1e-14

148:    test:
149:       suffix: 4
150:       args: -mat_solver_type cusparse -mat_type aijcusparse -mat_factor_type {{lu cholesky ilu icc}separate output}
151:       requires: cuda

153: TEST*/