1C***********************************************************************
2C LHS (Latin Hypercube Sampling) UNIX Library/Standalone.
3C Copyright (c) 2004, Sandia Corporation.  Under the terms of Contract
4C DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
5C retains certain rights in this software.
6C
7C This software is distributed under the GNU Lesser General Public License.
8C For more information, see the README file in the LHS directory.
9C***********************************************************************
10C     Last change:  SLD  11 Jul 101   11:15 am
11!LHS_EXPORT_DEC ATTRIBUTES DLLEXPORT::LHS_RUN
12      SUBROUTINE LHS_RUN(MAXVAR,MAXOBS,MAXNAM,IError,
13     x LSTNAME,INDXNAM,PTVALST,NUMNAM,SMATX,NUMVAR,RMATX,RFLAG)
14c     LHS_RUN runs the LHS sample engine
15c     LHS_RUN inputs:
16c        MAXVAR - Maximum number of variables, integer
17c        MAXOBS - Maximum number of observation, integer
18c        MAXNAM - Maximum number of variable names including "same as"
19c        RFLAG - flag to indicate whether the rank matrix is being passed in
20c     LHS_RUN returns:
21c        LSTDNAM - a list of distribution names (type character)
22c        INDXNAM - integer array containing index number(position)
23c                  of names in sample data
24c        PTVALST - an array of the associated point values (type real)
25c        SMATX - the sample data matrix dimension (MAXVAR,MAXOBS)
26c        RFLAG - flag to indicate if no rank data needs to be passed (RFLAG=0),
27c                if rank data will be set (RFLAG = 1),
28c                if rank data will be obtained (RFLAG = 2),
29c                or if rank data will both be "set" and "get" (RFLAG=3)
30c        RMATX - the rank data matrix dimension (MAXVAR,MAXOBS)
31c        NUMNAM - the integer number of variable names including "same as"
32c        NUMVAR- the number of variables for which observations
33c                are generated (with no duplication for same as)
34c        IError - integer error flag
35c           If an error is found, IError = 1 is returned
36c
37      USE KILLFILE
38      USE PARMS
39c     PARAMS provides:   NMAX
40      USE CPARAM
41c     CPARAMS provides:  ISeed,N,NV,NREP,IVarNm,PValue,IDIST,List
42      USE InByCall
43cc    InByCall provides: LINT,LPREP,LFILES,LDIST,NNames
44      USE CSAMP
45c     CSAMP provides: X and XSAVE arrays
46      USE CRANK
47c     CRANK provides: XV array
48
49      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
50      DOUBLE PRECISION :: SMATX(MAXVAR,MAXOBS)
51      DOUBLE PRECISION :: RMATX(MAXVAR,MAXOBS)
52      DOUBLE PRECISION :: PTVALST(MAXNAM)
53      INTEGER :: INDXNAM(MAXNAM)
54      CHARACTER(LEN=16) :: LSTNAME(MAXNAM)
55      INTEGER :: RFLAG
56
57C====================
58      CHARACTER*11 TIMRES
59C====================
60c     next lines added 1-17-96 to convert date/time routine to LF90
61      integer :: dt(8)
62      character (len=5) :: zone
63      character * 8 DATRES
64      character * 2 Chr2
65C
66C     -- STATEMENT FUNCTION
67      LOC(I,J)=(J-1)*N+I
68      NUMNAM = NNames
69      NUMVAR = NV
70c
71c     Check that LHS_INIT or LHS_INIT_MEM, at least one LHS_xDIST,
72c        and LHS_PREP have been called prior calling LHS_RUN
73      I3 = LINIT + LDIST + LPREP
74      IF (I3 /= 3) THEN
75         KLLERR = .TRUE.
76         IError = 1
77         WRITE(*,9022)
78         WRITE(4,9022)
79         WRITE(99,9022)
80         Return
81      END IF
82
83c
84c     Following section from Standard LHS.for routine
85C
86C     -- If only constant and same as distribution types were found,
87C     -- write the output file here and terminate gracefully without
88C     -- sampling.
89C
90      If (NV == 0) Then
91C        -- Write the sample to Unit 1
92         Call SamOut(1)
93cc         If(KLLERR) Return SamOut does not have any error conditions  sld01
94         Write (4,*) ' '
95         Write (4,*) '   Only Constant and Same As distribution ',
96     1      'types were found in the LHS input.'
97         Write (4,*) '   Thus, no sampling is being performed ',
98     1      'by this LHS program run.'
99         Write (4,*) '   The requested constants were written into ',
100     1      'the LHS output file'
101         Write (4,*) '   for processing by other codes.'
102         Write (4,*)
103CCCC         Call FileOC(0)       sld: this is done in LHS_CLOSE
104cc         If(KLLERR) Return     FILEOC has no error conditions         sld01
105ccc         KLLERR = .TRUE.
106ccc
107ccc      IError = 1
108csld asked Greg if this really should be treated as an error
109ccc---sld, Answer was no.
110         RETURN
111      End If
112c
113c     End of section from standard LHS routine
114c
115c
116c
117c     Following section from standard LHS routine
118C
119C     -- THIS LOOP IS EXECUTED ONCE FOR EACH REPETITION REQUESTED.
120C     -- UNITS 7, 8, AND 9 HAVE BEEN DEFINED IN SUBROUTINE RDPAR
121C
122      DO 1000 IREP=1,NREP
123         REWIND 7
124         REWIND 8
125         REWIND 9
126c         IF (IREP .GT. 1) CALL BANNER(IREP)
127         IF (IREP .GT. 1) THEN
128            CALL BANNER(IREP)
129cc            If(KLLERR) Return    BANNER has no error conditions       sld01
130         End If
131C
132C        -- GENERATE THE DISTRIBUTION REQUESTED FOR VARIABLE J
133         DO 100 J=1,NV
134            IDT=IDIST(J)
135            IF(IDT.EQ.1)THEN
136               CALL BETA(J)
137               If(KLLERR) Return
138            ELSE IF (IDT .EQ. 2  .OR.  IDT .EQ. 3  .OR.
139     1               (IDT .GE. 27  .AND.  IDT .LE. 35) ) THEN
140               CALL NORMAL(J,IDT)
141               If(KLLERR) Return
142            ELSE IF(IDT.GE.4.AND.IDT.LE.7)THEN
143               CALL UNIFRM(J,IDT)
144               If(KLLERR) Return
145            ELSE IF(IDT.EQ.8)THEN
146               CALL TRIANG(J)
147               If(KLLERR) Return
148            ELSE IF(IDT .GE. 9  .AND.  IDT .LE. 11) THEN
149               CALL CUMULC(J,IDT)
150               If(KLLERR) Return
151            ELSE IF (IDT .EQ. 12  .OR.  IDT .EQ. 13) THEN
152               CALL CUMULD(J)
153               If(KLLERR) Return
154            ELSE IF (IDT .EQ. 14) THEN
155               CALL POISON(J)
156               If(KLLERR) Return
157            ELSE IF (IDT .EQ. 15) THEN
158               CALL GEOM(J)
159               If(KLLERR) Return
160            ELSE IF (IDT .EQ. 16) THEN
161               CALL BINOM(J)
162               If(KLLERR) Return
163            ELSE IF (IDT .EQ. 17) THEN
164               CALL NBINOM(J)
165               If(KLLERR) Return
166            ELSE IF (IDT .EQ. 18) THEN
167               CALL HGEOM(J)
168               If(KLLERR) Return
169            ELSE IF (IDT .EQ. 19  .OR.  IDT .EQ. 25  .OR.
170     1               IDT .EQ. 26) THEN
171               CALL EXPON(J,IDT)
172               If(KLLERR) Return
173            ELSE IF (IDT .EQ. 20) THEN
174               CALL WEIBUL(J)
175               If(KLLERR) Return
176            ELSE IF (IDT .EQ. 21) THEN
177               CALL PARETO(J)
178                     If(KLLERR) Return
179            ELSE IF (IDT .EQ. 22) THEN
180               CALL GAMMA(J)
181               If(KLLERR) Return
182            ELSE IF (IDT .EQ. 23) THEN
183               CALL IGAUS(J)
184               If(KLLERR) Return
185            ELSE IF (IDT .EQ. 24) THEN
186               CALL ENTRPY(J)
187               If(KLLERR) Return
188            ELSE IF (IDT .EQ. 36) THEN
189               CALL GUMBEL(J)
190               If(KLLERR) Return
191            ELSE IF (IDT .EQ. 37) THEN
192               CALL FRECHET(J)
193               If(KLLERR) Return
194            ENDIF
195  100    CONTINUE
196C
197C        -- IF A RANDOM SAMPLE HAS BEEN GENERATED IT MUST BE SORTED FROM
198C        -- FROM SMALLEST TO LARGEST ON EACH VARIABLE AS PART OF THE
199C        -- STRUCTURING TO REDUCE THE POSSIBILITY OF SPURIOUS CORRELATIONS
200C====================
201      Call date_and_time (DATRES,TIMRES,zone,dt)
202      Chr2 = TIMRES(8:9)
203C
204      If (IPrint > 0 ) Write (*,9010) 'END OF SAMPLING TIME:  ',
205     1                                dt(5), dt(6), dt(7), Chr2
206
207C====================
208         IF (N .NE. 1) THEN
209            IF (IRS .EQ. 1) THEN
210               DO 250 I=1,NV
211                  DO 240 J=1,N
212                     XV(J)=X(LOC(J,I))
213  240             CONTINUE
214                  CALL SIFT (XV,N)
215cc                 If(KLLERR) Return      SIFT has no error conditions  sld01
216                  DO 245 J=1,N
217                     X(LOC(J,I))=XV(J)
218  245             CONTINUE
219  250          CONTINUE
220            END IF
221            NNV=N*NV
222            DO 260 III=1, NNV
223               XSAVE(III)=X(III)
224  260       CONTINUE
225C
226C           -- SUBROUTINE MIX PAIRS THE SAMPLE OBSERVATIONS TO MATCH THE
227C           -- DESIRED CORRELATION STRUCTURE
228C====================
229            Call date_and_time (DATRES,TIMRES,zone,dt)
230            Chr2 = TIMRES(8:9)
231C
232            If (IPrint > 0 ) Write (*,9010) 'END OF SIFT TIME:      ',
233     1                                dt(5), dt(6), dt(7), Chr2
234
235cc          This is used in most cases when the ranks are not pre-specified
236            IF ((RFLAG .EQ. 0).OR.(RFLAG .EQ. 2)) THEN
237            CALL MIX
238            END IF
239cc          This is used in the case where one is sampling incrementally
240cc          and has predefined ranks
241            IF ((RFLAG .EQ. 1).OR.(RFLAG .EQ. 3)) THEN
242cc            WRITE (*,*) 'We are in the loop'
243
244              DO 280 J=1,NV
245                DO 270 I=1,N
246                  RXV(I)=RMATX(J,I)
247                  X(LOC(I,J))=RXV(I)
248cc                WRITE (*,*) 'RXV(I) ', RXV(I), 'I ', I, 'J ', J
249  270           CONTINUE
250  280         CONTINUE
251              DO 290 J=1,NV
252                DO 285 I=1,N
253                  K=INT(X(LOC(I,J))+0.01)
254                  X(LOC(I,J))=XSAVE(LOC(K,J))
255  285           CONTINUE
256  290          CONTINUE
257            END IF
258
259            If(KLLERR) Return
260C
261            Call date_and_time (DATRES,TIMRES,zone,dt)
262            Chr2 = TIMRES(8:9)
263C
264            If (IPrint > 0 ) Write (*,9010) 'END OF MIX TIME:       ',
265     1                                dt(5), dt(6), dt(7), Chr2
266C====================
267         END IF
268C
269ccc
270ccc    Following section modified from LHS.for
271ccc    (i.e. SAMSTOR call added and SamOut only called
272ccc          for user defined output file)
273ccc    Call SAMSTOR to store matrix, point values, variable names and
274ccc    positions in matrix for return to user
275       CALL SamStor(MAXVAR,MAXOBS,MAXNAM,LSTNAME,INDXNAM,
276     x              PTVALST,SMATX)
277ccc
278ccc only call SamOUt if unit 1 is user defined output file
279      IF (IScrh1 == 2) THEN
280C        -- THE SAMPLE IS WRITTEN OUT TO UNIT 1
281         Call SamOut(IRep)
282         If(KLLERR) Return
283      END IF
284ccc   End of modifications
285c
286c        -- Prepare the XSave matrix to reflect the post-MIX values
287         NNV=N*NV
288         DO 300 III=1, NNV
289            XSAVE(III)=X(III)
290  300    CONTINUE
291C
292C        -- SUBROUTINE DATOUT OUTPUTS THE LATEST SAMPLE AND CORRESPONDING
293C        -- RANKS.  SUBROUTINE HSTOUT OUTPUTS HISTOGRAMS FOR THE CURRENT
294C        -- SAMPLE.  SUBROUTINE COROUT OUTPUTS RAW AND RANK CORRELATIONS
295C        -- FOR THE CURRENT SAMPLE.
296c         IF (IDATA .NE. 0) CALL DATOUT
297         IF (IDATA .NE. 0) THEN
298            CALL DATOUT
299            If(KLLERR) Return
300         End If
301c         IF (IHIST .NE. 0) CALL HSTOUT
302         IF (IHIST .NE. 0) THEN
303            CALL HSTOUT
304            If(KLLERR) Return
305         End If
306c         IF (ICORR .NE. 0) CALL COROUT
307         IF (ICORR .NE. 0) THEN
308            CALL COROUT
309            If(KLLERR) Return
310         End IF
311
312         IF ((RFLAG .EQ. 2).OR.(RFLAG.EQ.3)) THEN
313            DO 620 J=1,NV
314               DO 610 I=1,N
315 610              XV(I)=X(LOC(I,J))
316                  CALL RANKER
317               DO 620 I=1,N
318                  RMATX(J,I)=RXV(I)
319 620        CONTINUE
320         End IF
321
322C====================
323      Call date_and_time (DATRES,TIMRES,zone,dt)
324      Chr2 = TIMRES(8:9)
325c
326      If (IPrint > 0 ) Write (*,9010) 'END TIME:              ',
327     1                                dt(5), dt(6), dt(7), Chr2
328C====================
329C
330 1000 CONTINUE
331cccc  End of standard LHS section
332c
333c     re-initialize subroutine called tracking switches
334c     except for LRUN in case LHS_COROUT is called
335      LINIT = 0
336      LPREP = 0
337      LRUN = 1
338      LFILES = 0
339      LDIST = 0
340      IScrh1 = 0
341      IScrh6 = 0
342c
343      Return
344c
345 9010 FORMAT (5x,A,I2,":",I2,":",I2,".",A)
346
347 9022 FORMAT('1',5X,' LHS_INIT or LHS_INIT_MEM,'/,5X,
348     x' at least one distribution call ',/,8X,
349     x'(i.e. LHS_DIST, LHS_UDIST or LHS_SDIST) and',/,5X,
350     x'LHS_PREP must be called prior to calling LHS_RUN')
351      END SUBROUTINE
352