← Back to branch summary

~chaffra/+junk/trilinos

« back to all changes in this revision

Viewing changes to packages/trilinoscouplings/examples/epetraext/EpetraExt_petsc.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme, Christophe Prud'homme, Johannes Ring
  • Date: 2009-12-13 12:53:22 UTC
  • mfrom: (5.1.2 sid)
  • Revision ID: james.westby@ubuntu.com-20091213125322-in0nrdjc55deqsw9
Tags: 10.0.3.dfsg-1
[Christophe Prud'homme]
* New upstream release

[Johannes Ring]
* debian/patches/libname.patch: Add prefix 'libtrilinos_' to all
  libraries. 
* debian/patches/soname.patch: Add soversion to libraries.
* debian/watch: Update download URL.
* debian/control:
  - Remove python-numeric from Build-Depends (virtual package).
  - Remove automake and autotools from Build-Depends and add cmake to
    reflect switch to CMake.
  - Add python-support to Build-Depends.
* debian/rules: 
  - Cleanup and updates for switch to CMake.

Show diffs side-by-side

added added

removed removed

Lines of Context:
packages/trilinoscouplings/examples/epetraext/EpetraExt_petsc.cpp
 
1
#include "ml_config.h"
 
2
#include "EpetraExt_config.h"
 
3
#if defined(HAVE_PETSC) && defined(HAVE_ML_EPETRA) && defined(HAVE_ML_TEUCHOS) && defined(HAVE_ML_AZTECOO) && defined(HAVE_MPI)
 
4
#include "petscksp.h"
 
5
#include "EpetraExt_PETScAIJMatrix.h"
 
6
#include "Epetra_Vector.h"
 
7
#include "Epetra_Map.h"
 
8
#include "ml_MultiLevelPreconditioner.h"
 
9
#include "AztecOO.h"
 
10
#include "Epetra_LinearProblem.h"
 
11
 
 
12
/* 
 
13
   This example demonstrates how to apply a Trilinos preconditioner to a PETSc
 
14
   linear system.
 
15
 
 
16
   For information on configuring and building Trilinos with the PETSc aij
 
17
   interface enabled, please see EpetraExt's doxygen documentation at
 
18
   http://trilinos.sandia.gov/packages/epetraext, development version
 
19
   or release 9.0 or later.
 
20
 
 
21
   The PETSc matrix is a 2D 5-point Laplace operator stored in AIJ format.
 
22
   This matrix is wrapped as an Epetra_PETScAIJMatrix, and an ML AMG
 
23
   preconditioner is created for it.  The associated linear system is solved twice,
 
24
   the first time using AztecOO's preconditioned CG, the second time using PETSc's.
 
25
 
 
26
   To invoke this example, use something like:
 
27
 
 
28
       mpirun -np 5 ./TrilinosCouplings_petsc.exe -m 150 -n 150 -petsc_smoother -ksp_monitor_true_residual
 
29
*/
 
30
 
 
31
static char help[] = "Demonstrates how to solve a PETSc linear system with KSP\
 
32
and a Trilinos AMG preconditioner.  In particular, it shows how to wrap a PETSc\
 
33
AIJ matrix as an Epetra matrix, create the AMG preconditioner, and wrap it as a\
 
34
shell preconditioner for a PETSc Krylov method.\
 
35
Input parameters include:\n\
 
36
  -random_exact_sol : use a random exact solution vector\n\
 
37
  -view_exact_sol   : write exact solution vector to stdout\n\
 
38
  -m <mesh_x>       : number of mesh points in x-direction\n\
 
39
  -n <mesh_n>       : number of mesh points in y-direction\n\n";
 
40
 
 
41
extern PetscErrorCode ShellApplyML(void*,Vec,Vec);
 
42
 
 
43
#undef __FUNCT__
 
44
#define __FUNCT__ "main"
 
45
int main(int argc,char **args)
 
46
{
 
47
  Vec            x,b,u;  /* approx solution, RHS, exact solution */
 
48
  Mat            A;        /* linear system matrix */
 
49
  KSP            ksp;     /* linear solver context */
 
50
  KSP            kspSmoother=0;  /* solver context for PETSc fine grid smoother */
 
51
  PC pc;
 
52
  PetscRandom    rctx;     /* random number generator context */
 
53
  PetscReal      norm;     /* norm of solution error */
 
54
  PetscInt       i,j,Ii,J,Istart,Iend,its;
 
55
  PetscInt       m = 50,n = 50; /* #mesh points in x & y directions, resp. */
 
56
  PetscErrorCode ierr;
 
57
  PetscTruth     flg;
 
58
  PetscScalar    v,one = 1.0,neg_one = -1.0;
 
59
  PetscInt rank=0;
 
60
  MPI_Comm comm;
 
61
 
 
62
  PetscInitialize(&argc,&args,(char *)0,help);
 
63
  ierr = PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);CHKERRQ(ierr);
 
64
  ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr);
 
