Actual source code: ex7.c

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

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

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

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

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

 22:   PetscFunctionBegin;
 23:   /* printf("MatShellMult...user should implement this routine without using a matrix\n"); */
 24:   PetscCall(MatShellGetContext(As, &P));
 25:   PetscCall(MatMult(P, x, y));
 26:   PetscFunctionReturn(PETSC_SUCCESS);
 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;

 41:   PetscFunctionBeginUser;
 42:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 43:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 44:   PetscCall(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:   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
 52:   PetscCall(PetscObjectSetName((PetscObject)x, "Solution"));
 53:   PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
 54:   PetscCall(VecSetFromOptions(x));
 55:   PetscCall(VecDuplicate(x, &b));
 56:   PetscCall(VecDuplicate(x, &u));

 58:   /* Create matrix P, to be used as preconditioner */
 59:   PetscCall(MatCreate(PETSC_COMM_WORLD, &P));
 60:   PetscCall(MatSetSizes(P, PETSC_DECIDE, PETSC_DECIDE, n, n));
 61:   PetscCall(MatSetFromOptions(P));
 62:   PetscCall(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:     PetscCall(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:   PetscCall(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:   PetscCall(MatSetValues(P, 1, &i, 2, col, value, INSERT_VALUES));
 83:   PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
 84:   PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));

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

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

 94:   /* Check As is a linear operator: As*(ax + y) = a As*x + As*y */
 95:   PetscCall(MatIsLinear(As, 10, &flg));
 96:   PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Shell matrix As is non-linear! Use '-info |grep MatIsLinear' to get detailed report");

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

101:   PetscCall(MatSetOption(As, MAT_SYMMETRIC, PETSC_TRUE));
102:   PetscCall(MatMultTranspose(As, u, x));
103:   PetscCall(VecAXPY(x, -1.0, b));
104:   PetscCall(VecNorm(x, NORM_INFINITY, &norm));
105:   PetscCheck(norm <= PETSC_SMALL, PetscObjectComm((PetscObject)As), PETSC_ERR_PLIB, "Error ||A x-A^T x||_\\infty: %1.6e", (double)norm);
106:   PetscCall(MatSetOption(As, MAT_HERMITIAN, PETSC_TRUE));
107:   PetscCall(MatMultHermitianTranspose(As, u, x));
108:   PetscCall(VecAXPY(x, -1.0, b));
109:   PetscCall(VecNorm(x, NORM_INFINITY, &norm));
110:   PetscCheck(norm <= PETSC_SMALL, PetscObjectComm((PetscObject)As), PETSC_ERR_PLIB, "Error ||A x-A^H x||_\\infty: %1.6e", (double)norm);

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

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

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

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

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

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

141:   PetscCall(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*/