Actual source code: ex7.c


  2: static char help[] = "Illustrate how to solves a matrix-free linear system with KSP.\n\n";

  4: /*
  5:   Note: modified from ~src/ksp/ksp/tutorials/ex1.c
  6: */
  7: #include <petscksp.h>

  9: /*
 10:    MatShellMult - Computes the matrix-vector product, y = As x.

 12:    Input Parameters:
 13:    As - the matrix-free matrix
 14:    x  - vector

 16:    Output Parameter:
 17:    y - vector
 18:  */
 19: PetscErrorCode MyMatShellMult(Mat As, Vec x, Vec y)
 20: {
 21:   Mat P;

 23:   /* printf("MatShellMult...user should implement this routine without using a matrix\n"); */
 24:   MatShellGetContext(As, &P);
 25:   MatMult(P, x, y);
 26:   return 0;
 27: }

 29: int main(int argc, char **args)
 30: {
 31:   Vec         x, b, u; /* approx solution, RHS, exact solution */
 32:   Mat         P, As;   /* preconditioner matrix, linear system (matrix-free) */
 33:   KSP         ksp;     /* linear solver context */
 34:   PC          pc;      /* preconditioner context */
 35:   PetscReal   norm;    /* norm of solution error */
 36:   PetscInt    i, n = 100, col[3], its;
 37:   PetscMPIInt size;
 38:   PetscScalar one = 1.0, value[3];
 39:   PetscBool   flg;

 42:   PetscInitialize(&argc, &args, (char *)0, help);
 43:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 44:   PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);

 46:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 47:          Compute the matrix and right-hand-side vector that define
 48:          the linear system, As x = b.
 49:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 50:   /* Create vectors */
 51:   VecCreate(PETSC_COMM_WORLD, &x);
 52:   PetscObjectSetName((PetscObject)x, "Solution");
 53:   VecSetSizes(x, PETSC_DECIDE, n);
 54:   VecSetFromOptions(x);
 55:   VecDuplicate(x, &b);
 56:   VecDuplicate(x, &u);

 58:   /* Create matrix P, to be used as preconditioner */
 59:   MatCreate(PETSC_COMM_WORLD, &P);
 60:   MatSetSizes(P, PETSC_DECIDE, PETSC_DECIDE, n, n);
 61:   MatSetFromOptions(P);
 62:   MatSetUp(P);

 64:   value[0] = -1.0;
 65:   value[1] = 2.0;
 66:   value[2] = -1.0;
 67:   for (i = 1; i < n - 1; i++) {
 68:     col[0] = i - 1;
 69:     col[1] = i;
 70:     col[2] = i + 1;
 71:     MatSetValues(P, 1, &i, 3, col, value, INSERT_VALUES);
 72:   }
 73:   i      = n - 1;
 74:   col[0] = n - 2;
 75:   col[1] = n - 1;
 76:   MatSetValues(P, 1, &i, 2, col, value, INSERT_VALUES);
 77:   i        = 0;
 78:   col[0]   = 0;
 79:   col[1]   = 1;
 80:   value[0] = 2.0;
 81:   value[1] = -1.0;
 82:   MatSetValues(P, 1, &i, 2, col, value, INSERT_VALUES);
 83:   MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
 84:   MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);

 86:   /* Set exact solution */
 87:   VecSet(u, one);

 89:   /* Create a matrix-free matrix As, P is used as a data context in MyMatShellMult() */
 90:   MatCreateShell(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, n, P, &As);
 91:   MatSetFromOptions(As);
 92:   MatShellSetOperation(As, MATOP_MULT, (void (*)(void))MyMatShellMult);

 94:   /* Check As is a linear operator: As*(ax + y) = a As*x + As*y */
 95:   MatIsLinear(As, 10, &flg);

 98:   /* Compute right-hand-side vector. */
 99:   MatMult(As, u, b);

101:   MatSetOption(As, MAT_SYMMETRIC, PETSC_TRUE);
102:   MatMultTranspose(As, u, x);
103:   VecAXPY(x, -1.0, b);
104:   VecNorm(x, NORM_INFINITY, &norm);
106:   MatSetOption(As, MAT_HERMITIAN, PETSC_TRUE);
107:   MatMultHermitianTranspose(As, u, x);
108:   VecAXPY(x, -1.0, b);
109:   VecNorm(x, NORM_INFINITY, &norm);

112:   /* Create the linear solver and set various options */
113:   KSPCreate(PETSC_COMM_WORLD, &ksp);
114:   KSPSetOperators(ksp, As, P);

116:   /* Set linear solver defaults for this problem (optional). */
117:   KSPGetPC(ksp, &pc);
118:   PCSetType(pc, PCNONE);
119:   KSPSetTolerances(ksp, 1.e-5, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);

121:   /* Set runtime options */
122:   KSPSetFromOptions(ksp);

124:   /* Solve linear system */
125:   KSPSolve(ksp, b, x);

127:   /* Check the error */
128:   VecAXPY(x, -1.0, u);
129:   VecNorm(x, NORM_2, &norm);
130:   KSPGetIterationNumber(ksp, &its);
131:   PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its);

133:   /* Free work space. */
134:   VecDestroy(&x);
135:   VecDestroy(&u);
136:   VecDestroy(&b);
137:   MatDestroy(&P);
138:   MatDestroy(&As);
139:   KSPDestroy(&ksp);

141:   PetscFinalize();
142:   return 0;
143: }

145: /*TEST

147:    test:
148:       args: -ksp_monitor_short -ksp_max_it 10
149:    test:
150:       suffix: 2
151:       args: -ksp_monitor_short -ksp_max_it 10

153: TEST*/