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

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

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

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

 64:   PetscCall(MatGetFactor(A[0], solvertype, facttype, &F));
 65:   /* test mumps options */
 66:   PetscCall(MatMumpsSetIcntl(F, 7, 5));
 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));

 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:
 77:     PetscCall(MatLUFactorSymbolic(F, A[0], perm, iperm, &info));
 78:     break;
 79:   case MAT_FACTOR_ILU:
 80:     PetscCall(MatILUFactorSymbolic(F, A[0], perm, iperm, &info));
 81:     break;
 82:   case MAT_FACTOR_ICC:
 83:     PetscCall(MatICCFactorSymbolic(F, A[0], perm, &info));
 84:     break;
 85:   case MAT_FACTOR_CHOLESKY:
 86:     PetscCall(MatCholeskyFactorSymbolic(F, A[0], perm, &info));
 87:     break;
 88:   default:
 89:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not for factor type %s", factortype);
 90:   }

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

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

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

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

129: /*TEST

131:    test:

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

138:    testset:
139:       nsize: 2
140:       requires: mumps
141:       args: -mat_solver_type mumps
142:       output_file: output/ex28_3.out

144:       test:
145:         suffix: 3
146:         requires: !__float128
147:       test:
148:         suffix: 3_fp128
149:         requires: __float128
150:         args: -tol 1e-14

152:    test:
153:       suffix: 4
154:       args: -mat_solver_type cusparse -mat_type aijcusparse -mat_factor_type {{lu cholesky ilu icc}separate output}
155:       requires: cuda

157: TEST*/