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