65
 
 
66
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
 
67
  ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);CHKERRQ(ierr);
 
68
  ierr = MatSetType(A, MATAIJ);CHKERRQ(ierr);
 
69
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
 
70
  ierr = MatMPIAIJSetPreallocation(A,5,PETSC_NULL,5,PETSC_NULL);CHKERRQ(ierr);
 
71
  PetscObjectGetComm( (PetscObject)A, &comm);
 
72
  ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
 
73
  if (!rank) printf("Matrix has %d (%dx%d) rows\n",m*n,m,n);
 
74
 
 
75
  ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);
 
76
 
 
77
  for (Ii=Istart; Ii<Iend; Ii++) { 
 
78
    v = -1.0; i = Ii/n; j = Ii - i*n;  
 
79
    if (i>0)   {J = Ii - n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
 
80
    if (i<m-1) {J = Ii + n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
 
81
    if (j>0)   {J = Ii - 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
 
82
    if (j<n-1) {J = Ii + 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
 
83
    v = 4.0; ierr = MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);CHKERRQ(ierr);
 
84
  }
 
85
 
 
86
  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
 
87
  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
 
88
 
 
89
  /* Wrap the PETSc matrix as an Epetra_PETScAIJMatrix. This is lightweight,
 
90
     i.e., no deep data copies. */
 
91
  Epetra_PETScAIJMatrix epA(A);
 
92
 
 
93
  ierr = VecCreate(PETSC_COMM_WORLD,&u);CHKERRQ(ierr);
 
94
  ierr = VecSetSizes(u,PETSC_DECIDE,m*n);CHKERRQ(ierr);
 
95
  ierr = VecSetFromOptions(u);CHKERRQ(ierr);
 
96
  ierr = VecDuplicate(u,&b);CHKERRQ(ierr); 
 
97
 
 
98
  ierr = VecDuplicate(b,&x);CHKERRQ(ierr);
 
99
 
 
100
  ierr = PetscOptionsHasName(PETSC_NULL,"-random_exact_sol",&flg);CHKERRQ(ierr);
 
101
  if (flg) {
 
102
    ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
 
103
    ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
 
104
    ierr = VecSetRandom(u,rctx);CHKERRQ(ierr);
 
105
    ierr = PetscRandomDestroy(rctx);CHKERRQ(ierr);
 
106
  } else {
 
107
    ierr = VecSet(u,one);CHKERRQ(ierr);
 
108
  }
 
109
  ierr = MatMult(A,u,b);CHKERRQ(ierr);
 
110
  ierr = VecNorm(b,NORM_2,&norm);CHKERRQ(ierr);
 
111
  if (rank==0) printf("||b|| = %f\n",norm);
 
112
 
 
113
  /* Copy the PETSc vectors to Epetra vectors. */
 
114
  Epetra_Vector epu(epA.Map()), epb(epA.Map());
 
115
  PetscScalar *vals;
 
116
  ierr = VecGetArray(u,&vals);CHKERRQ(ierr);
 
117
  PetscInt length;
 
118
  ierr = VecGetLocalSize(u,&length);CHKERRQ(ierr);
 
119
  for (int j=0; j<length; j++) epu[j] = vals[j];
 
120
  ierr = VecRestoreArray(u,&vals);CHKERRQ(ierr);
 
121
  epA.Multiply(false, epu, epb);
 
122
 
 
123
  /* Check norms of the Epetra and PETSc vectors. */
 
124
  epu.Norm2(&norm);
 
125
  if (rank == 0) printf("||epetra u||_2 = %f\n",norm);
 
126
  ierr = VecNorm(u,NORM_2,&norm);CHKERRQ(ierr);
 
127
  if (rank == 0) printf("||petsc u||_2  = %f\n",norm);
 
128
  epb.Norm2(&norm);
 
129
  if (rank == 0) printf("||epetra b||_2 = %f\n",norm);
 
130
  ierr = VecNorm(b,NORM_2,&norm);CHKERRQ(ierr);
 
131
  if (rank == 0) printf("||petsc b||_2  = %f\n",norm);
 
132
 
 
133
  /* Create the ML AMG preconditioner. */
 
134
 
 
135
  /* Parameter list that holds options for AMG preconditioner. */
 
136
  Teuchos::ParameterList mlList;
 
137
  /* Set recommended defaults for Poisson-like problems. */
 
138
  ML_Epetra::SetDefaults("SA",mlList);
 
139
  /* Specify how much information ML prints to screen.
 
140
     0 is the minimum (no output), 10 is the maximum. */
 
141
  mlList.set("ML output",10);
 
142
  /* Set the fine grid smoother.  PETSc will be much faster for any
 
143
     smoother requiring row access, e.g., SOR.  For any smoother whose
 
144
     kernel is a matvec, Trilinos/PETSc performance should be comparable,
 
145
     as Trilinos simply calls the PETSc matvec.
 
146
 
 
147
     To use a PETSc smoother, create a KSP object, set the KSP type to
 
148
     KSPRICHARDSON, and set the desired smoother as the KSP preconditioner.
 
149
     It is important that you call KSPSetInitialGuessNonzero.  Otherwise, the
 
150
     post-smoother phase will incorrectly ignore the current approximate
 
151
     solution.  The KSP pointer must be cast to void* and passed to ML via
 
152
     the parameter list.
 
153
 
 
154
     You are responsible for freeing the KSP object.
 
155
  */
 
156
  ierr = PetscOptionsHasName(PETSC_NULL,"-petsc_smoother",&flg);CHKERRQ(ierr);
 
157
  if (flg) {
 
158
    ierr = KSPCreate(comm,&kspSmoother);CHKERRQ(ierr);
 
159
    ierr = KSPSetOperators(kspSmoother,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
 
160
    ierr = KSPSetType(kspSmoother,KSPRICHARDSON);CHKERRQ(ierr);
 
161
    ierr = KSPSetTolerances(kspSmoother, 1e-12, 1e-50, 1e7,1);
 
162
    ierr = KSPSetInitialGuessNonzero(kspSmoother,PETSC_TRUE);CHKERRQ(ierr);
 
163
    ierr = KSPGetPC(kspSmoother,&pc);CHKERRQ(ierr);
 
164
    ierr = PCSetType(pc, PCSOR);CHKERRQ(ierr);
 
165
    ierr = PCSORSetSymmetric(pc,SOR_LOCAL_SYMMETRIC_SWEEP);CHKERRQ(ierr);
 
166
    ierr = PCSetFromOptions(pc);CHKERRQ(ierr);
 
167
    ierr = KSPSetUp(kspSmoother);CHKERRQ(ierr);
 
168
    mlList.set("smoother: type (level 0)","petsc");
 
169
    mlList.set("smoother: petsc ksp (level 0)",(void*)kspSmoother);
 
170
  } else {
 
171
    mlList.set("smoother: type (level 0)","symmetric Gauss-Seidel");
 
172
  }
 
173
 
 
174
  /* how many fine grid pre- or post-smoothing sweeps to do */
 
175
  mlList.set("smoother: sweeps (level 0)",2);
 
176
 
 
177
  ML_Epetra::MultiLevelPreconditioner *Prec = new ML_Epetra::MultiLevelPreconditioner(epA,mlList);
 
178
 
 
179
  /* Trilinos CG */
 
180
  epu.PutScalar(0.0);
 
181
  Epetra_LinearProblem Problem(&epA, &epu, &epb);
 
182
  AztecOO solver(Problem);
 
183
  solver.SetAztecOption(AZ_solver, AZ_cg);
 
184
  solver.SetAztecOption(AZ_output, 1);
 
185
  solver.SetAztecOption(AZ_conv, AZ_noscaled);
 
186
  solver.SetPrecOperator(Prec);
 
187
  solver.Iterate(30, 1e-12);
 
188
 
 
189
  /* PETSc CG */
 
190
  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
 
191
  ierr = KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
 
192
  ierr = KSPSetTolerances(ksp,1e-12,1.e-50,PETSC_DEFAULT,
 
193
                          PETSC_DEFAULT);CHKERRQ(ierr);
 
194
  ierr = KSPSetType(ksp,KSPCG);CHKERRQ(ierr);
 
195
 
 
196
  /* Wrap the ML preconditioner as a PETSc shell preconditioner. */
 
197
  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
 
198
  ierr = PCSetType(pc,PCSHELL);CHKERRQ(ierr);
 
199
  ierr = PCShellSetApply(pc,ShellApplyML);CHKERRQ(ierr);
 
200
  ierr = PCShellSetContext(pc,(void*)Prec);CHKERRQ(ierr);
 
201
  ierr = PCShellSetName(pc,"ML AMG");CHKERRQ(ierr);
 
202
 
 
203
  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
 
204
 
 
205
  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
 
206
 
 
207
  ierr = KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
 
208
 
 
209
  ierr = VecAXPY(x,neg_one,u);CHKERRQ(ierr);
 
210
  ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
 
211
  ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
 
212
 
 
213
  ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A iterations %D\n",
 
214
                     norm,its);CHKERRQ(ierr);
 
215
 
 
216
  ierr = KSPDestroy(ksp);CHKERRQ(ierr);
 
217
  ierr = VecDestroy(u);CHKERRQ(ierr);  ierr = VecDestroy(x);CHKERRQ(ierr);
 
218
  ierr = VecDestroy(b);CHKERRQ(ierr);  ierr = MatDestroy(A);CHKERRQ(ierr);
 
219
 
 
220
  if (kspSmoother) {ierr = KSPDestroy(kspSmoother);CHKERRQ(ierr);}
 
221
  if (Prec) delete Prec;
 
222
 
 
223
  ierr = PetscFinalize();CHKERRQ(ierr);
 
224
  return 0;
 
225
} /*main*/
 
226
 
 
227
/* ***************************************************************** */
 
228
 
 
229
PetscErrorCode ShellApplyML(void *ctx,Vec x,Vec y)
 
230
{
 
231
  PetscErrorCode  ierr;
 
232
  ML_Epetra::MultiLevelPreconditioner *mlp = (ML_Epetra::MultiLevelPreconditioner*)ctx;
 
233
 
 
234
  /* Wrap x and y as Epetra_Vectors. */
 
235
  PetscScalar *xvals,*yvals;
 
236
  ierr = VecGetArray(x,&xvals);CHKERRQ(ierr);
 
237
  Epetra_Vector epx(View,mlp->OperatorDomainMap(),xvals);
 
238
  ierr = VecGetArray(y,&yvals);CHKERRQ(ierr);
 
239
  Epetra_Vector epy(View,mlp->OperatorRangeMap(),yvals);
 
240
 
 
241
  /* Apply ML. */
 
242
  mlp->ApplyInverse(epx,epy);
 
243
  
 
244
  /* Clean up and return. */
 
245
  ierr = VecRestoreArray(x,&xvals);CHKERRQ(ierr);
 
246
  ierr = VecRestoreArray(y,&yvals);CHKERRQ(ierr);
 
247
  return 0;
 
248
} /*ShellApplyML*/
 
249
#else
 
250
 
 
251
#include <stdlib.h>
 
252
#include <stdio.h>
 
253
#ifdef HAVE_MPI
 
254
#include "mpi.h"
 
255
#endif
 
256
 
 
257
int main(int argc, char *argv[])
 
258
{
 
259
#ifdef HAVE_MPI
 
260
  MPI_Init(&argc,&argv);
 
261
#endif
 
262
 
 
263
  printf("Please configure EpetraExt with:");
 
264
#if !defined(HAVE_PETSC)
 
265
  printf("--enable-petsc");
 
266
#endif
 
267
#if !defined(HAVE_ML_EPETRA)
 
268
  printf("--enable-epetra");
 
269
#endif
 
270
#if !defined(HAVE_ML_TEUCHOS)
 
271
  printf("--enable-teuchos");
 
272
#endif
 
273
#if !defined(HAVE_ML_AZTECOO)
 
274
  printf("--enable-aztecoo");
 
275
#endif
 
276
 
 
277
#ifdef HAVE_MPI
 
278
  MPI_Finalize();
 
279
#endif
 
280
  
 
281
  return(EXIT_SUCCESS);
 
282
}
 
283
#endif