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