C file: USRLC
C
C this file contains the test routines supplied with the microscope
C package. lower case letters are assumed to be available. if only
C upper case letters are available use the file usruc instead of this
C one. for details see the MICROSCOPE manual.
SUBROUTINE SUBUSR
INTEGER IABS
DOUBLE PRECISION AA, OPA
INTEGER K, KK, NM1, NOMAX
INTEGER IOUT
LOGICAL OK, LOGCNT, FIRST
INTEGER OUTPUT, LINES, WIDTH, ILP
INTEGER IDSPLA, IPRMPT
LOGICAL LSCRN
COMMON / SCREEN / OUTPUT, LINES, WIDTH, ILP
COMMON / SCREEN / IDSPLA, IPRMPT, LSCRN
INTEGER INPUTD, GRAPHD, HELPD
INTEGER RECORD, RSTRTD
COMMON / IO / INPUTD, GRAPHD, HELPD
COMMON / IO / RECORD, RSTRTD
INTEGER LCHN
COMMON / LOG / LCHN
DOUBLE PRECISION ETA
INTEGER IROUND, N
LOGICAL ADD
COMMON / USER / ETA, IROUND, N, ADD
C NOMAX IS THE CURRENT MAXIMUM NUMBER OF VERSIONS OF THE TRIAL FUNCTION
DATA NOMAX /15/
C MARK THAT SCREEN DISPLAY HAS BEEN DESTROYED:
LSCRN = .FALSE.
C BEGIN OF MAIN LOOP:
100 CONTINUE
CALL BLSCRN(OUTPUT)
CALL PCURSR(OUTPUT,1,1)
C PREPARE LOGGING:
IOUT = OUTPUT
LOGCNT = .FALSE.
200 CONTINUE
WRITE (IOUT,28000)
IF (N.GE.1) GO TO 300
WRITE (IOUT,30000)
GO TO 500
300 CONTINUE
IF (N.GT.6) GO TO 400
OPA = 1.0D0 + ETA
NM1 = N-1
IF (.NOT.ADD) WRITE (IOUT,10000) NM1
IF (ADD) WRITE (IOUT,12000) NM1
GO TO 500
400 CONTINUE
IF (N.EQ.7) WRITE (IOUT,32000)
IF (N.EQ.8) WRITE (IOUT,22000)
IF (N.EQ.9) WRITE (IOUT,24000)
IF (N.EQ.10) WRITE (IOUT,14000)
IF (N.EQ.11) WRITE (IOUT,16000)
IF (N.EQ.12) WRITE (IOUT,30000)
IF (N.EQ.13) WRITE (IOUT,18000)
IF (N.EQ.14) WRITE (IOUT,20000)
IF (N.EQ.15) WRITE (IOUT,20000)
GO TO 500
500 CONTINUE
WRITE (IOUT,34000) N,ETA
IF (IROUND.GT.0) WRITE (IOUT,36000) IROUND
IF (IROUND.EQ.0) WRITE (IOUT,38000)
600 CONTINUE
IF (LOGCNT) GO TO 1900
WRITE (OUTPUT,40000)
700 CONTINUE
CALL SIREAD(INPUTD,K,OK)
IF (OK) GO TO 800
WRITE (OUTPUT,42000)
GO TO 700
800 CONTINUE
IF (K.GT.(-4)) GO TO 900
WRITE (OUTPUT,44000)
GO TO 600
900 CONTINUE
IF (K.GT.(-3)) GO TO 1200
1000 CONTINUE
WRITE (OUTPUT,46000)
CALL SRREAD(INPUTD,AA,OK)
IF (OK) GO TO 1100
WRITE (OUTPUT,42000)
GO TO 1000
1100 CONTINUE
ETA = AA
GO TO 100
1200 CONTINUE
IF (K.GT.(-2)) GO TO 1300
ADD = .NOT.ADD
GO TO 100
1300 CONTINUE
IF (K.GT.(-1)) GO TO 1600
WRITE (OUTPUT,48000)
1400 CONTINUE
CALL SIREAD(INPUTD,KK,OK)
IF (OK) GO TO 1500
WRITE (OUTPUT,42000)
GO TO 1400
1500 CONTINUE
IROUND = IABS(KK)
GO TO 100
1600 CONTINUE
IF (K.GT.0) GO TO 1700
WRITE (OUTPUT,50000)
IF (LCHN.EQ.0) GO TO 1900
LOGCNT = .TRUE.
IOUT = LCHN
WRITE (IOUT,26000)
GO TO 200
1700 CONTINUE
IF (K.GT.NOMAX) GO TO 1800
N = K
GO TO 100
1800 CONTINUE
WRITE (OUTPUT,52000) NOMAX
GO TO 600
1900 CONTINUE
RETURN
10000 FORMAT(33H The current trial function is s(,I1,7H,eta,x))
12000 FORMAT(33H The current trial function is s(,I1,16H,eta,x) + exp(x)
X)
14000 FORMAT(47H The current trial function is the cubic spline/
X 31H interpolant of the exponential)
16000 FORMAT(47H The current trial function is the error in the/
X 44H cubic spline interpolant of the exponential)
18000 FORMAT(56H The current trial function is the linear function eta*x
X)
20000 FORMAT(42H The current trial function is abs(x)**eta)
22000 FORMAT(30H The current trial function is/
X44H f(x,y) = eta*abs(x)*x**2 + (1-eta)*abs(y)*y)
24000 FORMAT(30H The current trial function is/
X 30H f(x,y) = x**2*abs(x)*y*abs(y))
26000 FORMAT(49H Leaving User Subroutine with parameters set to: )
28000 FORMAT(42H Example for user intervention subroutine,/
X 50H this version supplied with the MICROSCOPE package/
X 38H for details see the MICROSCOPE manual/)
30000 FORMAT(//48H The current trial function is the zero function)
32000 FORMAT(48H The current trial function is f(x) = exp(eta*x))
34000 FORMAT(//37H The current values of N and eta are ,I3,5H and ,
X 1PD12.4)
36000 FORMAT(/13H Rounding to ,I2,26H digits is being simulated)
38000 FORMAT(/31H No rounding is being simulated)
40000 FORMAT(/55H Type a positive integer to select a new trial function
X/32H 0 to leave options as they are
X/38H -1 to change rounding characteristics
X/48H -2 to change the addition of a term (0 < N < 7)
X/30H -3 to change the value of eta/)
42000 FORMAT(34H Number not recognized - try again)
44000 FORMAT(23H Number must be .GE. -3)
46000 FORMAT(23H type new value of eta:)
48000 FORMAT(57H Give number of digits to which results shall be rounded
X,/26H or 0 to turn off rounding)
50000 FORMAT(//27H to obtain a screen display,
X 30H use the GO, RS, or FO command)
52000 FORMAT(29H Version number must be .LT. ,I3)
END
DOUBLE PRECISION FUNCTION F(X)
DOUBLE PRECISION PHI, SPLN, ROUND, DEXP
DOUBLE PRECISION DABS
DOUBLE PRECISION X(1), XX, YY
DOUBLE PRECISION ETA
INTEGER IROUND, N
LOGICAL ADD
COMMON / USER / ETA, IROUND, N, ADD
C
XX = X(1)
YY = X(2)
IF (N.GT.0) GO TO 100
F = 0.0D0
GO TO 400
100 CONTINUE
IF (N.GT.6) GO TO 200
F = 0.0D0
IF (XX.GE.0.0D0) F = ETA*PHI(N,X)
IF (.NOT.ADD) GO TO 400
F = F + DEXP(XX)
GO TO 400
200 CONTINUE
IF (N.EQ.7) F = DEXP(ETA*XX)
IF (N.EQ.8) F = ETA*DABS(XX)*XX**2+(1.0D0-ETA)*DABS(YY)*YY
IF (N.EQ.9) F = XX**2*DABS(XX)*YY*DABS(YY)
IF (N.EQ.10) F = SPLN(XX)
IF (N.EQ.11) F = DEXP(XX) - SPLN(XX)
IF (N.EQ.12) F = 0.0D0
IF (N.EQ.13) F = ETA*XX
IF (N.NE.14) GO TO 300
F = 0.0D0
IF (XX.EQ.0.0D0) GO TO 400
F = DABS(XX)**ETA
IF (XX.LT.0.0D0) F = F*(-1.0D0)**(1+IFIX(SNGL(ETA)))
300 CONTINUE
IF (N.NE.15) GO TO 400
F = 0.0D0
IF (XX.EQ.0.0D0) GO TO 400
F = DABS(XX)**ETA
IF (XX.LT.0.0D0) F = F*(-1.0D0)**(IFIX(SNGL(ETA)))
400 CONTINUE
IF (IROUND.GT.0) F = ROUND(F,IROUND)
RETURN
END
DOUBLE PRECISION FUNCTION PHI(M,X)
DOUBLE PRECISION X(1), DI, XX
INTEGER I, M, N
N = M - 1
XX = X(1)
IF (N.GE.0) GO TO 100
PHI = 0.0D0
GO TO 400
100 CONTINUE
IF (N.GT.0) GO TO 200
PHI = 1.0D0
GO TO 400
200 CONTINUE
PHI = 1.0D0
DO 300 I = 1,N
DI = I
PHI = PHI*XX/DI
300 CONTINUE
400 CONTINUE
RETURN
END
SUBROUTINE RAND(X)
INTEGER MOD, INT
DOUBLE PRECISION X
INTEGER J, K, M, IX
INTEGER IRAND
REAL RM, XR, SEED
DATA SEED /0.0/
DATA K,J,M,RM /5701,3612,566927,566927.0/
IX = INT(SEED*RM)
IRAND = MOD(J*IX+K,M)
XR = IRAND
SEED = (XR+0.5)/RM
X = SEED
RETURN
END
DOUBLE PRECISION FUNCTION ROUND(X,N)
C SIMULATE ROUNDING TO N DIGITS
DOUBLE PRECISION X, EPS1, EPS2, FACTOR
INTEGER N
CALL RAND(EPS1)
CALL RAND(EPS2)
FACTOR = (1.0D0+(EPS1+EPS1-1.0D0)*10.0D0**(-N))
ROUND = X*FACTOR+(EPS2+EPS2-1.0D0)*10.0D0**(-N)
RETURN
END
DOUBLE PRECISION FUNCTION SPLN(X)
DOUBLE PRECISION DEXP
DOUBLE PRECISION E, X, A0, A1
DOUBLE PRECISION A2, E2, A3
E = DEXP(1.0D0)
E2 = E*E
IF (X.GT.1.0D0) GO TO 100
C X IS .LE. 1
A0 = 1.0D0
A1 = (-2.0D0*E2+12.0D0*E-9.0D0)/7.0D0
A2 = 0.0D0
A3 = (2.0D0*E**2-5.0D0*E+2.0D0)/7.0D0
GO TO 200
100 CONTINUE
C X IS .GT. 1
A0 = (5.0D0*E2-16.0D0*E+12.0D0)/7.0D0
A1 = (-17.0D0*E2+60.0D0*E-24.0D0)/7.0D0
A2 = (15.0D0*E2-48.0D0*E+15.0D0)/7.0D0
A3 = (-3.0D0*E2+11.0D0*E-3.0D0)/7.0D0
200 CONTINUE
SPLN = ((A3*X+A2)*X+A1)*X+A0
RETURN
END
.