V3FIT
fdjac_mod.f
1  MODULE fdjac_mod
2  USE stel_kinds
3 
4  INTEGER, PARAMETER :: flag_singletask = -1, flag_cleanup = -100
5  INTEGER :: mp, np, ncntp, max_processors, num_lm_params
6  REAL(rprec), DIMENSION(:), POINTER :: xp, wap
7  REAL(rprec) :: eps
8  LOGICAL, DIMENSION(:), ALLOCATABLE :: flip
9 
10  CONTAINS
11 
12  SUBROUTINE fdjac2_mp(fcn, m, n, x, fvec, fjac, ldfjac,
13  1 iflag, ncnt, epsfcn, fnorm_min, x_min, wa)
14  USE stel_kinds
15  USE mpi_params
16  USE mpi_inc
17  IMPLICIT NONE
18 C-----------------------------------------------
19 C D u m m y A r g u m e n t s
20 C-----------------------------------------------
21  INTEGER, INTENT(in) :: m, n, ldfjac, ncnt
22  INTEGER :: iflag
23  REAL(rprec), INTENT(in) :: epsfcn
24  REAL(rprec) :: x(n)
25  REAL(rprec), INTENT(in) :: fvec(m)
26  REAL(rprec) :: fjac(ldfjac,n)
27  REAL(rprec) :: fnorm_min, x_min(n), wa(m)
28  EXTERNAL fcn
29 #if defined(MPI_OPT)
30 C-----------------------------------------------
31 C L o c a l V a r i a b l e s
32 C-----------------------------------------------
33  INTEGER :: i, j, iproc_min, num_split, joff
34  INTEGER :: iflag_array(1)
35  REAL(rprec) :: epsmch, dpmpar, temp, enorm, cur_norm
36  REAL(rprec), PARAMETER :: one = 1, zero = 0
37  REAL(rprec), DIMENSION(:), ALLOCATABLE :: buffer, h,
38  1 fnorm_array !mpi stuff
39  REAL(rprec), ALLOCATABLE :: fjac_overflow(:,:), fnorm_overflow(:)
40  INTEGER :: istat, ierr_flag, ikey
41  LOGICAL :: lflip
42  LOGICAL, ALLOCATABLE :: flip_overflow(:)
43 C-----------------------------------------------
44 C E x t e r n a l F u n c t i o n s
45 C-----------------------------------------------
46  EXTERNAL dpmpar, enorm
47 C-----------------------------------------------
48 c
49 c SUBROUTINE fdjac2
50 c
51 c this SUBROUTINE computes a forward-difference approximation
52 c to the m by n jacobian matrix associated with a specified
53 c problem of m functions in n variables.
54 c
55 c the SUBROUTINE statement is
56 c
57 c SUBROUTINE fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa)
58 c
59 c WHERE
60 c
61 c fcn is the name of the user-supplied SUBROUTINE which
62 c calculates the functions. fcn must be declared
63 c in an EXTERNAL statement in the user calling
64 c program, and should be written as follows.
65 c
66 c SUBROUTINE fcn(m,n,x,fvec,iflag)
67 c INTEGER m,n,iflag
68 c REAL(rprec) x(n),fvec(m)
69 c ----------
70 c calculate the functions at x and
71 c RETURN this vector in fvec.
72 c ----------
73 c RETURN
74 c END
75 c
76 c the value of iflag should not be changed by fcn unless
77 c the user wants to terminate execution of fdjac2.
78 c in this CASE set iflag to a negative INTEGER.
79 c
80 c m is a positive INTEGER input variable set to the number
81 c of functions.
82 c
83 c n is a positive INTEGER input variable set to the number
84 c of variables. n must not exceed m.
85 c
86 c x is an input array of length n.
87 c
88 c fvec is an input array of length m which must contain the
89 c functions evaluated at x.
90 c
91 c fjac is an output m by n array which contains the
92 c approximation to the jacobian matrix evaluated at x.
93 c
94 c ldfjac is a positive INTEGER input variable not less than m
95 c which specifies the leading DIMENSION of the array fjac.
96 c
97 c iflag is an INTEGER variable which can be used to terminate
98 c the execution of fdjac2. see description of fcn.
99 c
100 c epsfcn is an input variable used in determining a suitable
101 c step length for the forward-difference approximation. this
102 c approximation assumes that the relative errors in the
103 c functions are of the order of epsfcn. If epsfcn is less
104 c than the machine precision, it is assumed that the relative
105 c errors in the functions are of the order of the machine
106 c precision.
107 c
108 c wa is a work array of length m.
109 c
110 c subprograms called
111 c
112 c user-supplied ...... fcn
113 c
114 c MINpack-supplied ... dpmpar
115 c
116 c fortran-supplied ... abs,max,sqrt
117 c
118 c argonne national laboratory. minpack project. march 1980.
119 c burton s. garbow, kenneth e. hillstrom, jorge j. more
120 c
121 c **********
122  ALLOCATE (fnorm_array(n), h(n), buffer(m), stat=istat)
123  IF (istat .ne. 0) stop 'Allocation error in fdjac2_mp'
124 c
125 c epsmch is the machine precision.
126 c
127  epsmch = dpmpar(1)
128 c
129  eps = sqrt(max(epsfcn,epsmch))
130  cur_norm = enorm(m,fvec)
131 
132  h(1:n) = eps*abs(x(1:n))
133  WHERE (h .eq. zero) h=eps
134  WHERE (flip) h = -h
135  ierr_flag = 0
136 
137 
138 ! NOTE: if numprocs > n, num_split = 0 and we go to the overflow
139 ! ("left-over") loop for num_split=mod(n,numprocs) == n
140  num_split = n/numprocs
141  mpi_comm_workers = mpi_comm_world
142  joff = 1
143 
144  DO i = 1, num_split
145  j = myid + joff
146  temp = x(j)
147  x(j) = temp + h(j)
148 
149  iflag = j
150  CALL fcn(m, n, x, wa, iflag, ncnt)
151 
152  x(j) = temp
153  buffer(:m) = (wa(:m) - fvec(:m))/h(j)
154 
155 !
156 ! STORE FNORM OF PERTURBED STATE (X + H)
157 !
158  temp = enorm(m,wa)
159  IF (temp > cur_norm) lflip = .not. flip(j)
160 
161 ! WRITE ITERATION, PROCESSOR, CHISQ TO SCREEN
162  WRITE (6, '(2x,i6,8x,i3,7x,1es12.4)') ncnt+j,
163  1 myid+1, temp**2
164  CALL flush(6)
165 
166  CALL mpi_allgather(buffer, m, mpi_real8,
167  1 fjac(1,joff), m, mpi_real8, mpi_comm_world, ierr_mpi)
168  CALL mpi_allgather(temp, 1, mpi_real8,
169  1 fnorm_array(joff), 1, mpi_real8, mpi_comm_world, ierr_mpi)
170  CALL mpi_allgather(lflip, 1, mpi_logical,
171  1 flip(joff), 1, mpi_logical, mpi_comm_world, ierr_mpi)
172 
173  joff = joff + numprocs
174 
175  IF (ierr_mpi .ne. 0) GO TO 100
176  END DO
177 
178 !
179 ! account for left-over variables OR special case when numprocs>n:
180 ! must redefine the MPIC_COMM_WORDERS group to consist of MOD(n,numprocs) only
181 !
182  num_split = mod(n, numprocs) !overflow...
183  IF (num_split .ne. 0) THEN
184 
185  ALLOCATE(fjac_overflow(m,numprocs), fnorm_overflow(numprocs),
186  1 flip_overflow(numprocs), stat = istat)
187  IF (istat .ne. 0) THEN
188  print *,' istat = ',istat,' in fdjac_mod, myid: ', myid
189  stop
190  END IF
191  j = myid + joff
192  IF (j .le. n) THEN
193  temp = x(j)
194  x(j) = temp + h(j)
195  ikey = worker_split_key+1
196  ELSE
197  temp = 0
198  ikey = mpi_undefined
199  END IF
200 
201 ! Set new WORKERS communicator to include first MOD(n,numprocs) only
202  CALL mpi_comm_split(mpi_comm_world, ikey, worker_id,
203  1 mpi_comm_workers, ierr_mpi)
204  IF (ierr_mpi .ne. 0)
205  1 stop 'MPI_COMM_SPLIT error in fdjac_mp_queue'
206 
207 
208  IF (j .le. n) THEN !(MPI_COMM_WORKERS ONLY ACCESS HERE)
209  iflag = j
210  CALL fcn(m, n, x, wa, iflag, ncnt)
211  x(j) = temp
212  buffer(:m) = (wa(:m) - fvec(:m))/h(j)
213  temp = enorm(m,wa)
214  IF (temp > cur_norm) lflip = .not. flip(j)
215 
216  WRITE (6, '(2x,i6,8x,i3,7x,1es12.4)') ncnt+j,
217  1 myid+1, temp**2
218 
219  CALL mpi_comm_free(mpi_comm_workers, ierr_mpi)
220  ELSE !(ALL OTHER PROCESSORS GO HERE)
221  buffer = 0; lflip = .true.
222  END IF
223 
224  CALL mpi_allgather(buffer, m, mpi_real8,
225  1 fjac_overflow, m, mpi_real8, mpi_comm_world, ierr_mpi)
226  CALL mpi_allgather(temp, 1, mpi_real8,
227  1 fnorm_overflow, 1, mpi_real8, mpi_comm_world, ierr_mpi)
228  CALL mpi_allgather(lflip, 1, mpi_logical,
229  1 flip_overflow, 1, mpi_logical, mpi_comm_world, ierr_mpi)
230 
231  fjac(:,joff:n) = fjac_overflow(:,1:num_split)
232  fnorm_array(joff:n) = fnorm_overflow(1:num_split)
233  flip(joff:n) = flip_overflow(1:num_split)
234 
235  DEALLOCATE(fjac_overflow, fnorm_overflow, flip_overflow)
236 
237  END IF
238 
239 !
240 ! Find processor with minimum fnorm_min value and broadcast wa (=fvec_min), x_min, fnorm_min
241 !
242  fnorm_min = minval(fnorm_array)
243  iflag_array = minloc(fnorm_array)
244  iproc_min = iflag_array(1)
245  IF (iproc_min .le. 0 .or. iproc_min .gt. n) THEN
246  print *,'IPROC_MIN=',iproc_min,' out of range in fdjac2_mp'
247  stop
248  END IF
249  wa(:) = fjac(:, iproc_min)*h(iproc_min) + fvec(:)
250  x_min(:) = x(:)
251  x_min(iproc_min) = x(iproc_min) + h(iproc_min)
252 
253  DEALLOCATE (h, fnorm_array, buffer)
254 !
255 ! Do any special cleanup now for iflag = flag_cleanup
256 ! Barrier appears in fcn call
257 !
258  ikey = flag_cleanup
259  mpi_comm_workers = mpi_comm_world
260  CALL fcn(m, n, x, wa, ikey, ncnt)
261 
262  RETURN
263 
264  100 CONTINUE
265  WRITE (6, *) ' MPI_BCAST error in FDJAC2_MP: IERR=', ierr_mpi,
266  1 ' PROCESSOR: ',myid
267 #endif
268  END SUBROUTINE fdjac2_mp
269 
270 
271  SUBROUTINE fdjac2_mp_queue(fcn, m, n, x, fvec, fjac, ldfjac,
272  1 iflag, ncnt, epsfcn, fnorm_min, x_min, wa)
273  USE stel_kinds
274  USE mpi_params
275  USE mpi_inc
276  IMPLICIT NONE
277 C-----------------------------------------------
278 C D u m m y A r g u m e n t s
279 C-----------------------------------------------
280  INTEGER, INTENT(in) :: m, n, ldfjac, ncnt
281  INTEGER :: iflag
282  REAL(rprec), INTENT(in) :: epsfcn
283  REAL(rprec) :: x(n)
284  REAL(rprec), INTENT(in) :: fvec(m)
285  REAL(rprec) :: fjac(ldfjac,n)
286  REAL(rprec) :: fnorm_min, x_min(n), wa(m)
287  EXTERNAL fcn
288 #if defined(MPI_OPT)
289 C-----------------------------------------------
290 C L o c a l V a r i a b l e s
291 C-----------------------------------------------
292  INTEGER :: i, j, iproc_min
293  INTEGER :: iflag_array(1)
294  REAL(rprec) :: epsmch, dpmpar, temp, enorm, cur_norm
295  REAL(rprec), PARAMETER :: one = 1, zero = 0
296  REAL(rprec), DIMENSION(:), ALLOCATABLE :: buffer, h,
297  1 fnorm_array !mpi stuff
298  INTEGER :: status(MPI_STATUS_size) !mpi stuff
299  INTEGER :: numsent, sender, istat, ikey !mpi stuff
300  INTEGER :: anstype, column, ierr_flag !mpi stuff
301 C-----------------------------------------------
302 C E x t e r n a l F u n c t i o n s
303 C-----------------------------------------------
304  EXTERNAL dpmpar, enorm
305 C-----------------------------------------------
306 c
307 c SUBROUTINE fdjac2
308 c
309 c this SUBROUTINE computes a forward-difference approximation
310 c to the m by n jacobian matrix ASSOCIATED with a specified
311 c problem of m functions in n variables.
312 c
313 c the SUBROUTINE statement is
314 c
315 c SUBROUTINE fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa)
316 c
317 c WHERE
318 c
319 c fcn is the name of the user-supplied SUBROUTINE which
320 c calculates the functions. fcn must be declared
321 c in an EXTERNAL statement in the user calling
322 c program, and should be written as follows.
323 c
324 c SUBROUTINE fcn(m,n,x,fvec,iflag)
325 c INTEGER m,n,iflag
326 c REAL(rprec) x(n),fvec(m)
327 c ----------
328 c calculate the functions at x and
329 c RETURN this vector in fvec.
330 c ----------
331 c RETURN
332 c END
333 c
334 c the value of iflag should not be changed by fcn unless
335 c the user wants to terminate execution of fdjac2.
336 c in this case set iflag to a negative integer.
337 c
338 c m is a positive INTEGER input variable set to the number
339 c of functions.
340 c
341 c n is a positive INTEGER input variable set to the number
342 c of variables. n must not exceed m.
343 c
344 c x is an input array of length n.
345 c
346 c fvec is an input array of length m which must contain the
347 c functions evaluated at x.
348 c
349 c fjac is an output m by n array which contains the
350 c approximation to the jacobian matrix evaluated at x.
351 c
352 c ldfjac is a positive INTEGER input variable not less than m
353 c which specifies the leading DIMENSION of the array fjac.
354 c
355 c iflag is an INTEGER variable which can be used to terminate
356 c the execution of fdjac2. see description of fcn.
357 c
358 c epsfcn is an input variable used in determining a suitable
359 c step length for the forward-difference approximation. this
360 c approximation assumes that the relative errors in the
361 c functions are of the order of epsfcn. If epsfcn is less
362 c than the machine precision, it is assumed that the relative
363 c errors in the functions are of the order of the machine
364 c precision.
365 c
366 c wa is a work array of length m.
367 c
368 c subprograms called
369 c
370 c user-supplied ...... fcn
371 c
372 c MINpack-supplied ... dpmpar
373 c
374 c fortran-supplied ... abs,max,sqrt
375 c
376 c argonne national laboratory. minpack project. march 1980.
377 c burton s. garbow, kenneth e. hillstrom, jorge j. more
378 c
379 c **********
380  ALLOCATE (buffer(m), fnorm_array(n), h(n), stat=istat)
381  IF (istat .ne. 0) stop 'Allocation error in fdjac2_mp'
382  ierr_flag = 0
383 c
384 c epsmch is the machine precision.
385 c
386  epsmch = dpmpar(1)
387 c
388  eps = sqrt(max(epsfcn,epsmch))
389 
390 c Set up WORKER COMMUNICATOR
391  IF (myid .eq. master) THEN
392  ikey = mpi_undefined
393  ELSE
394  ikey = worker_split_key+1
395  END IF
396  CALL mpi_comm_split(mpi_comm_world, ikey, worker_id,
397  1 mpi_comm_workers, ierr_mpi)
398  IF (ierr_mpi .ne. 0) stop 'MPI_COMM_SPLIT error in fdjac_mp_queue'
399 
400 c
401 c Begin parallel MPI master/worker version of the function call. A
402 c bank queue or master/worker MPI algorithm is used here to achieve
403 c parallel load balancing in the somewhat uneven work load involved
404 c in calculating the Jacobian function needed in the Levenberg-Marquardt
405 c optimization. This model is based on the example given in Chapt. 5 of
406 c "The LAM Companion to Using MPI" by Zdzislaw Meglicki (online see:
407 c http://carpanta.dc.fi.udc.es/docs/mpi/mpi-lam/mpi.html). These
408 c modifications were made by D.A. Spong 8/23/2000.
409 c
410 c ****Master portion of the code****
411 c
412  IF (myid .eq. master) THEN
413  numsent = 0 !numsent is a counter used to track how many
414  !jobs have been sent to workers
415  cur_norm = enorm(m,fvec)
416 c
417 c calculate the displacements
418 c
419  h(1:n) = eps*abs(x(1:n))
420  WHERE (h .eq. zero) h=eps
421  WHERE (flip) h = -h
422 
423 c Send forward difference displacements from master to each
424 c worker process. Tag these with the column number (j)
425 c
426  DO j = 1,min(numprocs-1,n)
427  numsent = numsent+1
428  temp = x(numsent)
429  x(numsent) = temp + h(numsent)
430  CALL mpi_send(x, n, mpi_real8, j,
431  1 numsent, mpi_comm_world, ierr_mpi)
432  IF (ierr_mpi .ne. 0) stop 'MPI_SEND error(1) in fdjac2'
433  x(numsent) = temp
434  END DO !j = 1,MIN(numprocs-1,n)
435 c
436 c Looping through the columns, collect answers from the workers.
437 c As answers are received, new uncalculated columns are sent
438 c out to these same workers.
439 c
440  DO j = 1,n
441  CALL mpi_recv(wa, m, mpi_real8,
442  1 mpi_any_source, mpi_any_tag,
443  2 mpi_comm_world, status, ierr_mpi)
444  IF (ierr_mpi .ne. 0) stop 'MPI_RECV error(1) in fdjac2'
445  sender = status(mpi_source)
446  anstype = status(mpi_tag) ! column is tag value
447  IF (anstype .gt. n) stop 'ANSTYPE > N IN FDJAC2'
448 
449  fjac(:m,anstype) = (wa(:m) - fvec(:m))/h(anstype)
450 !
451 ! STORE FNORM OF PERTURBED STATE (X + H)
452 !
453  temp = enorm(m,wa)
454  fnorm_array(anstype) = temp
455  IF (temp > cur_norm) flip(anstype) = .not. flip(anstype)
456 
457 ! WRITE ITERATION, PROCESSOR, CHISQ TO SCREEN
458  WRITE (6, '(2x,i6,8x,i3,7x,1es12.4)') ncnt+anstype,
459  1 sender, temp**2
460 
461 c
462 c If more columns are left, then send another column to the worker(sender)
463 c that just sent in an answer
464 c
465  IF (numsent .lt. n) THEN
466  numsent = numsent+1
467  temp = x(numsent)
468  x(numsent) = temp + h(numsent)
469 
470  CALL mpi_send(x, n, mpi_real8,
471  1 sender, numsent, mpi_comm_world, ierr_mpi)
472  IF (ierr_mpi .ne. 0) stop 'MPI_SEND error(2) in fdjac2'
473  x(numsent) = temp
474 
475  ELSE ! Tell workers that there is no more work to do
476 
477  CALL mpi_send(mpi_bottom, 0, mpi_real8,
478  1 sender, 0, mpi_comm_world, ierr_mpi)
479  IF (ierr_mpi .ne. 0) stop 'MPI_end error(3) in fdjac2'
480  ENDIF ! IF( myid .eq. master ) THEN
481  END DO ! DO j = 1,n
482 c
483 c ****Worker portion of the code****
484 c Skip this when processor id exceeds work to be done
485 c
486  ELSE IF (myid .le. n) THEN ! IF( myid .ne. master )
487 c
488 c Otherwise accept the next available column, check the tag,
489 c and if the tag is non-zero call subroutine fcn.
490 c If the tag is zero, there are no more columns
491 c and worker skips to the end.
492 c
493  90 CALL mpi_recv(x, n, mpi_real8, master,
494  1 mpi_any_tag, mpi_comm_world, status, ierr_mpi)
495  IF (ierr_mpi .ne. 0) stop 'MPI_RECV error(2) in fdjac2'
496 
497  column = status(mpi_tag)
498  IF (column .ne. 0) THEN
499  iflag = column
500 c Call the chisq fcn for the portion of displacement vector which
501 c was just received. Note that WA stores the local fvec_min array
502  CALL fcn(m, n, x, wa, iflag, ncnt)
503 
504  IF (iflag.ne.0 .and. ierr_flag.eq.0) ierr_flag = iflag
505 c
506 c Send this function evaluation back to the master process tagged
507 c with the column number so the master knows where to put it
508 c
509  CALL mpi_send(wa, m, mpi_real8, master,
510  1 column, mpi_comm_world, ierr_mpi)
511  IF (ierr_mpi .ne. 0) stop 'MPI_SEND error(4) in fdjac2'
512  GO TO 90 !Return to 90 and check if master process has sent any more jobs
513  END IF
514  ENDIF ! IF( myid .ne. master )
515 
516 !
517 ! Broadcast the fjac matrix
518 !
519  DO j=1,n
520  IF (myid .eq. master) buffer(:m) = fjac(:m,j)
521  CALL mpi_bcast(buffer, m, mpi_real8, master,
522  1 mpi_comm_world, ierr_mpi)
523  IF (ierr_mpi .ne. 0) GO TO 100
524  IF (myid .ne. master) fjac(:m,j) = buffer(:m)
525  END DO
526 
527 !
528 ! Find processor with minimum fnorm_min value and broadcast wa (=fvec_min), x_min, fnorm_min
529 !
530  IF (myid .eq. master) THEN
531  fnorm_min = minval(fnorm_array)
532  iflag_array = minloc(fnorm_array)
533  iproc_min = iflag_array(1)
534  IF (iproc_min .le. 0 .or. iproc_min .gt. n) THEN
535  print *,'IPROC_MIN=',iproc_min,' out of range in fdjac2_mp'
536  stop
537  END IF
538  wa(:) = fjac(:, iproc_min)*h(iproc_min) + fvec(:)
539  x_min(:) = x(:)
540  x_min(iproc_min) = x(iproc_min) + h(iproc_min)
541  END IF
542 
543  CALL mpi_bcast(fnorm_min, 1, mpi_real8, master,
544  1 mpi_comm_world, ierr_mpi)
545  IF (ierr_mpi .ne. 0) GOTO 100
546  CALL mpi_bcast(wa, m, mpi_real8, master, mpi_comm_world,ierr_mpi)
547  IF (ierr_mpi .ne. 0) GOTO 100
548  CALL mpi_bcast(x_min, n, mpi_real8, master, mpi_comm_world,
549  1 ierr_mpi)
550  IF (ierr_mpi .ne. 0) GOTO 100
551 c
552 c Original serial version of the Jacobian calculation:
553 c
554 c DO j = 1, n
555 c temp = x(j)
556 c h = eps*ABS(temp)
557 c IF (h .eq. zero) h = eps
558 c x(j) = temp + h
559 c CALL fcn (m, n, x, wa, iflag)
560 c IF (iflag .lt. 0) EXIT
561 c x(j) = temp
562 c fjac(:m,j) = (wa - fvec)/h
563 c END DO
564 
565  DEALLOCATE (h, buffer, fnorm_array)
566 !
567 ! Do any special cleanup now for iflag = flag_cleanup
568 ! Barrier appears in fcn call
569 !
570  column = flag_cleanup
571  CALL fcn(m, n, x, wa, column, ncnt)
572 
573 !
574 ! Reassign initial x value to all processors and perform error handling
575 ! This is necessary because other processors had x + h in them
576 !
577  CALL mpi_bcast(x, n, mpi_real8, master, mpi_comm_world, ierr_mpi)
578  IF (ierr_mpi .ne. 0) GOTO 100
579 
580  IF (ierr_flag .ne. 0) THEN
581  iflag = ierr_flag
582  CALL mpi_bcast(iflag, 1, mpi_integer, myid,
583  1 mpi_comm_world, ierr_mpi)
584  IF (ierr_mpi .ne. 0) GOTO 100
585  END IF
586 
587  IF (myid .ne. master)
588  1 CALL mpi_comm_free(mpi_comm_workers, ierr_mpi)
589 
590  RETURN
591 
592  100 CONTINUE
593  WRITE (6, *) ' MPI_BCAST error in FDJAC2_MP: IERR=', ierr_mpi,
594  1 ' PROCESSOR: ',myid
595 #endif
596  END SUBROUTINE fdjac2_mp_queue
597 
598 
599  SUBROUTINE fdjac2(fcn, m, n, x, fvec, fjac, ldfjac, iflag, ncnt,
600  1 epsfcn, wa, time, fnorm_min, x_min, fvec_min)
601  USE stel_kinds
602 #if !defined(MPI_OPT)
603  IMPLICIT NONE
604 C-----------------------------------------------
605 C D u m m y A r g u m e n t s
606 C-----------------------------------------------
607  INTEGER, INTENT(in) :: m, n, ldfjac, ncnt
608  INTEGER, TARGET :: iflag
609  REAL(rprec) :: epsfcn, time
610  REAL(rprec), TARGET :: x(n), wa(m)
611  REAL(rprec), INTENT(in) :: fvec(m)
612  REAL(rprec), INTENT(out) :: fjac(ldfjac,n)
613  REAL(rprec), INTENT(out) :: fnorm_min, x_min(n), fvec_min(m)
614 C-----------------------------------------------
615 C L o c a l P a r a m e t e r s
616 C-----------------------------------------------
617  REAL(rprec), PARAMETER :: zero = 0
618 C-----------------------------------------------
619 C L o c a l V a r i a b l e s
620 C-----------------------------------------------
621  INTEGER :: j, istat, iread, ic1, ic2, irate, count_max
622  REAL(rprec) :: epsmch, h, dpmpar, temp, cur_norm
623 C-----------------------------------------------
624 C E x t e r n a l F u n c t i o n s
625 C-----------------------------------------------
626  EXTERNAL fcn, dpmpar, multiprocess, fdjac_parallel
627  REAL(rprec), EXTERNAL :: enorm
628 C-----------------------------------------------
629 c
630 c SUBROUTINE fdjac2
631 c
632 c this SUBROUTINE computes a forward-difference approximation
633 c to the m by n jacobian matrix ASSOCIATED with a specified
634 c problem of m functions in n variables.
635 c
636 c Here
637 c
638 c fcn is the name of the user-supplied SUBROUTINE which
639 c calculates the functions. fcn must be declared
640 c in an external statement in the user calling
641 c program (see LMDIF1 for documentation), and should be written as follows:
642 c
643 c SUBROUTINE fcn(m,n,x,fvec,iflag,ncnt)
644 c INTEGER m,n,iflag
645 c REAL(rprec) x(n),fvec(m)
646 c ----------
647 c calculate the functions at x and
648 c RETURN this vector in fvec.
649 c ----------
650 c RETURN
651 c END
652 c
653 c fjac is an output m by n array which CONTAINS the
654 c approximation to the jacobian matrix evaluated at x.
655 c
656 c ldfjac is a positive INTEGER input variable not less than m
657 c which specifies the leading DIMENSION of the array fjac.
658 c
659 c iflag is an INTEGER variable which can be used to terminate
660 c the execution of fdjac2. see description of fcn.
661 c
662 c epsfcn is an input variable used in determining a suitable
663 c step length for the forward-difference approximation. this
664 c approximation assumes that the relative errors in the
665 c functions are of the order of epsfcn. IF epsfcn is less
666 c than the machine precision, it is assumed that the relative
667 c errors in the functions are of the order of the machine
668 c precision.
669 c
670 c wa is a work array of length m.
671 c
672 c subprograms called
673 c
674 c user-supplied ...... fcn
675 c
676 c MINpack-supplied ... dpmpar
677 c
678 c fortran-supplied ... ABS,max,sqrt
679 c
680 c argonne national laboratory. MINpack project. march 1980.
681 c burton s. garbow, kenneth e. hillstrom, jorge j. more
682 c
683 c **********
684 
685 c
686 c epsmch is the machine precision.
687 c
688  epsmch = dpmpar(1)
689 c
690  eps = sqrt(max(epsfcn,epsmch))
691 !
692 ! Load fdjac_mod values. pointers will automatically update targets...
693 ! Prepare for multi-processing...
694 !
695  mp = m
696  np = n
697  ncntp = ncnt
698  xp => x
699  wap => wa
700 
701 ! Find min chisq = fnorm**2 state for this jacobian evaluation
702 ! (Do NOT retain from previous evaluations, or could get into a non-converging loop...)
703 
704  fnorm_min = huge(fnorm_min)
705 
706  CALL system_clock(ic1, irate)
707  CALL multiprocess(n, max_processors, fdjac_parallel, fcn)
708  CALL system_clock(ic2, irate, count_max)
709  IF (ic2 .lt. ic1) ic2 = ic2 + count_max
710 
711 c DO j = 1, n
712 c temp = x(j)
713 c h = eps*ABS(temp)
714 c IF (h .eq. zero) h = eps
715 c x(j) = temp + h
716 c CALL fcn (m, n, x, wa, iflag, ncnt)
717 c IF (iflag .lt. 0) EXIT
718 c x(j) = temp
719 c fjac(:m,j) = (wa - fvec)/h
720 c END DO
721 
722 
723  cur_norm = enorm(m,fvec) ! where are we now?
724 
725  DO j = 1, n
726 
727  READ (j+1000, iostat=iread) istat, iflag, h, temp
728  IF (iread .ne. 0) THEN
729  WRITE (6, *) 'Error reading from file fort.', j+1000,
730  1 ' in fdjac2: IOSTAT = ', iread
731  iflag = -14
732  ELSE IF (j .ne. istat) THEN
733  WRITE (6, *) 'Wrong value for INDEX j READ in fdjac2'
734  iflag = -14
735  END IF
736 
737  IF (iflag .ne. 0) EXIT
738 
739 #if defined(CRAY)
740  DO k = 1, m
741  READ (j+1000) wa(k)
742  END DO
743 #else
744  READ (j+1000) wa
745 #endif
746  fjac(:m,j) = (wa - fvec)/h
747 
748  IF( temp > cur_norm) flip(j) = .not. flip(j) ! flip for next time
749 
750  IF( temp < fnorm_min) THEN
751  fnorm_min = temp
752  fvec_min = wa
753 #if defined(CRAY)
754  DO k = 1, n
755  READ (j+1000) x_min(k)
756  END DO
757 #else
758  READ (j+1000) x_min
759 #endif
760  END IF
761 
762  CLOSE (j+1000, status='delete') !!Needed to run correctly in multi-tasking...
763 
764  END DO
765 
766 !
767 ! Do any special cleanup now for iflag = flag_cleanup
768 !
769  iflag = flag_cleanup
770  CALL fcn(m, n, x, wa, iflag, ncnt)
771 
772  time = time + real(ic2 - ic1)/real(irate) !!Time in multi-process CALL
773 #endif
774  END SUBROUTINE fdjac2
775 
776 
777  END MODULE fdjac_mod
778 
mpi_inc
Umbrella module avoid multiple inlcudes of the mpif.h header.
Definition: mpi_inc.f:11