program tariff
! THIS PROGRAM SOLVES THE LINEAR ASYMMETRIC SPATIAL PRICE EQUILIBRIUM
! PROBLEM with tariffs using the modified projection method.
! Details of the algorithm and model, along with numerical results
! can be found in the paper, "Spatial price equilibrium models with
! discriminatory ad valorem tariffs: formulation and comparative
! computation using variational inequalities," Nagurney,
! Nicholson, and Bishop.
!
! This program is the initial parallel version implemented for
! the IBM RS/6000 SP using the MPI message passing library. This
! program is part of the Tariff Case Study.
! Read data in by subroutine 'getdata'
! The following parameter controls number of iterations in short runs:
parameter (niter = 100)
integer MAXEXTEN
parameter (MAXEXTEN = 500)
common /a/ nsupm, ndemm, ncomm, ncr
common /a1/ sslope(MAXEXTEN, 5), shxlin(MAXEXTEN), &
& supply(MAXEXTEN), supp(MAXEXTEN)
common /a2/ dslope(MAXEXTEN, 5), dhxlin(MAXEXTEN), &
& demand(MAXEXTEN), demp(MAXEXTEN)
common /a3/ gslope, thxlin(MAXEXTEN, MAXEXTEN)
common /a4/ ind1(MAXEXTEN, MAXEXTEN), ind2(MAXEXTEN, MAXEXTEN)
common /a5/ apatho(MAXEXTEN, MAXEXTEN), tar(MAXEXTEN, MAXEXTEN)
real*8 rtc, ti1, ti2, ti3, t1, t2, t3, to1, to2, to3
! MPI: Added the MPI library header file and
! data declarations to support MPI
! & data partitioning
!
include 'mpif.h'
integer Err, CommSize, Rank
integer Chunk, Extra, Lower, Upper
! MPI: Initialization
!
call mpi_init( Err )
call mpi_comm_size( MPI_COMM_WORLD, CommSize, Err )
call mpi_comm_rank( MPI_COMM_WORLD, Rank, Err )
! READ INPUT DATA
! read in number of supply markets, number of demand markets,
! number of commodities, and number of cross-terms in the price
! functions
ti1 = rtc()
call getdata
ti2 = rtc()
! COMPUTATION BEGINS
t1 = rtc()
npath = nsupm * ndemm
itcnt = 1
iict = 0
rho = .0001
! MPI: Partition the data -
! determine new bounds for demand index j
! using total number of processors &
! rank of *this* processor
!
! Number of columns per processor
Chunk = ndemm / CommSize
! Number of columns "left over"
Extra = mod( ndemm, CommSize )
! Lower bound to skips data assigned
! to lower ranked processors,
! including "extra" columns
Lower = Chunk * Rank + 1 + min( Extra, Rank )
! "extra" processors get 1 more column
if( Rank .lt. Extra ) then
Upper = Lower + Chunk
else
Upper = Lower + Chunk - 1
endif
! CONSTRUCT AN INITIAL FEASIBLE FLOW ASSIGNMENT
! MOST LIKELY CAN BE PARALLELIZED (parallelize on j)
! INITIALIZATIONS
!
! MPI: changed limits on demand index j
!
do 65 j = Lower, Upper
do 70 l = 1, nsupm
apatho(l, j) = 0.
70 continue
65 continue
! MPI: no change here -
! initialize entire demand vector
!
do 651 i=1, ndemm
demand(i) = 0.0
651 continue
do 652 i=1, nsupm
supply(i) = 0.0
652 continue
! (OUTER LOOP STARTS HERE)
! Exchange commenting on the following 2 statements for limited run:
do while (iict .lt. npath) ! production run
! do while (itcnt .lt. NITER) ! limited run
! CONSTRUCT NEW SUPPLY PRICE FUNCTIONS
! COMPUTE SUPPLY PRICES, DEMAND PRICES, AND TRANSPORTATION COSTS
do 1175 j = 1, nsupm
supp(j) = shxlin(j)
do 1176 k = 1, ncr
supp(j) = supp(j) + sslope(j, k) * supply(ind1(j, k))
1176 continue
1175 continue
! MPI: changed limits on demand index j
!
do 1185 j = Lower, Upper
demp(j) = dhxlin(j)
do 1187 k = 1, ncr
demp(j) = demp(j) + dslope(j, k) * demand(ind2(j, k))
1187 continue
1185 continue
do j = 1, nsupm
supply(j) = 0.
enddo
! MPI: reset of demand array moved
! to insure entire array resets
!
do j = 1, ndemm
demand(j) = 0.
enddo
! this is one of the two main iterative steps and can be parallelized
! MPI: changed limits on demand index j;
! removed the reset of demand array
!
do 2225 j = Lower, Upper
do 2245, k = 1, nsupm
test = (-supp(k) - thxlin(k,j)) * (1. + tar(k,j)) + demp(j)
apath = apatho(k,j) + (rho * (test))
if (apath .lt. 0.) apath = 0.
supply(k) = supply(k) + apath
demand(j) = demand(j) + apath
2245 continue
2225 continue
! MPI: accumulate and share supply & demand arrays
!
call mpi_allReduce( supply, supply, nsupm,
& MPI_REAL, MPI_SUM,
& MPI_COMM_WORLD, Err )
call mpi_allReduce( demand, demand, ndemm,
& MPI_REAL, MPI_SUM,
& MPI_COMM_WORLD, Err )
do 11175 j = 1, nsupm
supp(j) = shxlin(j)
do 11176 k = 1, ncr
supp(j) = supp(j) + sslope(j, k) * supply(ind1(j, k))
11176 continue
11175 continue
! MPI: Changed limits on demand index j
!
do 11185 j = Lower, Upper
demp(j) = dhxlin(j)
do 11187 k = 1, ncr
demp(j) = demp(j) + dslope(j, k) * demand(ind2(j, k))
11187 continue
11185 continue
do j = 1, nsupm
supply(j) = 0.
enddo
! MPI: reset of demand array moved
! to insure entire array resets
!
do j = 1, ndemm
demand(j) = 0.
enddo
! this iterative step can also be parallelized
! MPI: changed limits on demand index j;
! removed the reset of demand array
!
iict = 0
do 12225 j = Lower, Upper
do 12230 k = 1, nsupm
test =(-supp(k) - thxlin(k, j))*(1. + tar(k, j)) + demp(j)
apath = apatho(k, j) + rho * test
if (apath .lt. 0.) apath = 0.
if (abs(apath- apatho(k,j)) .le. .01) iict = iict+1
apatho(k, j) = apath
supply(k) = supply(k) + apath
demand(j) = demand(j) + apath
12230 continue
12225 continue
itcnt = itcnt + 1
! MPI: accumulate and share supply & demand arrays
!
call mpi_allReduce( supply, supply, nsupm,
& MPI_REAL, MPI_SUM,
& MPI_COMM_WORLD, Err )
call mpi_allReduce( demand, demand, ndemm,
& MPI_REAL, MPI_SUM,
& MPI_COMM_WORLD, Err )
! MPI: accumulate and share convergence test value
!
call mpi_allReduce( iict, iict, 1, MPI_INTEGER,
& MPI_SUM, MPI_COMM_WORLD, Err )
end do
! (END OF OUTER LOOP)
t2 = rtc()
! OUTPUT SECTION.
to1 = rtc()
write (*, 10) nsupm, ndemm, ncomm, ncr
10 format(' nsupm=',I5,' ndemm=',I5,' ncomm=',I5,' ncr=',I5/)
ti3 = ti2 - ti1
t3 = t2 - t1
write (*, 890) ti3, t3, itcnt, iict
890 format (' DATA INPUT WALLCLOCK TIME IN SECONDS=',f12.4,/ &
& ' COMPUTATION WALLCLOCK TIME IN SECONDS=',f12.4,// &
& ' itcnt=',i7,' iict=',i9)
itc = 0
! MPI: changed limits on demand index j
!
do j = Lower, Upper
do 250 k = 1, nsupm
if (apatho(k, j) .gt. 0.) itc = itc + 1
250 continue
enddo
! MPI: accumulate and share final shipment
! path count
!
call mpi_allReduce( itc, itc, 1, MPI_INTEGER, &
& MPI_SUM, MPI_COMM_WORLD, Err )
perc = (float(itc) / float(npath)) * 100.
write (*, 138) perc, itc
138 format (' perc=', f10.4, 2x, 'itc = ', i5/)
to2 = rtc()
to3 = to2 - to1
! MPI: added call to finalize routine
!
MPI_FINALIZE( Err )
write (*, 895) to3
895 format (' REPORT OUTPUT WALLCLOCK TIME IN SECONDS=',f12.4)
write (*, 899) (to2 - ti1)
899 format (' OVERALL WALLCLOCK TIME IN SECONDS=',f12.4)
stop
end