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

 20:   PetscInitialize(&argc, &args, (char *)0, help);
 21:   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:     MatCreate(PETSC_COMM_WORLD, &A[k]);
 26:     MatSetSizes(A[k], PETSC_DECIDE, PETSC_DECIDE, N, N);
 27:     MatSetFromOptions(A[k]);
 28:     MatSetUp(A[k]);
 29:     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:         MatSetValues(A[k], 1, &i, 2, col + 1, value + 1, INSERT_VALUES);
 40:       } else if (i == N - 1) {
 41:         MatSetValues(A[k], 1, &i, 2, col, value, INSERT_VALUES);
 42:       } else {
 43:         MatSetValues(A[k], 1, &i, 3, col, value, INSERT_VALUES);
 44:       }
 45:     }
 46:     MatAssemblyBegin(A[k], MAT_FINAL_ASSEMBLY);
 47:     MatAssemblyEnd(A[k], MAT_FINAL_ASSEMBLY);
 48:     MatSetOption(A[k], MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
 49:   }

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

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

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

 63:   MatGetFactor(A[0], solvertype, facttype, &F);
 64:   /* test mumps options */
 65: #if defined(PETSC_HAVE_MUMPS)
 66:   MatMumpsSetIcntl(F, 7, 5);
 67: #endif
 68:   PetscStrncpy(factortype, MatFactorTypes[facttype], sizeof(factortype));
 69:   PetscStrtoupper(solvertype);
 70:   PetscStrtoupper(factortype);
 71:   PetscPrintf(PETSC_COMM_WORLD, " %s %s:\n", solvertype, factortype);

 73:   MatFactorInfoInitialize(&info);
 74:   info.fill = 5.0;
 75:   MatGetOrdering(A[0], MATORDERINGNATURAL, &perm, &iperm);
 76:   switch (facttype) {
 77:   case MAT_FACTOR_LU:
 78:     MatLUFactorSymbolic(F, A[0], perm, iperm, &info);
 79:     break;
 80:   case MAT_FACTOR_ILU:
 81:     MatILUFactorSymbolic(F, A[0], perm, iperm, &info);
 82:     break;
 83:   case MAT_FACTOR_ICC:
 84:     MatICCFactorSymbolic(F, A[0], perm, &info);
 85:     break;
 86:   case MAT_FACTOR_CHOLESKY:
 87:     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:       MatLUFactorNumeric(F, A[k], &info);
 99:       break;
100:     case MAT_FACTOR_ICC:
101:     case MAT_FACTOR_CHOLESKY:
102:       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:     MatSolve(F, b, x);

111:     /* Check the residual */
112:     MatMult(A[k], x, u);
113:     VecAXPY(u, -1.0, b);
114:     VecNorm(u, NORM_INFINITY, &norm);
115:     if (norm > tol) 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++) MatDestroy(&A[k]);
120:   MatDestroy(&F);
121:   ISDestroy(&perm);
122:   ISDestroy(&iperm);
123:   VecDestroy(&x);
124:   VecDestroy(&b);
125:   VecDestroy(&u);
126:   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*/