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

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

 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:   }

 51:   /* Create vectors */
 52:   PetscCall(MatCreateVecs(A[0], &x, &b));
 53:   PetscCall(VecDuplicate(x, &u));

 55:   /* Set rhs vector b */
 56:   PetscCall(VecSet(b, 1.0));

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

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

 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:   }

 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:     }

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

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:   }

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: }

128: /*TEST

130:    test:

132:    test:
133:       suffix: 2
134:       args: -mat_solver_type superlu
135:       requires: superlu

137:    test:
138:       suffix: 3
139:       nsize: 2
140:       requires: mumps
141:       args: -mat_solver_type mumps

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

148: TEST*/