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