V3FIT
•All Classes Namespaces Files Functions Variables Enumerations Macros Pages
lmpar_mod.f
1  MODULE lmpar_mod
2  USE stel_kinds
3  INTEGER :: nscan, ldfjac
4  INTEGER, DIMENSION(:), POINTER :: ipvt
5  REAL(rprec) :: pnorm, fnorm1, delta, par, spread_ratio
6  REAL(rprec), DIMENSION(:), POINTER :: wa2p, wa3p, wa4p
7  REAL(rprec), DIMENSION(:), POINTER :: diag, qtf
8  REAL(rprec), DIMENSION(:,:), POINTER :: fjac
9  LOGICAL :: lfirst_lm
10 
11  CONTAINS
12 
13  SUBROUTINE levmarq_param_mp(x, wa1, wa2, wa3, wa4,
14  1 nfev, m, n, iflag, fcn)
15  USE fdjac_mod, ONLY: flag_cleanup
16  USE mpi_params
17  USE mpi_inc
18  IMPLICIT NONE
19 C-----------------------------------------------
20 C D u m m y A r g u m e n t s
21 C-----------------------------------------------
22  INTEGER, INTENT(in) :: n, m
23  INTEGER :: iflag, nfev
24  REAL(rprec), INTENT(in) :: x(n)
25  REAL(rprec) :: wa1(n), wa2(n), wa3(n), wa4(m)
26  EXTERNAL fcn
27 #if defined(MPI_OPT)
28 C-----------------------------------------------
29 C L o c a l P a r a m e t e r s
30 C-----------------------------------------------
31  REAL(rprec), PARAMETER :: zero = 0, one = 1
32  REAL(rprec), DIMENSION(11), PARAMETER :: factors =
33  1 (/ 1.0_dp, 0.5_dp, 0.25_dp, 0.128_dp, 0.75_dp,
34  2 1.25_dp, 1.5_dp, 0.9_dp, 1.1_dp, 1.75_dp, 2.1_dp /)
35 C-----------------------------------------------
36 C L o c a l V a r i a b l e s
37 C-----------------------------------------------
38  INTEGER :: iproc, iproc_min, nfact, num_lev, istat, ierr, j,
39  1 iflag_min(1)
40  INTEGER, ALLOCATABLE, DIMENSION(:) :: iflag_array
41  REAL(rprec) :: scale_factor
42  REAL(rprec), DIMENSION(:), ALLOCATABLE :: fnorm_array
43  CHARACTER(LEN=1) :: ext, low_mark
44 C-----------------------------------------------
45 C E x t e r n a l F u n c t i o n s
46 C-----------------------------------------------
47  REAL(rprec), EXTERNAL :: enorm
48 C-----------------------------------------------
49 ! Perform Numprocs different function calls (in parallel) to find the minimum norm (chi-sq).
50 ! MPI calls are used to determine which processor has the minimum norm and then the
51 ! information is sent to all processors using MPI Broadcasts (modifications made by D. A. Spong 8/27/2000).
52 !
53 ! Uses processors 0,...,numprocs-1 (which belong to the user-defined MPI_COMM_WORKERS communicator)
54 
55  nfact = SIZE(factors)
56  num_lev = numprocs
57  ALLOCATE (iflag_array(numprocs), fnorm_array(numprocs),
58  1 stat=istat)
59  IF (istat .ne. 0) stop 'Allocation error in levmarq_param_mp'
60 
61  iproc = myid
62  IF (lfirst_lm .and. num_lev > 2) THEN
63 !
64 ! Do an exponential spread the first call (lfirst_lm=.true.) to see where we are
65 !
66  scale_factor = 10._dp**(-iproc)
67 ! scale_factor = 10._dp**(one-iproc)
68  lfirst_lm = .false.
69  ELSE IF (num_lev > 2*nfact) THEN
70  scale_factor = ((iproc+1)*maxval(factors))/num_lev
71  ELSE IF (iproc .lt. nfact) THEN
72  scale_factor = factors(iproc+1)
73  ELSE
74  scale_factor = ((iproc-nfact)*minval(factors))/
75  1 (num_lev-nfact)
76  END IF
77 
78  delta = scale_factor*delta
79 
80  CALL lmpar (n, fjac, ldfjac, ipvt, diag, qtf,
81  1 delta, par, wa1, wa2, wa3, wa4)
82 !
83 ! store the direction p and x + p. calculate the norm of p.
84 !
85  IF (par .eq. zero) wa1 = wa1*scale_factor
86  wa1 = -wa1
87  wa2 = x + wa1
88  wa3 = diag*wa1
89  pnorm = enorm(n,wa3)
90 !
91 ! evaluate the function at x + p and calculate its norm.
92 ! Only do for 0 < myid < n processors (the MPI_COMM_WORKERS group),
93 ! which the v3post routines are expecting. To use other processors, must
94 ! pass the communicator id (maybe through hi bits of iflag, using IAND...)
95 !
96  mpi_comm_workers = mpi_comm_world
97  iflag = iproc
98  CALL fcn (m, n, wa2, wa4, iflag, nfev)
99 
100 !
101 ! Gather iflag information to all processors and check for iflag < 0
102 !
103  CALL mpi_allgather(iflag, 1, mpi_integer, iflag_array, 1,
104  1 mpi_integer, mpi_comm_world, ierr)
105  IF (ierr .ne. 0) stop 'MPI_ALLGATHER failed in levmarq_param_mp'
106 
107  fnorm1 = enorm(m,wa4)
108 
109 !
110 ! Find processor with minimum fnorm1 value
111 !
112  CALL mpi_allgather(fnorm1, 1, mpi_real8, fnorm_array, 1,
113  1 mpi_real8, mpi_comm_world, ierr)
114  IF (ierr .ne. 0) stop 'MPI_ALLGATHER failed in LMDIF'
115  iflag_min = minloc(fnorm_array)
116  iproc_min = iflag_min(1) - 1
117 
118 ! iflag = MINVAL(iflag_array)
119 ! IF (iflag .lt. 0) GOTO 100
120 
121  ext = ' '
122  low_mark = ' '
123  IF (myid .eq. master) ext = '*'
124  IF (iproc .eq. iproc_min) low_mark = '*'
125  iproc = iproc+1
126  WRITE(6, '(2x,i6,8x,i3,4x,2(3x,1es12.4,a),3x,1es12.4)')
127  1 iproc+nfev, iproc, fnorm1**2, low_mark, par, ext,
128  2 delta
129 
130  CALL flush(6)
131 
132  fnorm1 = minval(fnorm_array)
133 
134  100 CONTINUE
135 
136  DEALLOCATE (iflag_array, fnorm_array)
137 
138 !
139 ! Broadcast all relevant scalars and arrays from the
140 ! processor with minimum fnorm1 to the other processors,
141 ! overwriting their data. Note: diag, ipvt are same already on
142 ! ALL processors. wa3 is overwritten...
143 !
144  CALL mpi_bcast(pnorm,1,mpi_real8,iproc_min,
145  1 mpi_comm_world,ierr)
146  IF (ierr .ne. 0) GOTO 3000
147  CALL mpi_bcast(par,1,mpi_real8,iproc_min,
148  1 mpi_comm_world,ierr)
149  IF (ierr .ne. 0) GOTO 3000
150  CALL mpi_bcast(delta,1,mpi_real8,iproc_min,
151  1 mpi_comm_world,ierr)
152  IF (ierr .ne. 0) GOTO 3000
153  CALL mpi_bcast(wa1,n,mpi_real8,iproc_min,
154  1 mpi_comm_world,ierr)
155  IF (ierr .ne. 0) GOTO 3000
156  CALL mpi_bcast(wa2,n,mpi_real8,iproc_min,
157  1 mpi_comm_world,ierr)
158  IF (ierr .ne. 0) GOTO 3000
159  CALL mpi_bcast(wa4,m,mpi_real8,iproc_min,
160  1 mpi_comm_world,ierr)
161  IF (ierr .ne. 0) GOTO 3000
162 
163 !
164 ! BROADCAST JACOBIAN fjac(:n,j), j=1,n ROW BY ROW (CHANGED IN LMPAR) TO OTHER PROCESSES
165 !
166  DO j = 1, n
167  IF (myid .eq. iproc_min) wa3(:n) = fjac(:n,j)
168  CALL mpi_bcast(wa3, n, mpi_real8, iproc_min,
169  1 mpi_comm_world, ierr)
170  IF (ierr .ne. 0) GOTO 3000
171  IF (myid .ne. iproc_min) fjac(:n,j) = wa3(:n)
172  END DO
173 
174 !
175 ! CLEANUP AFTER LEVENBERG-MARQUARDT LOOP AS NEEDED (WA4 IS NOT CHANGED)
176 !
177  iflag = flag_cleanup
178  CALL fcn (m, n, x, wa4, iflag, nfev) !Contains Bcast Barrier
179 
180  nfev = nfev + num_lev
181 
182  RETURN
183 
184  3000 CONTINUE
185 
186  WRITE (6, *) 'MPI_BCAST error in LEVMARQ_PARAM_MP, ierr = ', ierr
187  stop
188 #endif
189  END SUBROUTINE levmarq_param_mp
190 
191  SUBROUTINE levmarq_param(x, wa1, wa2, wa3, wa4,
192  1 time, nfev, m, n, iflag, fcn)
193  USE fdjac_mod
194  IMPLICIT NONE
195 C-----------------------------------------------
196 C D u m m y A r g u m e n t s
197 C-----------------------------------------------
198  INTEGER :: nfev, n, m, iflag
199  REAL(rprec) :: time
200  REAL(rprec), TARGET :: x(n), wa1(n), wa2(n), wa3(n), wa4(m)
201  EXTERNAL fcn
202 #if !defined(MPI_OPT)
203 C-----------------------------------------------
204 C L o c a l P a r a m e t e r s
205 C-----------------------------------------------
206  REAL(rprec), PARAMETER :: zero=0
207 C-----------------------------------------------
208 C L o c a l V a r i a b l e s
209 C-----------------------------------------------
210  INTEGER :: j, istat, iread, ic1, ic2, irate, count_max,
211  1 jmin
212  REAL(rprec), DIMENSION(num_lm_params) ::
213  1 fnorm_min, pnorm_min, delta_min, par_min
214  CHARACTER(LEN=1) :: ext, low_mark
215 C-----------------------------------------------
216 C E x t e r n a l F u n c t i o n s
217 C-----------------------------------------------
218  EXTERNAL multiprocess, lmpar_parallel
219 C-----------------------------------------------
220 !
221 ! Initialize module variables for use in lmpar_parallel
222 !
223  xp => x
224  wap => wa1
225  wa2p => wa2
226  wa3p => wa3
227  wa4p => wa4
228  np = n
229  mp = m
230  ncntp = nfev
231 
232  CALL system_clock(ic1, irate)
233 
234  CALL multiprocess(num_lm_params, max_processors,
235  1 lmpar_parallel, fcn)
236 
237  CALL system_clock(ic2, irate, count_max)
238  IF (ic2 .lt. ic1) ic2 = ic2 + count_max
239 
240  nfev = nfev + num_lm_params
241 
242 !
243 ! Read in optimal wa1, wa2, wa4, par, delta, pnorm, fnorm1 value from file
244 !
245  DO j = 1, num_lm_params
246 
247  READ (j+1000, iostat=iread) istat, iflag, pnorm_min(j),
248  1 fnorm_min(j), par_min(j), delta_min(j)
249  IF (iread .ne. 0) THEN
250  WRITE (6, *) 'Error reading from file fort.', j+1000,
251  1 ' in levmarq_param', ' IOSTAT = ', iread
252  iflag = -15
253  ELSE IF (j .ne. istat) THEN
254  WRITE (6, *)
255  1 'Incorrect value READ in for INDEX j in levmarq_param'
256  iflag = -15
257  END IF
258 
259  IF (iflag .NE. 0) RETURN
260 
261  IF (j .EQ. 1) fnorm1 = fnorm_min(j)
262  IF (fnorm_min(j) .LE. fnorm1) THEN
263  jmin = j
264  fnorm1 = fnorm_min(jmin)
265  pnorm = pnorm_min(jmin)
266  par = par_min(jmin)
267  delta = delta_min(jmin)
268 #if defined(CRAY)
269  DO k = 1, n
270  READ (j+1000) wa1(k), wa2(k)
271  DO istat = 1, n
272  READ (j+1000) fjac(k, istat)
273  END DO
274  END DO
275  DO k = 1, m
276  READ (j+1000) wa4(k)
277  END DO
278 #else
279  READ (j+1000) wa1, wa2, wa4, fjac(1:n, 1:n)
280 #endif
281 
282  END IF
283 
284  CLOSE (j+1000, status='delete') !!Needed to run correctly in multi-tasking...
285 
286  END DO
287 
288  DO j = 1, num_lm_params
289  ext = ' '
290  low_mark = ' '
291  IF (j .eq. 1) ext = '*'
292  IF (j .eq. jmin) low_mark = '*'
293  WRITE (6, '(2x,i6,4x,2(3x,1es12.4,a),3x,1es12.4)') j+ncntp,
294  1 fnorm_min(j)**2, low_mark, par_min(j), ext, delta_min(j)
295  END DO
296 
297 !
298 ! Do any special cleanup now for IFLAG = flag_cleanup. WA4 LEFT UNCHANGED
299 !
300  iflag = flag_cleanup
301  CALL fcn(m, n, x, wa4, iflag, ncntp)
302 
303  time = time + real(ic2 - ic1)/real(irate) !!Time in multi-process CALL
304 #endif
305  END SUBROUTINE levmarq_param
306 
307  END MODULE lmpar_mod
mpi_inc
Umbrella module avoid multiple inlcudes of the mpif.h header.
Definition: mpi_inc.f:11
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17