Actual source code: ex28.c
petsc-3.4.2 2013-07-02
2: static char help[] = "Test procedural KSPSetFromOptions() or at runtime.\n\n";
4: /*T
5: Concepts: KSP^basic parallel example;
6: Processors: n
7: T*/
8: #include <petscksp.h>
12: int main(int argc,char **args)
13: {
14: Vec x, b, u; /* approx solution, RHS, exact solution */
15: Mat A; /* linear system matrix */
16: KSP ksp; /* linear solver context */
17: PC pc; /* preconditioner context */
18: PetscReal norm; /* norm of solution error */
20: PetscInt i,n = 10,col[3],its,rstart,rend,nlocal;
21: PetscScalar neg_one = -1.0,one = 1.0,value[3];
22: PetscBool TEST_PROCEDURAL=PETSC_FALSE;
24: PetscInitialize(&argc,&args,(char*)0,help);
25: PetscOptionsGetInt(NULL,"-n",&n,NULL);
26: PetscOptionsGetBool(NULL,"-procedural",&TEST_PROCEDURAL,NULL);
28: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
29: Compute the matrix and right-hand-side vector that define
30: the linear system, Ax = b.
31: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
33: /*
34: Create vectors. Note that we form 1 vector from scratch and
35: then duplicate as needed. For this simple case let PETSc decide how
36: many elements of the vector are stored on each processor. The second
37: argument to VecSetSizes() below causes PETSc to decide.
38: */
39: VecCreate(PETSC_COMM_WORLD,&x);
40: VecSetSizes(x,PETSC_DECIDE,n);
41: VecSetFromOptions(x);
42: VecDuplicate(x,&b);
43: VecDuplicate(x,&u);
45: /* Identify the starting and ending mesh points on each
46: processor for the interior part of the mesh. We let PETSc decide
47: above. */
49: VecGetOwnershipRange(x,&rstart,&rend);
50: VecGetLocalSize(x,&nlocal);
52: /* Create a tridiagonal matrix. See ../tutorials/ex23.c */
53: MatCreate(PETSC_COMM_WORLD,&A);
54: MatSetSizes(A,nlocal,nlocal,n,n);
55: MatSetFromOptions(A);
56: /* Assemble matrix */
57: if (!rstart) {
58: rstart = 1;
59: i = 0; col[0] = 0; col[1] = 1; value[0] = 2.0; value[1] = -1.0;
60: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
61: }
62: if (rend == n) {
63: rend = n-1;
64: i = n-1; col[0] = n-2; col[1] = n-1; value[0] = -1.0; value[1] = 2.0;
65: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
66: }
68: /* Set entries corresponding to the mesh interior */
69: value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
70: for (i=rstart; i<rend; i++) {
71: col[0] = i-1; col[1] = i; col[2] = i+1;
72: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
73: }
75: /* Assemble the matrix */
76: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
77: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
79: /* Set exact solution; then compute right-hand-side vector. */
80: VecSet(u,one);
81: MatMult(A,u,b);
83: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: Create the linear solver and set various options
85: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
86: KSPCreate(PETSC_COMM_WORLD,&ksp);
87: KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
89: /*
90: Set linear solver defaults for this problem (optional).
91: - By extracting the KSP and PC contexts from the KSP context,
92: we can then directly call any KSP and PC routines to set
93: various options.
94: - The following statements are optional; all of these
95: parameters could alternatively be specified at runtime via
96: KSPSetFromOptions();
97: */
98: if (TEST_PROCEDURAL) {
99: PetscMPIInt size;
100: KSPGetPC(ksp,&pc);
101: PCSetType(pc,PCREDUNDANT);
102: MPI_Comm_size(PETSC_COMM_WORLD,&size);
103: if (size < 3) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Num of processes %d must greater than 2",size);
104: PCRedundantSetNumber(pc,size-2);
105: } else {
106: KSPSetFromOptions(ksp);
107: }
108: /* Solve linear system */
109: KSPSolve(ksp,b,x);
111: if (TEST_PROCEDURAL) { /* View solver info; */
112: KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);
113: }
115: /* Check the error */
116: VecAXPY(x,neg_one,u);
117: VecNorm(x,NORM_2,&norm);
118: KSPGetIterationNumber(ksp,&its);
119: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %G, Iterations %D\n",norm,its);
121: /* Free work space. */
122: VecDestroy(&x); VecDestroy(&u);
123: VecDestroy(&b); MatDestroy(&A);
124: KSPDestroy(&ksp);
125: PetscFinalize();
126: return 0;
127: }