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