Actual source code: ex6.c


  2: static char help[] = "Solves a tridiagonal linear system with KSP. \n\
  3: It illustrates how to do one symbolic factorization and multiple numeric factorizations using same matrix structure. \n\n";

  5: #include <petscksp.h>
  6: int main(int argc, char **args)
  7: {
  8:   Vec         x, b, u; /* approx solution, RHS, exact solution */
  9:   Mat         A;       /* linear system matrix */
 10:   KSP         ksp;     /* linear solver context */
 11:   PC          pc;      /* preconditioner context */
 12:   PetscReal   norm;    /* norm of solution error */
 13:   PetscInt    i, col[3], its, rstart, rend, N = 10, num_numfac;
 14:   PetscScalar value[3];

 17:   PetscInitialize(&argc, &args, (char *)0, help);
 18:   PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL);

 20:   /* Create and assemble matrix. */
 21:   MatCreate(PETSC_COMM_WORLD, &A);
 22:   MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, N, N);
 23:   MatSetFromOptions(A);
 24:   MatSetUp(A);
 25:   MatGetOwnershipRange(A, &rstart, &rend);

 27:   value[0] = -1.0;
 28:   value[1] = 2.0;
 29:   value[2] = -1.0;
 30:   for (i = rstart; i < rend; i++) {
 31:     col[0] = i - 1;
 32:     col[1] = i;
 33:     col[2] = i + 1;
 34:     if (i == 0) {
 35:       MatSetValues(A, 1, &i, 2, col + 1, value + 1, INSERT_VALUES);
 36:     } else if (i == N - 1) {
 37:       MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
 38:     } else {
 39:       MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);
 40:     }
 41:   }
 42:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 43:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 44:   MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);

 46:   /* Create vectors */
 47:   MatCreateVecs(A, &x, &b);
 48:   VecDuplicate(x, &u);

 50:   /* Set exact solution; then compute right-hand-side vector. */
 51:   VecSet(u, 1.0);
 52:   MatMult(A, u, b);

 54:   /* Create the linear solver and set various options. */
 55:   KSPCreate(PETSC_COMM_WORLD, &ksp);
 56:   KSPGetPC(ksp, &pc);
 57:   PCSetType(pc, PCJACOBI);
 58:   KSPSetTolerances(ksp, 1.e-5, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
 59:   KSPSetOperators(ksp, A, A);
 60:   KSPSetFromOptions(ksp);

 62:   num_numfac = 1;
 63:   PetscOptionsGetInt(NULL, NULL, "-num_numfac", &num_numfac, NULL);
 64:   while (num_numfac--) {
 65:     /* An example on how to update matrix A for repeated numerical factorization and solve. */
 66:     PetscScalar one = 1.0;
 67:     PetscInt    i   = 0;
 68:     MatSetValues(A, 1, &i, 1, &i, &one, ADD_VALUES);
 69:     MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 70:     MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 71:     /* Update b */
 72:     MatMult(A, u, b);

 74:     /* Solve the linear system */
 75:     KSPSolve(ksp, b, x);

 77:     /* Check the solution and clean up */
 78:     VecAXPY(x, -1.0, u);
 79:     VecNorm(x, NORM_2, &norm);
 80:     KSPGetIterationNumber(ksp, &its);
 81:     if (norm > 100 * PETSC_MACHINE_EPSILON) PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its);
 82:   }

 84:   /* Free work space. */
 85:   VecDestroy(&x);
 86:   VecDestroy(&u);
 87:   VecDestroy(&b);
 88:   MatDestroy(&A);
 89:   KSPDestroy(&ksp);

 91:   PetscFinalize();
 92:   return 0;
 93: }

 95: /*TEST

 97:     test:
 98:       args: -num_numfac 2 -pc_type lu

100:     test:
101:       suffix: 2
102:       args: -num_numfac 2 -pc_type lu -pc_factor_mat_solver_type mumps
103:       requires: mumps

105:     test:
106:       suffix: 3
107:       nsize: 3
108:       args: -num_numfac 2 -pc_type lu -pc_factor_mat_solver_type mumps
109:       requires: mumps

111: TEST*/