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: #include <petsc/finclude/petsctao.h>
 15: module chwirut1fmodule
 16:   use petsctao
 17:   implicit none

 19:   PetscReal t(0:213)
 20:   PetscReal y(0:213)
 21:   PetscInt m, n

 23: contains

 25: ! --------------------------------------------------------------------
 26: !  FormFunction - Evaluates the function f(X) and gradient G(X)
 27: !
 28: !  Input Parameters:
 29: !  tao - the Tao context
 30: !  X   - input vector
 31: !  dummy - not used
 32: !
 33: !  Output Parameters:
 34: !  f - function vector
 35:   subroutine FormFunction(ta, x, f, dummy, ierr)

 37:     Tao ta
 38:     Vec x, f
 39:     PetscErrorCode ierr
 40:     PetscInt dummy

 42:     PetscInt i
 43:     PetscScalar, pointer, dimension(:)  :: x_v, f_v

 45:     ierr = 0

 47: !     Get pointers to vector data
 48:     PetscCall(VecGetArray(x, x_v, ierr))
 49:     PetscCall(VecGetArray(f, f_v, ierr))

 51: !     Compute F(X)
 52:     do i = 0, m - 1
 53:       f_v(i + 1) = y(i) - exp(-x_v(1)*t(i))/(x_v(2) + x_v(3)*t(i))
 54:     end do

 56: !     Restore vectors
 57:     PetscCall(VecRestoreArray(X, x_v, ierr))
 58:     PetscCall(VecRestoreArray(F, f_v, ierr))

 60:   end

 62:   subroutine FormStartingPoint(x)
 63:     Vec x
 64:     PetscScalar, pointer, dimension(:)  :: x_v
 65:     PetscErrorCode ierr

 67:     PetscCall(VecGetArray(x, x_v, ierr))
 68:     x_v(1) = 0.15
 69:     x_v(2) = 0.008
 70:     x_v(3) = 0.01
 71:     PetscCall(VecRestoreArray(x, x_v, ierr))
 72:   end

 74:   subroutine InitializeData()
 75:     integer i
 76:     i = 0
 77:     y(i) = 92.9000; t(i) = 0.5000; i = i + 1
 78:     y(i) = 78.7000; t(i) = 0.6250; i = i + 1
 79:     y(i) = 64.2000; t(i) = 0.7500; i = i + 1
 80:     y(i) = 64.9000; t(i) = 0.8750; i = i + 1
 81:     y(i) = 57.1000; t(i) = 1.0000; i = i + 1
 82:     y(i) = 43.3000; t(i) = 1.2500; i = i + 1
 83:     y(i) = 31.1000; t(i) = 1.7500; i = i + 1
 84:     y(i) = 23.6000; t(i) = 2.2500; i = i + 1
 85:     y(i) = 31.0500; t(i) = 1.7500; i = i + 1
 86:     y(i) = 23.7750; t(i) = 2.2500; i = i + 1
 87:     y(i) = 17.7375; t(i) = 2.7500; i = i + 1
 88:     y(i) = 13.8000; t(i) = 3.2500; i = i + 1
 89:     y(i) = 11.5875; t(i) = 3.7500; i = i + 1
 90:     y(i) = 9.4125; t(i) = 4.2500; i = i + 1
 91:     y(i) = 7.7250; t(i) = 4.7500; i = i + 1
 92:     y(i) = 7.3500; t(i) = 5.2500; i = i + 1
 93:     y(i) = 8.0250; t(i) = 5.7500; i = i + 1
 94:     y(i) = 90.6000; t(i) = 0.5000; i = i + 1
 95:     y(i) = 76.9000; t(i) = 0.6250; i = i + 1
 96:     y(i) = 71.6000; t(i) = 0.7500; i = i + 1
 97:     y(i) = 63.6000; t(i) = 0.8750; i = i + 1
 98:     y(i) = 54.0000; t(i) = 1.0000; i = i + 1
 99:     y(i) = 39.2000; t(i) = 1.2500; i = i + 1
