SUBROUTINE S8D(N,K,T,D,A,TIMES,AA)
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c This subroutine generates the random data for the small dense matrix
c version of Benchmark #8.
c
c Parameters:
c
c Provided by calling routine:
c N = The size of the A matrices
c K = The number of A matrices
c T = The length of the D array
c
c Returned by this routine:
c D = An array of length T, holding integers mod K
c A = K real NxN matrices of positive real numbers
c TIMES = An NxK long array of positive real numbers
c AA = An NxN matrix of positive real numbers
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
C--CVS variable declaration
TYPE CVS
sequence
character( 160 ) string
integer stringend
END TYPE CVS
C--CVS initilaize variables
TYPE( CVS ),save :: CVS_INFO =
$ CVS("BMARKGRP $Date: 2005/01/07 23:07:11 $ $Revision: 1.2 $" //
$ "$RCSfile: s8d.f,v $ $Name: rel_5 $", 0)
INTEGER T
c
c Which matrix to use in each step
INTEGER D(T)
c
c The matrices
REAL A(N,N,K)
c
c The scaling factor by which to multiply the columns of AA in order to
c obtain the A matrices
REAL TIMES(N,K)
c
c An array from which to generate the A matrices
REAL AA(N*N)
c
c Generate random numbers between 1 and K for the D array
CALL VRAND(0,D,T)
DO 5 I = 1,T
D(I) = MOD(D(I),K) + 1
5 CONTINUE
c
c Initialize a matrix AA from which to generate the A matrices
CALL RANDIV(N,N,AA)
c
c Initialize the scaling factors from which to derive the A matrices from
c the AA matrix
CALL RANDIV(N,K,TIMES)
c
c Generate the A matrices from the AA matrix and the scaling factors
DO 50 I = 1,N
DO 40 J = 1,N
DO 30 KK = 1,K
A(I,J,KK) = AA(I+N*J-N) * TIMES(J,KK)
30 CONTINUE
40 CONTINUE
50 CONTINUE
c
RETURN
END
SUBROUTINE RANDIV(M,N,A)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c Subroutine to produce M N-long lists of random numbers, with wide
c variance in their magnitudes, but whose sum is 1.0.
c
c Each list is done independently from the others, so we will talk
c about only one. The algorithm is to start with an N-long partition,
c whose sum of eventual numbers is known to be 1.0. The partition is
c split into two partitions, the dividing point being random. The sum
c is also split randomly, with one part going to one partition, the
c other part to the other partition. If both partitions contain
c more than one element, one is pushed onto the stack and the
c other is processed further. If either contains only one element,
c its sum is the final value of that element, so it is done and
c we work on the other partition without pushing this one on
c the stack. If both contain only one element, we pop a previously
c pushed partition off the stack and work on it. We then repeat
c the above until the stack is empty and there is no partition
c we are currently working on, which occurs after N-1 steps.
c Finally, we shuffle the positions of the elements.
c
c Notice that we require 3*(N-1) random numbers per list: one for
c splitting the partitions, one for splitting the sums, and one for
c the final shuffle. These random numbers will be computed as needed.
c
c Parameters:
c
c Provided by calling routine:
c M = Number of lists to generate
c N = Length of each list
c
c Returned by this routine:
c A = M by N Array to hold results
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
REAL A(M,N)
c
c Maximum allowed values of M and N
PARAMETER( MM = 100)
c
c The stack of the end points of the partitions
INTEGER STKEND(MM,MM)
c
c The stack of the starting points of the partitions
INTEGER STKST(MM,MM)
c
c The stack pointer (number of partitions found but not current)
INTEGER IPTR(MM)
c
c The starting and end points of the current partition
INTEGER ISTART(MM),IEND(MM)
EXTERNAL IRAND
c
c The sum of the elements in the current partition
REAL RSUM(MM)
c
c An array to hold randum real numbers between 0 and 1
REAL RND(MM, MM)
c
c Check values of M and N against MM
IF((M .GT. MM) .OR. (N .GT. MM)) THEN
PRINT 5,M,N,MM
5 FORMAT('Error in RANDIV. M = ',I5,' N = ',I5,
$ ' Maximum allowed = ',I6)
STOP
ENDIF
c
c Initalize the array of random real numbers
X = .5**32
DO 20 I = 1, N
DO 10 J = 1, M
RND(J,I) = X*IRAND(0)
10 CONTINUE
20 CONTINUE
c
c Initialize the sum to be 1.0, the partition to begin at the first
c element and end at the last, and the stack to be empty
DO 30 K = 1,M
STKEND(K,1) = N
ISTART(K) = 1
IEND(K) = N
IPTR(K) = 1
RSUM(K) = 1.0
30 CONTINUE
c
c We are guaranteed to finish in exactly N-1 steps
DO 50 ISTEP = 1,N-1
DO 40 K = 1,M
c
c Find the length of the current partition
IDIFF = IEND(K) - ISTART(K)
c
c Split the current sum
RAND1 = RSUM(K)*RND(K,ISTEP)
RAND2 = RSUM(K) - RAND1
c
c Default the next current sum to be the one for the first partition
TSUM = RAND1
c
c Push the sum for the other partition onto the stack
RND(K,IPTR(K)+1) = RAND2
c
c Calculate the end of the first partition by splitting the current
c partition
ITEND = (IRAND(0) * X)*IDIFF + ISTART(K)
c
c The start of the second partition is just past the end of the first
STKST(K,IPTR(K)+1) = ITEND + 1
c
c The end of the 2nd partition is the end of the current partition
STKEND(K,IPTR(K)+1) = IEND(K)
c
c The start of the first partition is the start of the current one
ITSTRT = ISTART(K)
c
c Assume both partitions are longer than 1
INC = 1
c
c Save the two sums in case the length of their partition is 1 -- will
c be overwritten later if not
A(K,ISTART(K)) = RAND1
A(K,IEND(K)) = RAND2
c
c Is the first partition of length 1? If so, then ...
IF (ITEND .EQ. ITSTRT) THEN
c
c ... we don't have to push onto the stack
INC = 0
c
c The next current start, end, and sum are those of the second partition
ITSTRT = ITSTRT + 1
ITEND = IEND(K)
TSUM = RAND2
c
c Unless the second partition is also of length 1, in which case ...
IF (IDIFF .EQ. 1) THEN
c
c ... we pop the next current start, end, and sum from the stack
INC = -1
ITSTRT = STKST(K,IPTR(K))
ITEND = STKEND(K,IPTR(K))
TSUM = RND(K,IPTR(K))
END IF
c
c If the first partition is not of length one but the second is, ...
ELSE IF (ITEND .EQ. IEND(K)-1) THEN
c
c ... the next current start, end, and sum are as assumed above, but
c we don't need to push the 2nd partition onto the stack
INC = 0
END IF
c
c Update the current start, end, and sum, as well as the stack pointer
ISTART(K) = ITSTRT
IEND(K) = ITEND
IPTR(K) = IPTR(K) + INC
RSUM(K) = TSUM
40 CONTINUE
50 CONTINUE
c
c Do the final shuffle
DO 70 I = N,2,-1
DO 60 K = 1,M
ITMP = (IRAND(0)*X*I) + 1
TMP = A(K,I)
A(K,I) = A(K,ITMP)
A(K,ITMP) = TMP
60 CONTINUE
70 CONTINUE
RETURN
END