Actual source code: ex28.c

  1: static char help[] = "Illustrate how to do one symbolic factorization and multiple numeric factorizations using same matrix 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, (char *)0, 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: #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));

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

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

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

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

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

130: /*TEST

132:    test:

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

139:    test:
140:       suffix: 3
141:       nsize: 2
142:       requires: mumps
143:       args: -mat_solver_type mumps

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

150: TEST*/