Actual source code: ex7.c
petsc-3.13.4 2020-08-01
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: 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: MatSetOption(As,MAT_SYMMETRIC,PETSC_TRUE);
94: MatMultTranspose(As,u,x);
95: VecAXPY(x,-1.0,b);
96: VecNorm(x,NORM_INFINITY,&norm);
97: if (norm > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)As),PETSC_ERR_PLIB,"Error ||A x-A^T x||_\\infty: %1.6e",norm);
98: MatSetOption(As,MAT_HERMITIAN,PETSC_TRUE);
99: MatMultHermitianTranspose(As,u,x);
100: VecAXPY(x,-1.0,b);
101: VecNorm(x,NORM_INFINITY,&norm);
102: if (norm > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)As),PETSC_ERR_PLIB,"Error ||A x-A^H x||_\\infty: %1.6e",norm);
104: /* Create the linear solver and set various options */
105: KSPCreate(PETSC_COMM_WORLD,&ksp);
106: KSPSetOperators(ksp,As,P);
108: /* Set linear solver defaults for this problem (optional). */
109: KSPGetPC(ksp,&pc);
110: PCSetType(pc,PCNONE);
111: KSPSetTolerances(ksp,1.e-5,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
113: /* Set runtime options */
114: KSPSetFromOptions(ksp);
116: /* Solve linear system */
117: KSPSolve(ksp,b,x);
119: /* Check the error */
120: VecAXPY(x,-1.0,u);
121: VecNorm(x,NORM_2,&norm);
122: KSPGetIterationNumber(ksp,&its);
123: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);
125: /* Free work space. */
126: VecDestroy(&x); VecDestroy(&u);
127: VecDestroy(&b); MatDestroy(&P);
128: MatDestroy(&As);
129: KSPDestroy(&ksp);
131: PetscFinalize();
132: return ierr;
133: }
135: /*TEST
137: test:
138: args: -ksp_monitor_short -ksp_max_it 10
139: test:
140: suffix: 2
141: args: -ksp_monitor_short -ksp_max_it 10
143: TEST*/