100:     y(i) = 29.3000; t(i) = 1.7500; i = i + 1
101:     y(i) = 21.4000; t(i) = 2.2500; i = i + 1
102:     y(i) = 29.1750; t(i) = 1.7500; i = i + 1
103:     y(i) = 22.1250; t(i) = 2.2500; i = i + 1
104:     y(i) = 17.5125; t(i) = 2.7500; i = i + 1
105:     y(i) = 14.2500; t(i) = 3.2500; i = i + 1
106:     y(i) = 9.4500; t(i) = 3.7500; i = i + 1
107:     y(i) = 9.1500; t(i) = 4.2500; i = i + 1
108:     y(i) = 7.9125; t(i) = 4.7500; i = i + 1
109:     y(i) = 8.4750; t(i) = 5.2500; i = i + 1
110:     y(i) = 6.1125; t(i) = 5.7500; i = i + 1
111:     y(i) = 80.0000; t(i) = 0.5000; i = i + 1
112:     y(i) = 79.0000; t(i) = 0.6250; i = i + 1
113:     y(i) = 63.8000; t(i) = 0.7500; i = i + 1
114:     y(i) = 57.2000; t(i) = 0.8750; i = i + 1
115:     y(i) = 53.2000; t(i) = 1.0000; i = i + 1
116:     y(i) = 42.5000; t(i) = 1.2500; i = i + 1
117:     y(i) = 26.8000; t(i) = 1.7500; i = i + 1
118:     y(i) = 20.4000; t(i) = 2.2500; i = i + 1
119:     y(i) = 26.8500; t(i) = 1.7500; i = i + 1
120:     y(i) = 21.0000; t(i) = 2.2500; i = i + 1
121:     y(i) = 16.4625; t(i) = 2.7500; i = i + 1
122:     y(i) = 12.5250; t(i) = 3.2500; i = i + 1
123:     y(i) = 10.5375; t(i) = 3.7500; i = i + 1
124:     y(i) = 8.5875; t(i) = 4.2500; i = i + 1
125:     y(i) = 7.1250; t(i) = 4.7500; i = i + 1
126:     y(i) = 6.1125; t(i) = 5.2500; i = i + 1
127:     y(i) = 5.9625; t(i) = 5.7500; i = i + 1
128:     y(i) = 74.1000; t(i) = 0.5000; i = i + 1
129:     y(i) = 67.3000; t(i) = 0.6250; i = i + 1
130:     y(i) = 60.8000; t(i) = 0.7500; i = i + 1
131:     y(i) = 55.5000; t(i) = 0.8750; i = i + 1
132:     y(i) = 50.3000; t(i) = 1.0000; i = i + 1
133:     y(i) = 41.0000; t(i) = 1.2500; i = i + 1
134:     y(i) = 29.4000; t(i) = 1.7500; i = i + 1
135:     y(i) = 20.4000; t(i) = 2.2500; i = i + 1
136:     y(i) = 29.3625; t(i) = 1.7500; i = i + 1
137:     y(i) = 21.1500; t(i) = 2.2500; i = i + 1
138:     y(i) = 16.7625; t(i) = 2.7500; i = i + 1
139:     y(i) = 13.2000; t(i) = 3.2500; i = i + 1
140:     y(i) = 10.8750; t(i) = 3.7500; i = i + 1
141:     y(i) = 8.1750; t(i) = 4.2500; i = i + 1
142:     y(i) = 7.3500; t(i) = 4.7500; i = i + 1
143:     y(i) = 5.9625; t(i) = 5.2500; i = i + 1
144:     y(i) = 5.6250; t(i) = 5.7500; i = i + 1
145:     y(i) = 81.5000; t(i) = .5000; i = i + 1
146:     y(i) = 62.4000; t(i) = .7500; i = i + 1
147:     y(i) = 32.5000; t(i) = 1.5000; i = i + 1
148:     y(i) = 12.4100; t(i) = 3.0000; i = i + 1
149:     y(i) = 13.1200; t(i) = 3.0000; i = i + 1
150:     y(i) = 15.5600; t(i) = 3.0000; i = i + 1
151:     y(i) = 5.6300; t(i) = 6.0000; i = i + 1
152:     y(i) = 78.0000; t(i) = .5000; i = i + 1
153:     y(i) = 59.9000; t(i) = .7500; i = i + 1
154:     y(i) = 33.2000; t(i) = 1.5000; i = i + 1
155:     y(i) = 13.8400; t(i) = 3.0000; i = i + 1
156:     y(i) = 12.7500; t(i) = 3.0000; i = i + 1
157:     y(i) = 14.6200; t(i) = 3.0000; i = i + 1
158:     y(i) = 3.9400; t(i) = 6.0000; i = i + 1
159:     y(i) = 76.8000; t(i) = .5000; i = i + 1
160:     y(i) = 61.0000; t(i) = .7500; i = i + 1
161:     y(i) = 32.9000; t(i) = 1.5000; i = i + 1
162:     y(i) = 13.8700; t(i) = 3.0000; i = i + 1
163:     y(i) = 11.8100; t(i) = 3.0000; i = i + 1
164:     y(i) = 13.3100; t(i) = 3.0000; i = i + 1
165:     y(i) = 5.4400; t(i) = 6.0000; i = i + 1
166:     y(i) = 78.0000; t(i) = .5000; i = i + 1
167:     y(i) = 63.5000; t(i) = .7500; i = i + 1
168:     y(i) = 33.8000; t(i) = 1.5000; i = i + 1
169:     y(i) = 12.5600; t(i) = 3.0000; i = i + 1
170:     y(i) = 5.6300; t(i) = 6.0000; i = i + 1
171:     y(i) = 12.7500; t(i) = 3.0000; i = i + 1
172:     y(i) = 13.1200; t(i) = 3.0000; i = i + 1
173:     y(i) = 5.4400; t(i) = 6.0000; i = i + 1
174:     y(i) = 76.8000; t(i) = .5000; i = i + 1
175:     y(i) = 60.0000; t(i) = .7500; i = i + 1
176:     y(i) = 47.8000; t(i) = 1.0000; i = i + 1
177:     y(i) = 32.0000; t(i) = 1.5000; i = i + 1
178:     y(i) = 22.2000; t(i) = 2.0000; i = i + 1
179:     y(i) = 22.5700; t(i) = 2.0000; i = i + 1
180:     y(i) = 18.8200; t(i) = 2.5000; i = i + 1
181:     y(i) = 13.9500; t(i) = 3.0000; i = i + 1
182:     y(i) = 11.2500; t(i) = 4.0000; i = i + 1
183:     y(i) = 9.0000; t(i) = 5.0000; i = i + 1
184:     y(i) = 6.6700; t(i) = 6.0000; i = i + 1
185:     y(i) = 75.8000; t(i) = .5000; i = i + 1
186:     y(i) = 62.0000; t(i) = .7500; i = i + 1
187:     y(i) = 48.8000; t(i) = 1.0000; i = i + 1
188:     y(i) = 35.2000; t(i) = 1.5000; i = i + 1
189:     y(i) = 20.0000; t(i) = 2.0000; i = i + 1
190:     y(i) = 20.3200; t(i) = 2.0000; i = i + 1
191:     y(i) = 19.3100; t(i) = 2.5000; i = i + 1
192:     y(i) = 12.7500; t(i) = 3.0000; i = i + 1
193:     y(i) = 10.4200; t(i) = 4.0000; i = i + 1
194:     y(i) = 7.3100; t(i) = 5.0000; i = i + 1
195:     y(i) = 7.4200; t(i) = 6.0000; i = i + 1
196:     y(i) = 70.5000; t(i) = .5000; i = i + 1
197:     y(i) = 59.5000; t(i) = .7500; i = i + 1
198:     y(i) = 48.5000; t(i) = 1.0000; i = i + 1
199:     y(i) = 35.8000; t(i) = 1.5000; i = i + 1
200:     y(i) = 21.0000; t(i) = 2.0000; i = i + 1
201:     y(i) = 21.6700; t(i) = 2.0000; i = i + 1
202:     y(i) = 21.0000; t(i) = 2.5000; i = i + 1
203:     y(i) = 15.6400; t(i) = 3.0000; i = i + 1
204:     y(i) = 8.1700; t(i) = 4.0000; i = i + 1
205:     y(i) = 8.5500; t(i) = 5.0000; i = i + 1
206:     y(i) = 10.1200; t(i) = 6.0000; i = i + 1
207:     y(i) = 78.0000; t(i) = .5000; i = i + 1
208:     y(i) = 66.0000; t(i) = .6250; i = i + 1
209:     y(i) = 62.0000; t(i) = .7500; i = i + 1
210:     y(i) = 58.0000; t(i) = .8750; i = i + 1
211:     y(i) = 47.7000; t(i) = 1.0000; i = i + 1
212:     y(i) = 37.8000; t(i) = 1.2500; i = i + 1
213:     y(i) = 20.2000; t(i) = 2.2500; i = i + 1
214:     y(i) = 21.0700; t(i) = 2.2500; i = i + 1
215:     y(i) = 13.8700; t(i) = 2.7500; i = i + 1
216:     y(i) = 9.6700; t(i) = 3.2500; i = i + 1
217:     y(i) = 7.7600; t(i) = 3.7500; i = i + 1
218:     y(i) = 5.4400; t(i) = 4.2500; i = i + 1
219:     y(i) = 4.8700; t(i) = 4.7500; i = i + 1
220:     y(i) = 4.0100; t(i) = 5.2500; i = i + 1
221:     y(i) = 3.7500; t(i) = 5.7500; i = i + 1
222:     y(i) = 24.1900; t(i) = 3.0000; i = i + 1
223:     y(i) = 25.7600; t(i) = 3.0000; i = i + 1
224:     y(i) = 18.0700; t(i) = 3.0000; i = i + 1
225:     y(i) = 11.8100; t(i) = 3.0000; i = i + 1
226:     y(i) = 12.0700; t(i) = 3.0000; i = i + 1
227:     y(i) = 16.1200; t(i) = 3.0000; i = i + 1
228:     y(i) = 70.8000; t(i) = .5000; i = i + 1
229:     y(i) = 54.7000; t(i) = .7500; i = i + 1
230:     y(i) = 48.0000; t(i) = 1.0000; i = i + 1
231:     y(i) = 39.8000; t(i) = 1.5000; i = i + 1
232:     y(i) = 29.8000; t(i) = 2.0000; i = i + 1
233:     y(i) = 23.7000; t(i) = 2.5000; i = i + 1
234:     y(i) = 29.6200; t(i) = 2.0000; i = i + 1
235:     y(i) = 23.8100; t(i) = 2.5000; i = i + 1
236:     y(i) = 17.7000; t(i) = 3.0000; i = i + 1
237:     y(i) = 11.5500; t(i) = 4.0000; i = i + 1
238:     y(i) = 12.0700; t(i) = 5.0000; i = i + 1
239:     y(i) = 8.7400; t(i) = 6.0000; i = i + 1
240:     y(i) = 80.7000; t(i) = .5000; i = i + 1
241:     y(i) = 61.3000; t(i) = .7500; i = i + 1
242:     y(i) = 47.5000; t(i) = 1.0000; i = i + 1
243:     y(i) = 29.0000; t(i) = 1.5000; i = i + 1
244:     y(i) = 24.0000; t(i) = 2.0000; i = i + 1
245:     y(i) = 17.7000; t(i) = 2.5000; i = i + 1
246:     y(i) = 24.5600; t(i) = 2.0000; i = i + 1
247:     y(i) = 18.6700; t(i) = 2.5000; i = i + 1
248:     y(i) = 16.2400; t(i) = 3.0000; i = i + 1
249:     y(i) = 8.7400; t(i) = 4.0000; i = i + 1
250:     y(i) = 7.8700; t(i) = 5.0000; i = i + 1
251:     y(i) = 8.5100; t(i) = 6.0000; i = i + 1
252:     y(i) = 66.7000; t(i) = .5000; i = i + 1
253:     y(i) = 59.2000; t(i) = .7500; i = i + 1
254:     y(i) = 40.8000; t(i) = 1.0000; i = i + 1
255:     y(i) = 30.7000; t(i) = 1.5000; i = i + 1
256:     y(i) = 25.7000; t(i) = 2.0000; i = i + 1
257:     y(i) = 16.3000; t(i) = 2.5000; i = i + 1
258:     y(i) = 25.9900; t(i) = 2.0000; i = i + 1
259:     y(i) = 16.9500; t(i) = 2.5000; i = i + 1
260:     y(i) = 13.3500; t(i) = 3.0000; i = i + 1
261:     y(i) = 8.6200; t(i) = 4.0000; i = i + 1
262:     y(i) = 7.2000; t(i) = 5.0000; i = i + 1
263:     y(i) = 6.6400; t(i) = 6.0000; i = i + 1
264:     y(i) = 13.6900; t(i) = 3.0000; i = i + 1
265:     y(i) = 81.0000; t(i) = .5000; i = i + 1
266:     y(i) = 64.5000; t(i) = .7500; i = i + 1
267:     y(i) = 35.5000; t(i) = 1.5000; i = i + 1
268:     y(i) = 13.3100; t(i) = 3.0000; i = i + 1
269:     y(i) = 4.8700; t(i) = 6.0000; i = i + 1
270:     y(i) = 12.9400; t(i) = 3.0000; i = i + 1
271:     y(i) = 5.0600; t(i) = 6.0000; i = i + 1
272:     y(i) = 15.1900; t(i) = 3.0000; i = i + 1
273:     y(i) = 14.6200; t(i) = 3.0000; i = i + 1
274:     y(i) = 15.6400; t(i) = 3.0000; i = i + 1
275:     y(i) = 25.5000; t(i) = 1.7500; i = i + 1
276:     y(i) = 25.9500; t(i) = 1.7500; i = i + 1
277:     y(i) = 81.7000; t(i) = .5000; i = i + 1
278:     y(i) = 61.6000; t(i) = .7500; i = i + 1
279:     y(i) = 29.8000; t(i) = 1.7500; i = i + 1
280:     y(i) = 29.8100; t(i) = 1.7500; i = i + 1
281:     y(i) = 17.1700; t(i) = 2.7500; i = i + 1
282:     y(i) = 10.3900; t(i) = 3.7500; i = i + 1
283:     y(i) = 28.4000; t(i) = 1.7500; i = i + 1
284:     y(i) = 28.6900; t(i) = 1.7500; i = i + 1
285:     y(i) = 81.3000; t(i) = .5000; i = i + 1
286:     y(i) = 60.9000; t(i) = .7500; i = i + 1
287:     y(i) = 16.6500; t(i) = 2.7500; i = i + 1
288:     y(i) = 10.0500; t(i) = 3.7500; i = i + 1
289:     y(i) = 28.9000; t(i) = 1.7500; i = i + 1
290:     y(i) = 28.9500; t(i) = 1.7500; i = i + 1

