Actual source code: ex7.c
petsc-3.12.0 2019-09-29
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/examples/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: PetscErrorCode ierr;
22: Mat P;
25: /* printf("MatShellMult...user should implement this routine without using a matrix\n"); */
26: MatShellGetContext(As,&P);
27: MatMult(P,x,y);
28: return(0);
29: }
31: int main(int argc,char **args)
32: {
33: Vec x, b, u; /* approx solution, RHS, exact solution */
34: Mat P,As; /* preconditioner matrix, linear system (matrix-free) */
35: KSP ksp; /* linear solver context */
36: PC pc; /* preconditioner context */
37: PetscReal norm; /* norm of solution error */
39: PetscInt i,n = 100,col[3],its;
40: PetscMPIInt size;
41: PetscScalar one = 1.0,value[3];
42: PetscBool flg;
44: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
45: MPI_Comm_size(PETSC_COMM_WORLD,&size);
46: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
48: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
49: Compute the matrix and right-hand-side vector that define
50: the linear system, As x = b.
51: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
52: /* Create vectors */
53: VecCreate(PETSC_COMM_WORLD,&x);
54: PetscObjectSetName((PetscObject) x, "Solution");
55: VecSetSizes(x,PETSC_DECIDE,n);
56: VecSetFromOptions(x);
57: VecDuplicate(x,&b);
58: VecDuplicate(x,&u);
60: /* Create matrix P, to be used as preconditioner */
61: MatCreate(PETSC_COMM_WORLD,&P);
62: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,n,n);
63: MatSetFromOptions(P);
64: MatSetUp(P);
66: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
67: for (i=1; i<n-1; i++) {
68: col[0] = i-1; col[1] = i; col[2] = i+1;
69: MatSetValues(P,1,&i,3,col,value,INSERT_VALUES);
70: }
71: i = n - 1; col[0] = n - 2; col[1] = n - 1;
72: MatSetValues(P,1,&i,2,col,value,INSERT_VALUES);
73: i = 0; col[0] = 0; col[1] = 1; value[0] = 2.0; value[1] = -1.0;
74: MatSetValues(P,1,&i,2,col,value,INSERT_VALUES);
75: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
76: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
78: /* Set exact solution */
79: VecSet(u,one);
81: /* Create a matrix-free matrix As, P is used as a data context in MyMatShellMult() */
82: MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,P,&As);
83: MatSetFromOptions(As);
84: MatShellSetOperation(As,MATOP_MULT,(void(*)(void))MyMatShellMult);
86: /* Check As is a linear operator: As*(ax + y) = a As*x + As*y */
87: MatIsLinear(As,10,&flg);
88: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Shell matrix As is non-linear! Use '-info |grep MatIsLinear' to get detailed report\n");
90: /* Compute right-hand-side vector. */
91: MatMult(As,u,b);
93: /* Create the linear solver and set various options */
94: KSPCreate(PETSC_COMM_WORLD,&ksp);
95: KSPSetOperators(ksp,As,P);
97: /* Set linear solver defaults for this problem (optional). */
98: KSPGetPC(ksp,&pc);
99: PCSetType(pc,PCNONE);
100: KSPSetTolerances(ksp,1.e-5,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
102: /* Set runtime options */
103: KSPSetFromOptions(ksp);
105: /* Solve linear system */
106: KSPSolve(ksp,b,x);
108: /* Check the error */
109: VecAXPY(x,-1.0,u);
110: VecNorm(x,NORM_2,&norm);
111: KSPGetIterationNumber(ksp,&its);
112: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);
114: /* Free work space. */
115: VecDestroy(&x); VecDestroy(&u);
116: VecDestroy(&b); MatDestroy(&P);
117: MatDestroy(&As);
118: KSPDestroy(&ksp);
120: PetscFinalize();
121: return ierr;
122: }
124: /*TEST
126: test:
127: args: -ksp_monitor_short -ksp_max_it 10
128: test:
129: suffix: 2
130: args: -ksp_monitor_short -ksp_max_it 10
132: TEST*/