Actual source code: chwirut1f.F90

  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 test_chwirut1.c
  9: !

 11: !
 12: ! ----------------------------------------------------------------------
 13: !
 14:       module chwirut1fmodule
 15:       use petsctao
 16: #include <petsc/finclude/petsctao.h>
 17:       PetscReal t(0:213)
 18:       PetscReal y(0:213)
 19:       PetscInt  m,n
 20:       end module chwirut1fmodule

 22:       program main
 23:       use chwirut1fmodule
 24: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 25: !                   Variable declarations
 26: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 27: !
 28: !  See additional variable declarations in the file chwirut1f.h

 30:       PetscErrorCode   ierr    ! used to check for functions returning nonzeros
 31:       Vec              x       ! solution vector
 32:       Vec              f       ! vector of functions
 33:       Tao              tao     ! Tao context
 34:       PetscInt         nhist
 35:       PetscMPIInt  size,rank    ! number of processes running
 36:       PetscReal      hist(100) ! objective value history
 37:       PetscReal      resid(100)! residual history
 38:       PetscReal      cnorm(100)! cnorm history
 39:       PetscInt      lits(100)   ! #ksp history
 40:       PetscInt oh
 41:       TaoConvergedReason reason

 43: !  Note: Any user-defined Fortran routines (such as FormGradient)
 44: !  MUST be declared as external.

 46:       external FormFunction

 48: !  Initialize TAO and PETSc
 49:       PetscCallA(PetscInitialize(ierr))

 51:       PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))
 52:       PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
 53:       PetscCheckA(size .eq. 1,PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only')

 55: !  Initialize problem parameters
 56:       m = 214
 57:       n = 3

 59: !  Allocate vectors for the solution and gradient
 60:       PetscCallA(VecCreateSeq(PETSC_COMM_SELF,n,x,ierr))
 61:       PetscCallA(VecCreateSeq(PETSC_COMM_SELF,m,f,ierr))

 63: !  The TAO code begins here

 65: !  Create TAO solver
 66:       PetscCallA(TaoCreate(PETSC_COMM_SELF,tao,ierr))
 67:       PetscCallA(TaoSetType(tao,TAOPOUNDERS,ierr))
 68: !  Set routines for function, gradient, and hessian evaluation

 70:       PetscCallA(TaoSetResidualRoutine(tao,f,FormFunction,0,ierr))

 72: !  Optional: Set initial guess
 73:       call InitializeData()
 74:       call FormStartingPoint(x)
 75:       PetscCallA(TaoSetSolution(tao, x, ierr))

 77: !  Check for TAO command line options
 78:       PetscCallA(TaoSetFromOptions(tao,ierr))
 79:       oh = 100
 80:       PetscCallA(TaoSetConvergenceHistory(tao,hist,resid,cnorm,lits,oh,PETSC_TRUE,ierr))
 81: !  SOLVE THE APPLICATION
 82:       PetscCallA(TaoSolve(tao,ierr))
 83:       PetscCallA(TaoGetConvergenceHistory(tao,nhist,ierr))
 84:       PetscCallA(TaoGetConvergedReason(tao, reason, ierr))
 85:       if (reason .le. 0) then
 86:          print *,'Tao failed.'
 87:          print *,'Try a different TAO method, adjust some parameters,'
 88:          print *,'or check the function evaluation routines.'
 89:       endif

 91: !  Free TAO data structures
 92:       PetscCallA(TaoDestroy(tao,ierr))

 94: !  Free PETSc data structures
 95:       PetscCallA(VecDestroy(x,ierr))
 96:       PetscCallA(VecDestroy(f,ierr))

 98:       PetscCallA(PetscFinalize(ierr))

100:       end

102: ! --------------------------------------------------------------------
103: !  FormFunction - Evaluates the function f(X) and gradient G(X)
104: !
105: !  Input Parameters:
106: !  tao - the Tao context
107: !  X   - input vector
108: !  dummy - not used
109: !
110: !  Output Parameters:
111: !  f - function vector

113:       subroutine FormFunction(tao, x, f, dummy, ierr)
114:       use chwirut1fmodule

116:       Tao        tao
117:       Vec              x,f
118:       PetscErrorCode   ierr
119:       PetscInt         dummy

121:       PetscInt         i
122:       PetscScalar, pointer, dimension(:)  :: x_v,f_v

124:       ierr = 0

126: !     Get pointers to vector data
127:       PetscCall(VecGetArrayF90(x,x_v,ierr))
128:       PetscCall(VecGetArrayF90(f,f_v,ierr))

130: !     Compute F(X)
131:       do i=0,m-1
132:          f_v(i+1) = y(i) - exp(-x_v(1)*t(i))/(x_v(2) + x_v(3)*t(i))
133:       enddo

135: !     Restore vectors
136:       PetscCall(VecRestoreArrayF90(X,x_v,ierr))
137:       PetscCall(VecRestoreArrayF90(F,f_v,ierr))

139:       end

141:       subroutine FormStartingPoint(x)
142:       use chwirut1fmodule

144:       Vec             x
145:       PetscScalar, pointer, dimension(:)  :: x_v
146:       PetscErrorCode  ierr

148:       PetscCall(VecGetArrayF90(x,x_v,ierr))
149:       x_v(1) = 0.15
150:       x_v(2) = 0.008
151:       x_v(3) = 0.01
152:       PetscCall(VecRestoreArrayF90(x,x_v,ierr))
153:       end

155:       subroutine InitializeData()
156:       use chwirut1fmodule

158:       integer i
159:       i=0
160:       y(i) =    92.9000;  t(i) =  0.5000; i=i+1
161:       y(i) =    78.7000;  t(i) =   0.6250; i=i+1
162:       y(i) =    64.2000;  t(i) =   0.7500; i=i+1
163:       y(i) =    64.9000;  t(i) =   0.8750; i=i+1
164:       y(i) =    57.1000;  t(i) =   1.0000; i=i+1
165:       y(i) =    43.3000;  t(i) =   1.2500; i=i+1
166:       y(i) =    31.1000;  t(i) =  1.7500; i=i+1
167:       y(i) =    23.6000;  t(i) =  2.2500; i=i+1
168:       y(i) =    31.0500;  t(i) =  1.7500; i=i+1
169:       y(i) =    23.7750;  t(i) =  2.2500; i=i+1
170:       y(i) =    17.7375;  t(i) =  2.7500; i=i+1
171:       y(i) =    13.8000;  t(i) =  3.2500; i=i+1
172:       y(i) =    11.5875;  t(i) =  3.7500; i=i+1
173:       y(i) =     9.4125;  t(i) =  4.2500; i=i+1
174:       y(i) =     7.7250;  t(i) =  4.7500; i=i+1
175:       y(i) =     7.3500;  t(i) =  5.2500; i=i+1
176:       y(i) =     8.0250;  t(i) =  5.7500; i=i+1
177:       y(i) =    90.6000;  t(i) =  0.5000; i=i+1
178:       y(i) =    76.9000;  t(i) =  0.6250; i=i+1
179:       y(i) =    71.6000;  t(i) = 0.7500; i=i+1
180:       y(i) =    63.6000;  t(i) =  0.8750; i=i+1
181:       y(i) =    54.0000;  t(i) =  1.0000; i=i+1
182:       y(i) =    39.2000;  t(i) =  1.2500; i=i+1
183:       y(i) =    29.3000;  t(i) = 1.7500; i=i+1
184:       y(i) =    21.4000;  t(i) =  2.2500; i=i+1
185:       y(i) =    29.1750;  t(i) =  1.7500; i=i+1
186:       y(i) =    22.1250;  t(i) =  2.2500; i=i+1
187:       y(i) =    17.5125;  t(i) =  2.7500; i=i+1
188:       y(i) =    14.2500;  t(i) =  3.2500; i=i+1
189:       y(i) =     9.4500;  t(i) =  3.7500; i=i+1
190:       y(i) =     9.1500;  t(i) =  4.2500; i=i+1
191:       y(i) =     7.9125;  t(i) =  4.7500; i=i+1
192:       y(i) =     8.4750;  t(i) =  5.2500; i=i+1
193:       y(i) =     6.1125;  t(i) =  5.7500; i=i+1
194:       y(i) =    80.0000;  t(i) =  0.5000; i=i+1
195:       y(i) =    79.0000;  t(i) =  0.6250; i=i+1
196:       y(i) =    63.8000;  t(i) =  0.7500; i=i+1
197:       y(i) =    57.2000;  t(i) =  0.8750; i=i+1
198:       y(i) =    53.2000;  t(i) =  1.0000; i=i+1
199:       y(i) =    42.5000;  t(i) =  1.2500; i=i+1
200:       y(i) =    26.8000;  t(i) =  1.7500; i=i+1
201:       y(i) =    20.4000;  t(i) =  2.2500; i=i+1
202:       y(i) =    26.8500;  t(i) =   1.7500; i=i+1
203:       y(i) =    21.0000;  t(i) =   2.2500; i=i+1
204:       y(i) =    16.4625;  t(i) =   2.7500; i=i+1
205:       y(i) =    12.5250;  t(i) =   3.2500; i=i+1
206:       y(i) =    10.5375;  t(i) =   3.7500; i=i+1
207:       y(i) =     8.5875;  t(i) =   4.2500; i=i+1
208:       y(i) =     7.1250;  t(i) =   4.7500; i=i+1
209:       y(i) =     6.1125;  t(i) =   5.2500; i=i+1
210:       y(i) =     5.9625;  t(i) =   5.7500; i=i+1
211:       y(i) =    74.1000;  t(i) =   0.5000; i=i+1
212:       y(i) =    67.3000;  t(i) =   0.6250; i=i+1
213:       y(i) =    60.8000;  t(i) =   0.7500; i=i+1
214:       y(i) =    55.5000;  t(i) =   0.8750; i=i+1
215:       y(i) =    50.3000;  t(i) =   1.0000; i=i+1
216:       y(i) =    41.0000;  t(i) =   1.2500; i=i+1
217:       y(i) =    29.4000;  t(i) =   1.7500; i=i+1
218:       y(i) =    20.4000;  t(i) =   2.2500; i=i+1
219:       y(i) =    29.3625;  t(i) =   1.7500; i=i+1
220:       y(i) =    21.1500;  t(i) =   2.2500; i=i+1
221:       y(i) =    16.7625;  t(i) =   2.7500; i=i+1
222:       y(i) =    13.2000;  t(i) =   3.2500; i=i+1
223:       y(i) =    10.8750;  t(i) =   3.7500; i=i+1
224:       y(i) =     8.1750;  t(i) =   4.2500; i=i+1
225:       y(i) =     7.3500;  t(i) =   4.7500; i=i+1
226:       y(i) =     5.9625;  t(i) =  5.2500; i=i+1
227:       y(i) =     5.6250;  t(i) =   5.7500; i=i+1
228:       y(i) =    81.5000;  t(i) =    .5000; i=i+1
229:       y(i) =    62.4000;  t(i) =    .7500; i=i+1
230:       y(i) =    32.5000;  t(i) =   1.5000; i=i+1
231:       y(i) =    12.4100;  t(i) =   3.0000; i=i+1
232:       y(i) =    13.1200;  t(i) =   3.0000; i=i+1
233:       y(i) =    15.5600;  t(i) =   3.0000; i=i+1
234:       y(i) =     5.6300;  t(i) =   6.0000; i=i+1
235:       y(i) =    78.0000;  t(i) =   .5000; i=i+1
236:       y(i) =    59.9000;  t(i) =    .7500; i=i+1
237:       y(i) =    33.2000;  t(i) =   1.5000; i=i+1
238:       y(i) =    13.8400;  t(i) =   3.0000; i=i+1
239:       y(i) =    12.7500;  t(i) =   3.0000; i=i+1
240:       y(i) =    14.6200;  t(i) =   3.0000; i=i+1
241:       y(i) =     3.9400;  t(i) =   6.0000; i=i+1
242:       y(i) =    76.8000;  t(i) =    .5000; i=i+1
243:       y(i) =    61.0000;  t(i) =    .7500; i=i+1
244:       y(i) =    32.9000;  t(i) =   1.5000; i=i+1
245:       y(i) =    13.8700;  t(i) = 3.0000; i=i+1
246:       y(i) =    11.8100;  t(i) =   3.0000; i=i+1
247:       y(i) =    13.3100;  t(i) =   3.0000; i=i+1
248:       y(i) =     5.4400;  t(i) =   6.0000; i=i+1
249:       y(i) =    78.0000;  t(i) =    .5000; i=i+1
250:       y(i) =    63.5000;  t(i) =    .7500; i=i+1
251:       y(i) =    33.8000;  t(i) =   1.5000; i=i+1
252:       y(i) =    12.5600;  t(i) =   3.0000; i=i+1
253:       y(i) =     5.6300;  t(i) =   6.0000; i=i+1
254:       y(i) =    12.7500;  t(i) =   3.0000; i=i+1
255:       y(i) =    13.1200;  t(i) =   3.0000; i=i+1
256:       y(i) =     5.4400;  t(i) =   6.0000; i=i+1
257:       y(i) =    76.8000;  t(i) =    .5000; i=i+1
258:       y(i) =    60.0000;  t(i) =    .7500; i=i+1
259:       y(i) =    47.8000;  t(i) =   1.0000; i=i+1
260:       y(i) =    32.0000;  t(i) =   1.5000; i=i+1
261:       y(i) =    22.2000;  t(i) =   2.0000; i=i+1
262:       y(i) =    22.5700;  t(i) =   2.0000; i=i+1
263:       y(i) =    18.8200;  t(i) =   2.5000; i=i+1
264:       y(i) =    13.9500;  t(i) =   3.0000; i=i+1
265:       y(i) =    11.2500;  t(i) =   4.0000; i=i+1
266:       y(i) =     9.0000;  t(i) =   5.0000; i=i+1
267:       y(i) =     6.6700;  t(i) =   6.0000; i=i+1
268:       y(i) =    75.8000;  t(i) =    .5000; i=i+1
269:       y(i) =    62.0000;  t(i) =    .7500; i=i+1
270:       y(i) =    48.8000;  t(i) =   1.0000; i=i+1
271:       y(i) =    35.2000;  t(i) =   1.5000; i=i+1
272:       y(i) =    20.0000;  t(i) =   2.0000; i=i+1
273:       y(i) =    20.3200;  t(i) =   2.0000; i=i+1
274:       y(i) =    19.3100;  t(i) =   2.5000; i=i+1
275:       y(i) =    12.7500;  t(i) =   3.0000; i=i+1
276:       y(i) =    10.4200;  t(i) =   4.0000; i=i+1
277:       y(i) =     7.3100;  t(i) =   5.0000; i=i+1
278:       y(i) =     7.4200;  t(i) =   6.0000; i=i+1
279:       y(i) =    70.5000;  t(i) =    .5000; i=i+1
280:       y(i) =    59.5000;  t(i) =    .7500; i=i+1
281:       y(i) =    48.5000;  t(i) =   1.0000; i=i+1
282:       y(i) =    35.8000;  t(i) =   1.5000; i=i+1
283:       y(i) =    21.0000;  t(i) =   2.0000; i=i+1
284:       y(i) =    21.6700;  t(i) =   2.0000; i=i+1
285:       y(i) =    21.0000;  t(i) =   2.5000; i=i+1
286:       y(i) =    15.6400;  t(i) =   3.0000; i=i+1
287:       y(i) =     8.1700;  t(i) =   4.0000; i=i+1
288:       y(i) =     8.5500;  t(i) =   5.0000; i=i+1
289:       y(i) =    10.1200;  t(i) =   6.0000; i=i+1
290:       y(i) =    78.0000;  t(i) =    .5000; i=i+1
291:       y(i) =    66.0000;  t(i) =    .6250; i=i+1
292:       y(i) =    62.0000;  t(i) =    .7500; i=i+1
293:       y(i) =    58.0000;  t(i) =    .8750; i=i+1
294:       y(i) =    47.7000;  t(i) =   1.0000; i=i+1
295:       y(i) =    37.8000;  t(i) =   1.2500; i=i+1
296:       y(i) =    20.2000;  t(i) =   2.2500; i=i+1
297:       y(i) =    21.0700;  t(i) =   2.2500; i=i+1
298:       y(i) =    13.8700;  t(i) =   2.7500; i=i+1
299:       y(i) =     9.6700;  t(i) =   3.2500; i=i+1
300:       y(i) =     7.7600;  t(i) =   3.7500; i=i+1
301:       y(i) =     5.4400;  t(i) =  4.2500; i=i+1
302:       y(i) =     4.8700;  t(i) =  4.7500; i=i+1
303:       y(i) =     4.0100;  t(i) =   5.2500; i=i+1
304:       y(i) =     3.7500;  t(i) =   5.7500; i=i+1
305:       y(i) =    24.1900;  t(i) =   3.0000; i=i+1
306:       y(i) =    25.7600;  t(i) =   3.0000; i=i+1
307:       y(i) =    18.0700;  t(i) =   3.0000; i=i+1
308:       y(i) =    11.8100;  t(i) =   3.0000; i=i+1
309:       y(i) =    12.0700;  t(i) =   3.0000; i=i+1
310:       y(i) =    16.1200;  t(i) =   3.0000; i=i+1
311:       y(i) =    70.8000;  t(i) =    .5000; i=i+1
312:       y(i) =    54.7000;  t(i) =    .7500; i=i+1
313:       y(i) =    48.0000;  t(i) =   1.0000; i=i+1
314:       y(i) =    39.8000;  t(i) =   1.5000; i=i+1
315:       y(i) =    29.8000;  t(i) =   2.0000; i=i+1
316:       y(i) =    23.7000;  t(i) =   2.5000; i=i+1
317:       y(i) =    29.6200;  t(i) =   2.0000; i=i+1
318:       y(i) =    23.8100;  t(i) =   2.5000; i=i+1
319:       y(i) =    17.7000;  t(i) =   3.0000; i=i+1
320:       y(i) =    11.5500;  t(i) =   4.0000; i=i+1
321:       y(i) =    12.0700;  t(i) =   5.0000; i=i+1
322:       y(i) =     8.7400;  t(i) =   6.0000; i=i+1
323:       y(i) =    80.7000;  t(i) =    .5000; i=i+1
324:       y(i) =    61.3000;  t(i) =    .7500; i=i+1
325:       y(i) =    47.5000;  t(i) =   1.0000; i=i+1
326:       y(i) =    29.0000;  t(i) =   1.5000; i=i+1
327:       y(i) =    24.0000;  t(i) =   2.0000; i=i+1
328:       y(i) =    17.7000;  t(i) =   2.5000; i=i+1
329:       y(i) =    24.5600;  t(i) =   2.0000; i=i+1
330:       y(i) =    18.6700;  t(i) =   2.5000; i=i+1
331:       y(i) =    16.2400;  t(i) =   3.0000; i=i+1
332:       y(i) =     8.7400;  t(i) =   4.0000; i=i+1
333:       y(i) =     7.8700;  t(i) =   5.0000; i=i+1
334:       y(i) =     8.5100;  t(i) =   6.0000; i=i+1
335:       y(i) =    66.7000;  t(i) =    .5000; i=i+1
336:       y(i) =    59.2000;  t(i) =    .7500; i=i+1
337:       y(i) =    40.8000;  t(i) =   1.0000; i=i+1
338:       y(i) =    30.7000;  t(i) =   1.5000; i=i+1
339:       y(i) =    25.7000;  t(i) =   2.0000; i=i+1
340:       y(i) =    16.3000;  t(i) =   2.5000; i=i+1
341:       y(i) =    25.9900;  t(i) =   2.0000; i=i+1
342:       y(i) =    16.9500;  t(i) =   2.5000; i=i+1
343:       y(i) =    13.3500;  t(i) =   3.0000; i=i+1
344:       y(i) =     8.6200;  t(i) =   4.0000; i=i+1
345:       y(i) =     7.2000;  t(i) =   5.0000; i=i+1
346:       y(i) =     6.6400;  t(i) =   6.0000; i=i+1
347:       y(i) =    13.6900;  t(i) =   3.0000; i=i+1
348:       y(i) =    81.0000;  t(i) =    .5000; i=i+1
349:       y(i) =    64.5000;  t(i) =    .7500; i=i+1
350:       y(i) =    35.5000;  t(i) =   1.5000; i=i+1
351:       y(i) =    13.3100;  t(i) =   3.0000; i=i+1
352:       y(i) =     4.8700;  t(i) =   6.0000; i=i+1
353:       y(i) =    12.9400;  t(i) =   3.0000; i=i+1
354:       y(i) =     5.0600;  t(i) =   6.0000; i=i+1
355:       y(i) =    15.1900;  t(i) =   3.0000; i=i+1
356:       y(i) =    14.6200;  t(i) =   3.0000; i=i+1
357:       y(i) =    15.6400;  t(i) =   3.0000; i=i+1
358:       y(i) =    25.5000;  t(i) =   1.7500; i=i+1
359:       y(i) =    25.9500;  t(i) =   1.7500; i=i+1
360:       y(i) =    81.7000;  t(i) =    .5000; i=i+1
361:       y(i) =    61.6000;  t(i) =    .7500; i=i+1
362:       y(i) =    29.8000;  t(i) =   1.7500; i=i+1
363:       y(i) =    29.8100;  t(i) =   1.7500; i=i+1
364:       y(i) =    17.1700;  t(i) =   2.7500; i=i+1
365:       y(i) =    10.3900;  t(i) =   3.7500; i=i+1
366:       y(i) =    28.4000;  t(i) =   1.7500; i=i+1
367:       y(i) =    28.6900;  t(i) =   1.7500; i=i+1
368:       y(i) =    81.3000;  t(i) =    .5000; i=i+1
369:       y(i) =    60.9000;  t(i) =    .7500; i=i+1
370:       y(i) =    16.6500;  t(i) =   2.7500; i=i+1
371:       y(i) =    10.0500;  t(i) =   3.7500; i=i+1
372:       y(i) =    28.9000;  t(i) =   1.7500; i=i+1
373:       y(i) =    28.9500;  t(i) =   1.7500; i=i+1

375:       end

377: !/*TEST
378: !
379: !   build:
380: !      requires: !complex
381: !
382: !   test:
383: !      args: -tao_monitor_short -tao_max_it 100 -tao_type pounders -tao_gatol 1.e-5
384: !      requires: !single
385: !
386: !TEST*/