xref: /petsc/src/tao/leastsquares/tutorials/chwirut2f.F90 (revision ccb4e88a40f0b86eaeca07ff64c64e4de2fae686)
1!  Program usage: mpiexec -n 1 chwirut1f [-help] [all TAO options]
2!
3!  Description:  This example demonstrates use of the TAO package to solve a
4!  nonlinear least-squares problem on a single processor.  We minimize the
5!  Chwirut function:
6!       sum_{i=0}^{n/2-1} ( alpha*(x_{2i+1}-x_{2i}^2)^2 + (1-x_{2i})^2)
7!
8!  The C version of this code is chwirut1.c
9!
10!!/*T
11!  Concepts: TAO^Solving an unconstrained minimization problem
12!  Routines: TaoCreate();
13!  Routines: TaoSetType();
14!  Routines: TaoSetResidualRoutine();
15!  Routines: TaoSetInitialVector();
16!  Routines: TaoSetFromOptions();
17!  Routines: TaoSolve();
18!  Routines: TaoDestroy();
19!  Processors: n
20!T*/
21
22!
23! ----------------------------------------------------------------------
24!
25#include "chwirut2f.h"
26
27! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
28!                   Variable declarations
29! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
30!
31!  See additional variable declarations in the file chwirut2f.h
32
33      PetscErrorCode   ierr    ! used to check for functions returning nonzeros
34      Vec              x       ! solution vector
35      Vec              f       ! vector of functions
36      Tao        tao     ! Tao context
37
38!  Note: Any user-defined Fortran routines (such as FormGradient)
39!  MUST be declared as external.
40
41      external FormFunction
42
43!  Initialize TAO and PETSc
44      call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
45      if (ierr .ne. 0) then
46         print*,'Unable to initialize PETSc'
47         stop
48      endif
49
50      call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
51      CHKERRA(ierr)
52      call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
53      CHKERRA(ierr)
54
55!  Initialize problem parameters
56      call InitializeData()
57
58      if (rank .eq. 0) then
59!  Allocate vectors for the solution and gradient
60         call VecCreateSeq(PETSC_COMM_SELF,n,x,ierr)
61         CHKERRA(ierr)
62         call VecCreateSeq(PETSC_COMM_SELF,m,f,ierr)
63         CHKERRA(ierr)
64
65!     The TAO code begins here
66
67!     Create TAO solver
68         call TaoCreate(PETSC_COMM_SELF,tao,ierr)
69         CHKERRA(ierr)
70         call TaoSetType(tao,TAOPOUNDERS,ierr)
71         CHKERRA(ierr)
72
73!     Set routines for function, gradient, and hessian evaluation
74         call TaoSetResidualRoutine(tao,f,                    &
75     &        FormFunction,0,ierr)
76         CHKERRA(ierr)
77
78!     Optional: Set initial guess
79         call FormStartingPoint(x)
80         call TaoSetInitialVector(tao, x, ierr)
81         CHKERRA(ierr)
82
83!     Check for TAO command line options
84         call TaoSetFromOptions(tao,ierr)
85         CHKERRA(ierr)
86!     SOLVE THE APPLICATION
87         call TaoSolve(tao,ierr)
88         CHKERRA(ierr)
89
90!     Free TAO data structures
91         call TaoDestroy(tao,ierr)
92         CHKERRA(ierr)
93
94!     Free PETSc data structures
95         call VecDestroy(x,ierr)
96         CHKERRA(ierr)
97         call VecDestroy(f,ierr)
98         CHKERRA(ierr)
99         call StopWorkers(ierr)
100         CHKERRA(ierr)
101
102      else
103         call TaskWorker(ierr)
104         CHKERRA(ierr)
105      endif
106
107      call PetscFinalize(ierr)
108      end
109
110! --------------------------------------------------------------------
111!  FormFunction - Evaluates the function f(X) and gradient G(X)
112!
113!  Input Parameters:
114!  tao - the Tao context
115!  X   - input vector
116!  dummy - not used
117!
118!  Output Parameters:
119!  f - function vector
120
121      subroutine FormFunction(tao, x, f, dummy, ierr)
122#include "chwirut2f.h"
123
124      Tao        tao
125      Vec              x,f
126      PetscErrorCode   ierr
127
128      PetscInt         i,checkedin
129      PetscInt         finished_tasks
130      PetscMPIInt      next_task
131      PetscMPIInt      status(MPI_STATUS_SIZE),tag,source
132      PetscInt         dummy
133
134! PETSc's VecGetArray acts differently in Fortran than it does in C.
135! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
136! will return an array of doubles referenced by x_array offset by x_index.
137!  i.e.,  to reference the kth element of X, use x_array(k + x_index).
138! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
139      PetscReal        f_v(0:1),x_v(0:1),fval(1)
140      PetscOffset      f_i,x_i
141
142      ierr = 0
143
144!     Get pointers to vector data
145      call VecGetArray(x,x_v,x_i,ierr)
146      CHKERRQ(ierr)
147      call VecGetArray(f,f_v,f_i,ierr)
148      CHKERRQ(ierr)
149
150!     Compute F(X)
151      if (size .eq. 1) then
152         ! Single processor
153         do i=0,m-1
154            call RunSimulation(x_v(x_i),i,f_v(i+f_i),ierr)
155         enddo
156      else
157         ! Multiprocessor main
158         next_task = zero
159         finished_tasks = 0
160         checkedin = 0
161
162         do while (finished_tasks .lt. m .or. checkedin .lt. size-1)
163            call MPI_Recv(fval,one,MPIU_SCALAR,MPI_ANY_SOURCE,               &
164     &           MPI_ANY_TAG,PETSC_COMM_WORLD,status,ierr)
165            tag = status(MPI_TAG)
166            source = status(MPI_SOURCE)
167            if (tag .eq. IDLE_TAG) then
168               checkedin = checkedin + 1
169            else
170               f_v(f_i+tag) = fval(1)
171               finished_tasks = finished_tasks + 1
172            endif
173            if (next_task .lt. m) then
174               ! Send task to worker
175               call MPI_Send(x_v(x_i),nn,MPIU_SCALAR,source,next_task,             &
176     &              PETSC_COMM_WORLD,ierr)
177               next_task = next_task + one
178            else
179               ! Send idle message to worker
180               call MPI_Send(x_v(x_i),nn,MPIU_SCALAR,source,IDLE_TAG,              &
181     &              PETSC_COMM_WORLD,ierr)
182            end if
183         enddo
184      endif
185
186!     Restore vectors
187      call VecRestoreArray(x,x_v,x_i,ierr)
188      CHKERRQ(ierr)
189      call VecRestoreArray(F,f_v,f_i,ierr)
190      CHKERRQ(ierr)
191      return
192      end
193
194      subroutine FormStartingPoint(x)
195#include "chwirut2f.h"
196
197      Vec             x
198      PetscReal       x_v(0:1)
199      PetscOffset     x_i
200      PetscErrorCode  ierr
201
202      call VecGetArray(x,x_v,x_i,ierr)
203      CHKERRQ(ierr)
204      x_v(x_i) = 0.15
205      x_v(x_i+1) = 0.008
206      x_v(x_i+2) = 0.01
207      call VecRestoreArray(x,x_v,x_i,ierr)
208      CHKERRQ(ierr)
209      return
210      end
211
212      subroutine InitializeData()
213#include "chwirut2f.h"
214
215      PetscInt i
216      i=0
217      y(i) =    92.9000;  t(i) =  0.5000; i=i+1
218      y(i) =    78.7000;  t(i) =   0.6250; i=i+1
219      y(i) =    64.2000;  t(i) =   0.7500; i=i+1
220      y(i) =    64.9000;  t(i) =   0.8750; i=i+1
221      y(i) =    57.1000;  t(i) =   1.0000; i=i+1
222      y(i) =    43.3000;  t(i) =   1.2500; i=i+1
223      y(i) =    31.1000;  t(i) =  1.7500; i=i+1
224      y(i) =    23.6000;  t(i) =  2.2500; i=i+1
225      y(i) =    31.0500;  t(i) =  1.7500; i=i+1
226      y(i) =    23.7750;  t(i) =  2.2500; i=i+1
227      y(i) =    17.7375;  t(i) =  2.7500; i=i+1
228      y(i) =    13.8000;  t(i) =  3.2500; i=i+1
229      y(i) =    11.5875;  t(i) =  3.7500; i=i+1
230      y(i) =     9.4125;  t(i) =  4.2500; i=i+1
231      y(i) =     7.7250;  t(i) =  4.7500; i=i+1
232      y(i) =     7.3500;  t(i) =  5.2500; i=i+1
233      y(i) =     8.0250;  t(i) =  5.7500; i=i+1
234      y(i) =    90.6000;  t(i) =  0.5000; i=i+1
235      y(i) =    76.9000;  t(i) =  0.6250; i=i+1
236      y(i) =    71.6000;  t(i) = 0.7500; i=i+1
237      y(i) =    63.6000;  t(i) =  0.8750; i=i+1
238      y(i) =    54.0000;  t(i) =  1.0000; i=i+1
239      y(i) =    39.2000;  t(i) =  1.2500; i=i+1
240      y(i) =    29.3000;  t(i) = 1.7500; i=i+1
241      y(i) =    21.4000;  t(i) =  2.2500; i=i+1
242      y(i) =    29.1750;  t(i) =  1.7500; i=i+1
243      y(i) =    22.1250;  t(i) =  2.2500; i=i+1
244      y(i) =    17.5125;  t(i) =  2.7500; i=i+1
245      y(i) =    14.2500;  t(i) =  3.2500; i=i+1
246      y(i) =     9.4500;  t(i) =  3.7500; i=i+1
247      y(i) =     9.1500;  t(i) =  4.2500; i=i+1
248      y(i) =     7.9125;  t(i) =  4.7500; i=i+1
249      y(i) =     8.4750;  t(i) =  5.2500; i=i+1
250      y(i) =     6.1125;  t(i) =  5.7500; i=i+1
251      y(i) =    80.0000;  t(i) =  0.5000; i=i+1
252      y(i) =    79.0000;  t(i) =  0.6250; i=i+1
253      y(i) =    63.8000;  t(i) =  0.7500; i=i+1
254      y(i) =    57.2000;  t(i) =  0.8750; i=i+1
255      y(i) =    53.2000;  t(i) =  1.0000; i=i+1
256      y(i) =    42.5000;  t(i) =  1.2500; i=i+1
257      y(i) =    26.8000;  t(i) =  1.7500; i=i+1
258      y(i) =    20.4000;  t(i) =  2.2500; i=i+1
259      y(i) =    26.8500;  t(i) =   1.7500; i=i+1
260      y(i) =    21.0000;  t(i) =   2.2500; i=i+1
261      y(i) =    16.4625;  t(i) =   2.7500; i=i+1
262      y(i) =    12.5250;  t(i) =   3.2500; i=i+1
263      y(i) =    10.5375;  t(i) =   3.7500; i=i+1
264      y(i) =     8.5875;  t(i) =   4.2500; i=i+1
265      y(i) =     7.1250;  t(i) =   4.7500; i=i+1
266      y(i) =     6.1125;  t(i) =   5.2500; i=i+1
267      y(i) =     5.9625;  t(i) =   5.7500; i=i+1
268      y(i) =    74.1000;  t(i) =   0.5000; i=i+1
269      y(i) =    67.3000;  t(i) =   0.6250; i=i+1
270      y(i) =    60.8000;  t(i) =   0.7500; i=i+1
271      y(i) =    55.5000;  t(i) =   0.8750; i=i+1
272      y(i) =    50.3000;  t(i) =   1.0000; i=i+1
273      y(i) =    41.0000;  t(i) =   1.2500; i=i+1
274      y(i) =    29.4000;  t(i) =   1.7500; i=i+1
275      y(i) =    20.4000;  t(i) =   2.2500; i=i+1
276      y(i) =    29.3625;  t(i) =   1.7500; i=i+1
277      y(i) =    21.1500;  t(i) =   2.2500; i=i+1
278      y(i) =    16.7625;  t(i) =   2.7500; i=i+1
279      y(i) =    13.2000;  t(i) =   3.2500; i=i+1
280      y(i) =    10.8750;  t(i) =   3.7500; i=i+1
281      y(i) =     8.1750;  t(i) =   4.2500; i=i+1
282      y(i) =     7.3500;  t(i) =   4.7500; i=i+1
283      y(i) =     5.9625;  t(i) =  5.2500; i=i+1
284      y(i) =     5.6250;  t(i) =   5.7500; i=i+1
285      y(i) =    81.5000;  t(i) =    .5000; i=i+1
286      y(i) =    62.4000;  t(i) =    .7500; i=i+1
287      y(i) =    32.5000;  t(i) =   1.5000; i=i+1
288      y(i) =    12.4100;  t(i) =   3.0000; i=i+1
289      y(i) =    13.1200;  t(i) =   3.0000; i=i+1
290      y(i) =    15.5600;  t(i) =   3.0000; i=i+1
291      y(i) =     5.6300;  t(i) =   6.0000; i=i+1
292      y(i) =    78.0000;  t(i) =   .5000; i=i+1
293      y(i) =    59.9000;  t(i) =    .7500; i=i+1
294      y(i) =    33.2000;  t(i) =   1.5000; i=i+1
295      y(i) =    13.8400;  t(i) =   3.0000; i=i+1
296      y(i) =    12.7500;  t(i) =   3.0000; i=i+1
297      y(i) =    14.6200;  t(i) =   3.0000; i=i+1
298      y(i) =     3.9400;  t(i) =   6.0000; i=i+1
299      y(i) =    76.8000;  t(i) =    .5000; i=i+1
300      y(i) =    61.0000;  t(i) =    .7500; i=i+1
301      y(i) =    32.9000;  t(i) =   1.5000; i=i+1
302      y(i) =    13.8700;  t(i) = 3.0000; i=i+1
303      y(i) =    11.8100;  t(i) =   3.0000; i=i+1
304      y(i) =    13.3100;  t(i) =   3.0000; i=i+1
305      y(i) =     5.4400;  t(i) =   6.0000; i=i+1
306      y(i) =    78.0000;  t(i) =    .5000; i=i+1
307      y(i) =    63.5000;  t(i) =    .7500; i=i+1
308      y(i) =    33.8000;  t(i) =   1.5000; i=i+1
309      y(i) =    12.5600;  t(i) =   3.0000; i=i+1
310      y(i) =     5.6300;  t(i) =   6.0000; i=i+1
311      y(i) =    12.7500;  t(i) =   3.0000; i=i+1
312      y(i) =    13.1200;  t(i) =   3.0000; i=i+1
313      y(i) =     5.4400;  t(i) =   6.0000; i=i+1
314      y(i) =    76.8000;  t(i) =    .5000; i=i+1
315      y(i) =    60.0000;  t(i) =    .7500; i=i+1
316      y(i) =    47.8000;  t(i) =   1.0000; i=i+1
317      y(i) =    32.0000;  t(i) =   1.5000; i=i+1
318      y(i) =    22.2000;  t(i) =   2.0000; i=i+1
319      y(i) =    22.5700;  t(i) =   2.0000; i=i+1
320      y(i) =    18.8200;  t(i) =   2.5000; i=i+1
321      y(i) =    13.9500;  t(i) =   3.0000; i=i+1
322      y(i) =    11.2500;  t(i) =   4.0000; i=i+1
323      y(i) =     9.0000;  t(i) =   5.0000; i=i+1
324      y(i) =     6.6700;  t(i) =   6.0000; i=i+1
325      y(i) =    75.8000;  t(i) =    .5000; i=i+1
326      y(i) =    62.0000;  t(i) =    .7500; i=i+1
327      y(i) =    48.8000;  t(i) =   1.0000; i=i+1
328      y(i) =    35.2000;  t(i) =   1.5000; i=i+1
329      y(i) =    20.0000;  t(i) =   2.0000; i=i+1
330      y(i) =    20.3200;  t(i) =   2.0000; i=i+1
331      y(i) =    19.3100;  t(i) =   2.5000; i=i+1
332      y(i) =    12.7500;  t(i) =   3.0000; i=i+1
333      y(i) =    10.4200;  t(i) =   4.0000; i=i+1
334      y(i) =     7.3100;  t(i) =   5.0000; i=i+1
335      y(i) =     7.4200;  t(i) =   6.0000; i=i+1
336      y(i) =    70.5000;  t(i) =    .5000; i=i+1
337      y(i) =    59.5000;  t(i) =    .7500; i=i+1
338      y(i) =    48.5000;  t(i) =   1.0000; i=i+1
339      y(i) =    35.8000;  t(i) =   1.5000; i=i+1
340      y(i) =    21.0000;  t(i) =   2.0000; i=i+1
341      y(i) =    21.6700;  t(i) =   2.0000; i=i+1
342      y(i) =    21.0000;  t(i) =   2.5000; i=i+1
343      y(i) =    15.6400;  t(i) =   3.0000; i=i+1
344      y(i) =     8.1700;  t(i) =   4.0000; i=i+1
345      y(i) =     8.5500;  t(i) =   5.0000; i=i+1
346      y(i) =    10.1200;  t(i) =   6.0000; i=i+1
347      y(i) =    78.0000;  t(i) =    .5000; i=i+1
348      y(i) =    66.0000;  t(i) =    .6250; i=i+1
349      y(i) =    62.0000;  t(i) =    .7500; i=i+1
350      y(i) =    58.0000;  t(i) =    .8750; i=i+1
351      y(i) =    47.7000;  t(i) =   1.0000; i=i+1
352      y(i) =    37.8000;  t(i) =   1.2500; i=i+1
353      y(i) =    20.2000;  t(i) =   2.2500; i=i+1
354      y(i) =    21.0700;  t(i) =   2.2500; i=i+1
355      y(i) =    13.8700;  t(i) =   2.7500; i=i+1
356      y(i) =     9.6700;  t(i) =   3.2500; i=i+1
357      y(i) =     7.7600;  t(i) =   3.7500; i=i+1
358      y(i) =     5.4400;  t(i) =  4.2500; i=i+1
359      y(i) =     4.8700;  t(i) =  4.7500; i=i+1
360      y(i) =     4.0100;  t(i) =   5.2500; i=i+1
361      y(i) =     3.7500;  t(i) =   5.7500; i=i+1
362      y(i) =    24.1900;  t(i) =   3.0000; i=i+1
363      y(i) =    25.7600;  t(i) =   3.0000; i=i+1
364      y(i) =    18.0700;  t(i) =   3.0000; i=i+1
365      y(i) =    11.8100;  t(i) =   3.0000; i=i+1
366      y(i) =    12.0700;  t(i) =   3.0000; i=i+1
367      y(i) =    16.1200;  t(i) =   3.0000; i=i+1
368      y(i) =    70.8000;  t(i) =    .5000; i=i+1
369      y(i) =    54.7000;  t(i) =    .7500; i=i+1
370      y(i) =    48.0000;  t(i) =   1.0000; i=i+1
371      y(i) =    39.8000;  t(i) =   1.5000; i=i+1
372      y(i) =    29.8000;  t(i) =   2.0000; i=i+1
373      y(i) =    23.7000;  t(i) =   2.5000; i=i+1
374      y(i) =    29.6200;  t(i) =   2.0000; i=i+1
375      y(i) =    23.8100;  t(i) =   2.5000; i=i+1
376      y(i) =    17.7000;  t(i) =   3.0000; i=i+1
377      y(i) =    11.5500;  t(i) =   4.0000; i=i+1
378      y(i) =    12.0700;  t(i) =   5.0000; i=i+1
379      y(i) =     8.7400;  t(i) =   6.0000; i=i+1
380      y(i) =    80.7000;  t(i) =    .5000; i=i+1
381      y(i) =    61.3000;  t(i) =    .7500; i=i+1
382      y(i) =    47.5000;  t(i) =   1.0000; i=i+1
383      y(i) =    29.0000;  t(i) =   1.5000; i=i+1
384      y(i) =    24.0000;  t(i) =   2.0000; i=i+1
385      y(i) =    17.7000;  t(i) =   2.5000; i=i+1
386      y(i) =    24.5600;  t(i) =   2.0000; i=i+1
387      y(i) =    18.6700;  t(i) =   2.5000; i=i+1
388      y(i) =    16.2400;  t(i) =   3.0000; i=i+1
389      y(i) =     8.7400;  t(i) =   4.0000; i=i+1
390      y(i) =     7.8700;  t(i) =   5.0000; i=i+1
391      y(i) =     8.5100;  t(i) =   6.0000; i=i+1
392      y(i) =    66.7000;  t(i) =    .5000; i=i+1
393      y(i) =    59.2000;  t(i) =    .7500; i=i+1
394      y(i) =    40.8000;  t(i) =   1.0000; i=i+1
395      y(i) =    30.7000;  t(i) =   1.5000; i=i+1
396      y(i) =    25.7000;  t(i) =   2.0000; i=i+1
397      y(i) =    16.3000;  t(i) =   2.5000; i=i+1
398      y(i) =    25.9900;  t(i) =   2.0000; i=i+1
399      y(i) =    16.9500;  t(i) =   2.5000; i=i+1
400      y(i) =    13.3500;  t(i) =   3.0000; i=i+1
401      y(i) =     8.6200;  t(i) =   4.0000; i=i+1
402      y(i) =     7.2000;  t(i) =   5.0000; i=i+1
403      y(i) =     6.6400;  t(i) =   6.0000; i=i+1
404      y(i) =    13.6900;  t(i) =   3.0000; i=i+1
405      y(i) =    81.0000;  t(i) =    .5000; i=i+1
406      y(i) =    64.5000;  t(i) =    .7500; i=i+1
407      y(i) =    35.5000;  t(i) =   1.5000; i=i+1
408      y(i) =    13.3100;  t(i) =   3.0000; i=i+1
409      y(i) =     4.8700;  t(i) =   6.0000; i=i+1
410      y(i) =    12.9400;  t(i) =   3.0000; i=i+1
411      y(i) =     5.0600;  t(i) =   6.0000; i=i+1
412      y(i) =    15.1900;  t(i) =   3.0000; i=i+1
413      y(i) =    14.6200;  t(i) =   3.0000; i=i+1
414      y(i) =    15.6400;  t(i) =   3.0000; i=i+1
415      y(i) =    25.5000;  t(i) =   1.7500; i=i+1
416      y(i) =    25.9500;  t(i) =   1.7500; i=i+1
417      y(i) =    81.7000;  t(i) =    .5000; i=i+1
418      y(i) =    61.6000;  t(i) =    .7500; i=i+1
419      y(i) =    29.8000;  t(i) =   1.7500; i=i+1
420      y(i) =    29.8100;  t(i) =   1.7500; i=i+1
421      y(i) =    17.1700;  t(i) =   2.7500; i=i+1
422      y(i) =    10.3900;  t(i) =   3.7500; i=i+1
423      y(i) =    28.4000;  t(i) =   1.7500; i=i+1
424      y(i) =    28.6900;  t(i) =   1.7500; i=i+1
425      y(i) =    81.3000;  t(i) =    .5000; i=i+1
426      y(i) =    60.9000;  t(i) =    .7500; i=i+1
427      y(i) =    16.6500;  t(i) =   2.7500; i=i+1
428      y(i) =    10.0500;  t(i) =   3.7500; i=i+1
429      y(i) =    28.9000;  t(i) =   1.7500; i=i+1
430      y(i) =    28.9500;  t(i) =   1.7500; i=i+1
431
432      return
433      end
434
435      subroutine TaskWorker(ierr)
436#include "chwirut2f.h"
437
438      PetscErrorCode ierr
439      PetscReal x(n),f(1)
440      PetscMPIInt tag
441      PetscInt index
442      PetscMPIInt status(MPI_STATUS_SIZE)
443
444      tag = IDLE_TAG
445      f   = 0.0
446      ! Send check-in message to rank-0
447      call MPI_Send(f,one,MPIU_SCALAR,zero,IDLE_TAG,PETSC_COMM_WORLD,ierr)
448      CHKERRQ(ierr)
449      do while (tag .ne. DIE_TAG)
450         call MPI_Recv(x,nn,MPIU_SCALAR,zero,MPI_ANY_TAG,PETSC_COMM_WORLD,     &
451     &        status,ierr)
452         CHKERRQ(ierr)
453         tag = status(MPI_TAG)
454         if (tag .eq. IDLE_TAG) then
455            call MPI_Send(f,one,MPIU_SCALAR,zero,IDLE_TAG,PETSC_COMM_WORLD,     &
456     &           ierr)
457            CHKERRQ(ierr)
458         else if (tag .ne. DIE_TAG) then
459            index = tag
460            ! Compute local part of residual
461            call RunSimulation(x,index,f(1),ierr)
462            CHKERRQ(ierr)
463
464            ! Return residual to rank-0
465            call MPI_Send(f,one,MPIU_SCALAR,zero,tag,PETSC_COMM_WORLD,ierr)
466            CHKERRQ(ierr)
467         end if
468      enddo
469      ierr = 0
470      return
471      end
472
473      subroutine RunSimulation(x,i,f,ierr)
474#include "chwirut2f.h"
475
476      PetscReal x(n),f
477      PetscInt i
478      PetscErrorCode ierr
479      f = y(i) - exp(-x(1)*t(i))/(x(2)+x(3)*t(i))
480      ierr = 0
481      return
482      end
483
484      subroutine StopWorkers(ierr)
485#include "chwirut2f.h"
486
487      integer checkedin
488      PetscMPIInt status(MPI_STATUS_SIZE)
489      PetscMPIInt source
490      PetscReal f(1),x(n)
491      PetscErrorCode ierr
492      PetscInt i
493
494      checkedin=0
495      do while (checkedin .lt. size-1)
496         call MPI_Recv(f,one,MPIU_SCALAR,MPI_ANY_SOURCE,MPI_ANY_TAG,         &
497     &        PETSC_COMM_WORLD,status,ierr)
498         CHKERRQ(ierr)
499         checkedin=checkedin+1
500         source = status(MPI_SOURCE)
501         do i=1,n
502           x(i) = 0.0
503         enddo
504         call MPI_Send(x,nn,MPIU_SCALAR,source,DIE_TAG,PETSC_COMM_WORLD,    &
505     &        ierr)
506         CHKERRQ(ierr)
507      enddo
508      ierr = 0
509      return
510      end
511
512!/*TEST
513!
514!   build:
515!      requires: !complex
516!
517!   test:
518!      nsize: 3
519!      args: -tao_smonitor -tao_max_it 100 -tao_type pounders -tao_gatol 1.e-5
520!      requires: !single
521!
522!
523!TEST*/
524