292:   end
293: end module chwirut1fmodule

295: program main
296:   use chwirut1fmodule
297:   implicit none
298: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
299: !                   Variable declarations
300: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
301: !
302: !  See additional variable declarations in the file chwirut1f.h

304:   PetscErrorCode ierr    ! used to check for functions returning nonzeros
305:   Vec x       ! solution vector
306:   Vec f       ! vector of functions
307:   Tao ta     ! Tao context
308:   PetscInt nhist
309:   PetscMPIInt size, rank    ! number of processes running
310:   PetscReal hist(100) ! objective value history
311:   PetscReal resid(100)! residual history
312:   PetscReal cnorm(100)! cnorm history
313:   PetscInt lits(100)   ! #ksp history
314:   PetscInt oh
315:   TaoConvergedReason reason

317: !  Initialize TAO and PETSc
318:   PetscCallA(PetscInitialize(ierr))

320:   PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
321:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
322:   PetscCheckA(size == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, 'This is a uniprocessor example only')

324: !  Initialize problem parameters
325:   m = 214
326:   n = 3

328: !  Allocate vectors for the solution and gradient
329:   PetscCallA(VecCreateSeq(PETSC_COMM_SELF, n, x, ierr))
330:   PetscCallA(VecCreateSeq(PETSC_COMM_SELF, m, f, ierr))

