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
8 LOGICAL,
DIMENSION(:),
ALLOCATABLE :: flip
12 SUBROUTINE fdjac2_mp(fcn, m, n, x, fvec, fjac, ldfjac,
13 1 iflag, ncnt, epsfcn, fnorm_min, x_min, wa)
21 INTEGER,
INTENT(in) :: m, n, ldfjac, ncnt
23 REAL(rprec),
INTENT(in) :: epsfcn
25 REAL(rprec),
INTENT(in) :: fvec(m)
26 REAL(rprec) :: fjac(ldfjac,n)
27 REAL(rprec) :: fnorm_min, x_min(n), wa(m)
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,
39 REAL(rprec),
ALLOCATABLE :: fjac_overflow(:,:), fnorm_overflow(:)
40 INTEGER :: istat, ierr_flag, ikey
42 LOGICAL,
ALLOCATABLE :: flip_overflow(:)
46 EXTERNAL dpmpar, enorm
122 ALLOCATE (fnorm_array(n), h(n), buffer(m), stat=istat)
123 IF (istat .ne. 0) stop
'Allocation error in fdjac2_mp'
129 eps = sqrt(max(epsfcn,epsmch))
130 cur_norm = enorm(m,fvec)
132 h(1:n) = eps*abs(x(1:n))
133 WHERE (h .eq. zero) h=eps
140 num_split = n/numprocs
141 mpi_comm_workers = mpi_comm_world
150 CALL fcn(m, n, x, wa, iflag, ncnt)
153 buffer(:m) = (wa(:m) - fvec(:m))/h(j)
159 IF (temp > cur_norm) lflip = .not. flip(j)
162 WRITE (6,
'(2x,i6,8x,i3,7x,1es12.4)') ncnt+j,
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)
173 joff = joff + numprocs
175 IF (ierr_mpi .ne. 0)
GO TO 100
182 num_split = mod(n, numprocs)
183 IF (num_split .ne. 0)
THEN
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
195 ikey = worker_split_key+1
202 CALL mpi_comm_split(mpi_comm_world, ikey, worker_id,
203 1 mpi_comm_workers, ierr_mpi)
205 1 stop
'MPI_COMM_SPLIT error in fdjac_mp_queue'
210 CALL fcn(m, n, x, wa, iflag, ncnt)
212 buffer(:m) = (wa(:m) - fvec(:m))/h(j)
214 IF (temp > cur_norm) lflip = .not. flip(j)
216 WRITE (6,
'(2x,i6,8x,i3,7x,1es12.4)') ncnt+j,
219 CALL mpi_comm_free(mpi_comm_workers, ierr_mpi)
221 buffer = 0; lflip = .true.
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)
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)
235 DEALLOCATE(fjac_overflow, fnorm_overflow, flip_overflow)
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'
249 wa(:) = fjac(:, iproc_min)*h(iproc_min) + fvec(:)
251 x_min(iproc_min) = x(iproc_min) + h(iproc_min)
253 DEALLOCATE (h, fnorm_array, buffer)
259 mpi_comm_workers = mpi_comm_world
260 CALL fcn(m, n, x, wa, ikey, ncnt)
265 WRITE (6, *)
' MPI_BCAST error in FDJAC2_MP: IERR=', ierr_mpi,
266 1
' PROCESSOR: ',myid
268 END SUBROUTINE fdjac2_mp
271 SUBROUTINE fdjac2_mp_queue(fcn, m, n, x, fvec, fjac, ldfjac,
272 1 iflag, ncnt, epsfcn, fnorm_min, x_min, wa)
280 INTEGER,
INTENT(in) :: m, n, ldfjac, ncnt
282 REAL(rprec),
INTENT(in) :: epsfcn
284 REAL(rprec),
INTENT(in) :: fvec(m)
285 REAL(rprec) :: fjac(ldfjac,n)
286 REAL(rprec) :: fnorm_min, x_min(n), wa(m)
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,
298 INTEGER :: status(MPI_STATUS_size)
299 INTEGER :: numsent, sender, istat, ikey
300 INTEGER :: anstype, column, ierr_flag
304 EXTERNAL dpmpar, enorm
380 ALLOCATE (buffer(m), fnorm_array(n), h(n), stat=istat)
381 IF (istat .ne. 0) stop
'Allocation error in fdjac2_mp'
388 eps = sqrt(max(epsfcn,epsmch))
391 IF (myid .eq. master)
THEN
394 ikey = worker_split_key+1
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'
412 IF (myid .eq. master)
THEN
415 cur_norm = enorm(m,fvec)
419 h(1:n) = eps*abs(x(1:n))
420 WHERE (h .eq. zero) h=eps
426 DO j = 1,min(numprocs-1,n)
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'
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)
447 IF (anstype .gt. n) stop
'ANSTYPE > N IN FDJAC2'
449 fjac(:m,anstype) = (wa(:m) - fvec(:m))/h(anstype)
454 fnorm_array(anstype) = temp
455 IF (temp > cur_norm) flip(anstype) = .not. flip(anstype)
458 WRITE (6,
'(2x,i6,8x,i3,7x,1es12.4)') ncnt+anstype,
465 IF (numsent .lt. n)
THEN
468 x(numsent) = temp + h(numsent)
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'
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'
486 ELSE IF (myid .le. n)
THEN
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'
497 column = status(mpi_tag)
498 IF (column .ne. 0)
THEN
502 CALL fcn(m, n, x, wa, iflag, ncnt)
504 IF (iflag.ne.0 .and. ierr_flag.eq.0) ierr_flag = iflag
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'
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)
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'
538 wa(:) = fjac(:, iproc_min)*h(iproc_min) + fvec(:)
540 x_min(iproc_min) = x(iproc_min) + h(iproc_min)
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,
550 IF (ierr_mpi .ne. 0)
GOTO 100
565 DEALLOCATE (h, buffer, fnorm_array)
570 column = flag_cleanup
571 CALL fcn(m, n, x, wa, column, ncnt)
577 CALL mpi_bcast(x, n, mpi_real8, master, mpi_comm_world, ierr_mpi)
578 IF (ierr_mpi .ne. 0)
GOTO 100
580 IF (ierr_flag .ne. 0)
THEN
582 CALL mpi_bcast(iflag, 1, mpi_integer, myid,
583 1 mpi_comm_world, ierr_mpi)
584 IF (ierr_mpi .ne. 0)
GOTO 100
587 IF (myid .ne. master)
588 1
CALL mpi_comm_free(mpi_comm_workers, ierr_mpi)
593 WRITE (6, *)
' MPI_BCAST error in FDJAC2_MP: IERR=', ierr_mpi,
594 1
' PROCESSOR: ',myid
596 END SUBROUTINE fdjac2_mp_queue
599 SUBROUTINE fdjac2(fcn, m, n, x, fvec, fjac, ldfjac, iflag, ncnt,
600 1 epsfcn, wa, time, fnorm_min, x_min, fvec_min)
602 #if !defined(MPI_OPT)
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)
617 REAL(rprec),
PARAMETER :: zero = 0
621 INTEGER :: j, istat, iread, ic1, ic2, irate, count_max
622 REAL(rprec) :: epsmch, h, dpmpar, temp, cur_norm
626 EXTERNAL fcn, dpmpar, multiprocess, fdjac_parallel
627 REAL(rprec),
EXTERNAL :: enorm
690 eps = sqrt(max(epsfcn,epsmch))
704 fnorm_min = huge(fnorm_min)
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
723 cur_norm = enorm(m,fvec)
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
732 ELSE IF (j .ne. istat)
THEN
733 WRITE (6, *)
'Wrong value for INDEX j READ in fdjac2'
737 IF (iflag .ne. 0)
EXIT
746 fjac(:m,j) = (wa - fvec)/h
748 IF( temp > cur_norm) flip(j) = .not. flip(j)
750 IF( temp < fnorm_min)
THEN
755 READ (j+1000) x_min(k)
762 CLOSE (j+1000, status=
'delete')
770 CALL fcn(m, n, x, wa, iflag, ncnt)
772 time = time + real(ic2 - ic1)/real(irate)
774 END SUBROUTINE fdjac2