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