332: !  The TAO code begins here

334: !  Create TAO solver
335:   PetscCallA(TaoCreate(PETSC_COMM_SELF, ta, ierr))
336:   PetscCallA(TaoSetType(ta, TAOPOUNDERS, ierr))
337: !  Set routines for function, gradient, and hessian evaluation

339:   PetscCallA(TaoSetResidualRoutine(ta, f, FormFunction, 0, ierr))

341: !  Optional: Set initial guess
342:   call InitializeData()
343:   call FormStartingPoint(x)
344:   PetscCallA(TaoSetSolution(ta, x, ierr))

346: !  Check for TAO command line options
347:   PetscCallA(TaoSetFromOptions(ta, ierr))
348:   oh = 100
349:   PetscCallA(TaoSetConvergenceHistory(ta, hist, resid, cnorm, lits, oh, PETSC_TRUE, ierr))
350: !  SOLVE THE APPLICATION
351:   PetscCallA(TaoSolve(ta, ierr))
352:   PetscCallA(TaoGetConvergenceHistory(ta, nhist, ierr))
353:   PetscCallA(TaoGetConvergedReason(ta, reason, ierr))
354:   if (reason%v <= 0) then
355:     print *, 'Tao failed.'
356:     print *, 'Try a different TAO method, adjust some parameters,'
357:     print *, 'or check the function evaluation routines.'
358:   end if

360: !  Free TAO data structures
361:   PetscCallA(TaoDestroy(ta, ierr))

363: !  Free PETSc data structures
364:   PetscCallA(VecDestroy(x, ierr))
365:   PetscCallA(VecDestroy(f, ierr))

367:   PetscCallA(PetscFinalize(ierr))

369: end

371: !/*TEST
372: !
373: !   build:
374: !      requires: !complex
375: !
376: !   test:
377: !      args: -tao_monitor_short -tao_max_it 100 -tao_type pounders -tao_gatol 1.e-5
378: !      requires: !single
379: !
380: !TEST*/