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;
41: PetscInitialize(&argc,&args,(char*)0,help);
42: MPI_Comm_size(PETSC_COMM_WORLD,&size);
43: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
45: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
46: Compute the matrix and right-hand-side vector that define
47: the linear system, As x = b.
48: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
49: /* Create vectors */
50: VecCreate(PETSC_COMM_WORLD,&x);
51: PetscObjectSetName((PetscObject) x, "Solution");
52: VecSetSizes(x,PETSC_DECIDE,n);
53: VecSetFromOptions(x);
54: VecDuplicate(x,&b);
55: VecDuplicate(x,&u);
57: /* Create matrix P, to be used as preconditioner */
58: MatCreate(PETSC_COMM_WORLD,&P);
59: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,n,n);
60: MatSetFromOptions(P);
61: MatSetUp(P);
63: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
64: for (i=1; i<n-1; i++) {
65: col[0] = i-1; col[1] = i; col[2] = i+1;
66: MatSetValues(P,1,&i,3,col,value,INSERT_VALUES);
67: }
68: i = n - 1; col[0] = n - 2; col[1] = n - 1;
69: MatSetValues(P,1,&i,2,col,value,INSERT_VALUES);
70: i = 0; col[0] = 0; col[1] = 1; value[0] = 2.0; value[1] = -1.0;
71: MatSetValues(P,1,&i,2,col,value,INSERT_VALUES);
72: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
73: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
75: /* Set exact solution */
76: VecSet(u,one);
78: /* Create a matrix-free matrix As, P is used as a data context in MyMatShellMult() */
79: MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,P,&As);
80: MatSetFromOptions(As);
81: MatShellSetOperation(As,MATOP_MULT,(void(*)(void))MyMatShellMult);
83: /* Check As is a linear operator: As*(ax + y) = a As*x + As*y */
84: MatIsLinear(As,10,&flg);
87: /* Compute right-hand-side vector. */
88: MatMult(As,u,b);
90: MatSetOption(As,MAT_SYMMETRIC,PETSC_TRUE);
91: MatMultTranspose(As,u,x);
92: VecAXPY(x,-1.0,b);
93: VecNorm(x,NORM_INFINITY,&norm);
95: MatSetOption(As,MAT_HERMITIAN,PETSC_TRUE);
96: MatMultHermitianTranspose(As,u,x);
97: VecAXPY(x,-1.0,b);
98: VecNorm(x,NORM_INFINITY,&norm);
101: /* Create the linear solver and set various options */
102: KSPCreate(PETSC_COMM_WORLD,&ksp);
103: KSPSetOperators(ksp,As,P);
105: /* Set linear solver defaults for this problem (optional). */
106: KSPGetPC(ksp,&pc);
107: PCSetType(pc,PCNONE);
108: KSPSetTolerances(ksp,1.e-5,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
110: /* Set runtime options */
111: KSPSetFromOptions(ksp);
113: /* Solve linear system */
114: KSPSolve(ksp,b,x);
116: /* Check the error */
117: VecAXPY(x,-1.0,u);
118: VecNorm(x,NORM_2,&norm);
119: KSPGetIterationNumber(ksp,&its);
120: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);
122: /* Free work space. */
123: VecDestroy(&x)); PetscCall(VecDestroy(&u);
124: VecDestroy(&b)); PetscCall(MatDestroy(&P);
125: MatDestroy(&As);
126: KSPDestroy(&ksp);
128: PetscFinalize();
129: return 0;
130: }
132: /*TEST
134: test:
135: args: -ksp_monitor_short -ksp_max_it 10
136: test:
137: suffix: 2
138: args: -ksp_monitor_short -ksp_max_it 10
140: TEST*/