Actual source code: ex2.c
1: /*$Id: ex2.c,v 1.26 2001/08/07 21:30:50 bsmith Exp $*/
3: static char help[] = "Demonstrates running several independent tasks in PETSc.\n\n";
5: /*T
6: Concepts: SLES^solving linear equations
7: Processors: n
9: Comments: Demonstrates how to use PetscSetCommWorld() to tell a subset of
10: processors (in this case each individual processor) to run
11: as if it was all the processors that PETSc is using. ADVANCED
12: example, not for beginning PETSc users.
14: T*/
16: /*
17: Include "petscsles.h" so that we can use SLES solvers. Note that this file
18: automatically includes:
19: petsc.h - base PETSc routines petscvec.h - vectors
20: petscsys.h - system routines petscmat.h - matrices
21: petscis.h - index sets petscksp.h - Krylov subspace methods
22: petscviewer.h - viewers petscpc.h - preconditioners
23: */
24: #include petscsles.h
26: EXTERN int slesex(int,char**);
30: int main(int argc,char **argv)
31: {
32: MPI_Init(&argc,&argv);
33: slesex(argc,argv);
34: MPI_Finalize();
35: return 0;
36: }
40: int slesex(int argc,char **args)
41: {
42: Vec x,b,u; /* approx solution, RHS, exact solution */
43: Mat A; /* linear system matrix */
44: SLES sles; /* linear solver context */
45: PC pc; /* preconditioner context */
46: KSP ksp; /* Krylov subspace method context */
47: PetscReal norm; /* norm of solution error */
48: int i,j,I,J,Istart,Iend,ierr,m = 8,n = 7,its;
49: PetscScalar v,one = 1.0,none = -1.0;
51: PetscSetCommWorld(PETSC_COMM_SELF);
52: PetscInitialize(&argc,&args,(char *)0,help);
54: PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);
55: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
57: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
58: Compute the matrix and right-hand-side vector that define
59: the linear system, Ax = b.
60: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
61: /*
62: Create parallel matrix, specifying only its global dimensions.
63: When using MatCreate(), the matrix format can be specified at
64: runtime. Also, the parallel partitioning of the matrix is
65: determined by PETSc at runtime.
66: */
67: MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n,&A);
68: MatSetFromOptions(A);
70: /*
71: Currently, all PETSc parallel matrix formats are partitioned by
72: contiguous chunks of rows across the processors. Determine which
73: rows of the matrix are locally owned.
74: */
75: MatGetOwnershipRange(A,&Istart,&Iend);
77: /*
78: Set matrix elements for the 2-D, five-point stencil in parallel.
79: - Each processor needs to insert only elements that it owns
80: locally (but any non-local elements will be sent to the
81: appropriate processor during matrix assembly).
82: - Always specify global row and columns of matrix entries.
83: */
84: for (I=Istart; I<Iend; I++) {
85: v = -1.0; i = I/n; j = I - i*n;
86: if (i>0) {J = I - n; MatSetValues(A,1,&I,1,&J,&v,INSERT_VALUES);}
87: if (i<m-1) {J = I + n; MatSetValues(A,1,&I,1,&J,&v,INSERT_VALUES);}
88: if (j>0) {J = I - 1; MatSetValues(A,1,&I,1,&J,&v,INSERT_VALUES);}
89: if (j<n-1) {J = I + 1; MatSetValues(A,1,&I,1,&J,&v,INSERT_VALUES);}
90: v = 4.0; MatSetValues(A,1,&I,1,&I,&v,INSERT_VALUES);
91: }
93: /*
94: Assemble matrix, using the 2-step process:
95: MatAssemblyBegin(), MatAssemblyEnd()
96: Computations can be done while messages are in transition,
97: by placing code between these two statements.
98: */
99: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
100: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
102: /*
103: Create parallel vectors.
104: - When using VecCreate() VecSetSizes() and VecSetFromOptions(),
105: we specify only the vector's global
106: dimension; the parallel partitioning is determined at runtime.
107: - Note: We form 1 vector from scratch and then duplicate as needed.
108: */
109: VecCreate(PETSC_COMM_WORLD,&u);
110: VecSetSizes(u,PETSC_DECIDE,m*n);
111: VecSetFromOptions(u);
112: VecDuplicate(u,&b);
113: VecDuplicate(b,&x);
115: /*
116: Set exact solution; then compute right-hand-side vector.
117: */
118: VecSet(&one,u);
119: MatMult(A,u,b);
121: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: Create the linear solver and set various options
123: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125: /*
126: Create linear solver context
127: */
128: SLESCreate(PETSC_COMM_WORLD,&sles);
130: /*
131: Set operators. Here the matrix that defines the linear system
132: also serves as the preconditioning matrix.
133: */
134: SLESSetOperators(sles,A,A,DIFFERENT_NONZERO_PATTERN);
136: /*
137: Set linear solver defaults for this problem (optional).
138: - By extracting the KSP and PC contexts from the SLES context,
139: we can then directly directly call any KSP and PC routines
140: to set various options.
141: - The following four statements are optional; all of these
142: parameters could alternatively be specified at runtime via
143: SLESSetFromOptions();
144: */
145: SLESGetKSP(sles,&ksp);
146: SLESGetPC(sles,&pc);
147: PCSetType(pc,PCJACOBI);
148: KSPSetTolerances(ksp,1.e-7,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
150: /*
151: Set runtime options, e.g.,
152: -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
153: These options will override those specified above as long as
154: SLESSetFromOptions() is called _after_ any other customization
155: routines.
156: */
157: SLESSetFromOptions(sles);
159: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160: Solve the linear system
161: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
163: SLESSolve(sles,b,x,&its);
165: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166: Check solution and clean up
167: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169: /*
170: Check the error
171: */
172: VecAXPY(&none,u,x);
173: VecNorm(x,NORM_2,&norm);
174: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A iterations %d\n",norm,its);
176: /*
177: Free work space. All PETSc objects should be destroyed when they
178: are no longer needed.
179: */
180: SLESDestroy(sles);
181: VecDestroy(u); VecDestroy(x);
182: VecDestroy(b); MatDestroy(A);
183: PetscFinalize();
184: return 0;
185: }