PROGRAM ibquad ! ! Use DOPT to find good incomplete-block designs for fitting ! quadratic surfaces. ! USE d_optimal_design IMPLICIT NONE INTEGER, PARAMETER :: maxf=6, maxb=10, maxcol=maxf*(maxf+3)/2, & kmax=maxcol+maxb, nrmax=kmax*(kmax-1)/2, & mxcand=3**maxf, nmax=100 INTEGER :: nfact, nblock, blksiz(maxb), n, kin, kfull, nrbar, & ncand, picked(nmax), ifault, point, i, j, ind(maxcol), & pos, nrep, design(nmax), i1, i2 REAL (dp) :: d(kmax), rbar(nrmax), lndet, xx(kmax), & tol(kmax), zpz(mxcand,maxb), wk(kmax), dbest, stdet REAL :: x(mxcand,maxcol) LOGICAL, PARAMETER :: rstart = .TRUE. INTEGER, PARAMETER :: in(maxb) = (/ (0,i=1,maxb) /) INTEGER, ALLOCATABLE :: seed(:) CHARACTER (LEN=1) :: bel REAL (dp), PARAMETER :: zero = 0.0_dp WRITE(*, *)' +++ Incomplete blocks for quad. surfaces +++' bel = CHAR(7) 10 WRITE(*, *) 'Enter no. of factors: ' READ(*, *) nfact IF (nfact > maxf .OR. nfact < 2) THEN WRITE(*, 900) bel 900 FORMAT(' ', a, '** Illegal value entered **') GO TO 10 END IF WRITE(*, *) 'Enter no. of blocks: ' READ(*, *) nblock IF (nblock > maxb .OR. nblock < 0) THEN WRITE(*, 900) bel GO TO 10 END IF WRITE(*, *)'Enter size of each block: ' READ(*, *) blksiz(1:nblock) n = SUM( blksiz(1:nblock) ) IF (n > nmax) THEN WRITE(*, *) bel, '** Sorry, design is too large **' GO TO 10 END IF WRITE(*, *)'How many tries?: ' READ(*, *) nrep CALL RANDOM_SEED(size=i) ALLOCATE( seed(i) ) WRITE(*, *)'Enter', i, ' integers for random no. seeds: ' READ(*, *) seed CALL RANDOM_SEED(put=seed) kin = nfact * (nfact + 3) / 2 kfull = kin + nblock IF (n < kfull) THEN WRITE(*, *) bel, '** Design too small to fit model **' GO TO 10 END IF nrbar = kfull * (kfull - 1) / 2 ncand = 3**nfact ! Open a file for the output ! Default file name: DESIGN.OPT OPEN(10, FILE='DESIGN.OPT') WRITE(10, 1000) nfact 1000 FORMAT('Output from DOPT for fitting quadratic surfaces'/ & i2, ' factors at 2 levels') IF (nblock .GT. 1) THEN WRITE(10, '(a, 10i4)') 'Block sizes: ', blksiz ELSE WRITE(10, '(a, i5)') 'No. of experimental runs = ', n END IF ! Generate the candidate design points. ! ind(i) stores the level (-1,0,1) of the i-th factor. point = 1 ind(1:nfact) = -1 GO TO 55 30 i = 1 40 ind(i) = ind(i) + 1 IF (ind(i) > 1) THEN IF (i == nfact) GO TO 90 ind(1:i) = -1 i = i + 1 GO TO 40 END IF ! ! Calculate row of X. ! 55 x(point,1:nfact) = ind(1:nfact) pos = nfact + 1 DO i = 1, nfact DO j = 1, i x(point,pos) = ind(i) * ind(j) pos = pos + 1 END DO END DO point = point + 1 GO TO 30 90 dbest = zero DO i = 1, nrep CALL dopt(x, mxcand, ncand, kin, n, nblock, in, blksiz, kfull, & rstart, nrbar, d, rbar, picked, lndet, xx, tol, zpz, wk, & ifault) IF (ifault /= 0) THEN WRITE(*, *) bel, 'IFAULT = ', ifault ELSE WRITE(*, *) i, ' Log Det. = ', lndet END IF IF (lndet > dbest) THEN design(1:n) = picked(1:n) dbest = lndet END IF END DO WRITE(*, *) stdet = EXP(dbest - kfull*LOG(REAL(n, KIND=dp))) WRITE(*, *) 'Max. log det. = ', dbest, ' Std. det = ', stdet WRITE(*, *) 'Design:' WRITE(10, *) 'Max. log det. = ', dbest, ' Std. det = ', stdet WRITE(10, *) 'Design:' i2 = 0 DO i = 1, nblock i1 = i2 + 1 i2 = i2 + blksiz(i) WRITE(*, 960) i, (j, j=1,nfact) WRITE(10, 960) i, (j, j=1,nfact) 960 FORMAT(' BLOCK', i5 / t15, 'F A C T O R' / ' ', 15I5) DO j = i1, i2 WRITE(*, 970) x(design(j), 1:nfact) WRITE(10, 970) x(design(j), 1:nfact) 970 FORMAT(' ', 15F5.0) END DO END DO WRITE(*, *) GO TO 10 STOP END PROGRAM ibquad