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