home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 June / SIMTEL_0692.cdr / msdos / fortran / linpkdrv.arc / SEX.FOR < prev    next >
Encoding:
Text File  |  1984-01-14  |  10.1 KB  |  395 lines
  1. C     MAIN PROGRAM
  2.       INTEGER LUNIT
  3. C     ALLOW 5000 UNDERFLOWS.
  4. C     CALL TRAPS(0,0,5001,0,0)
  5. C
  6. C     OUTPUT UNIT NUMBER
  7. C
  8.       LUNIT = 6
  9. C
  10.       CALL SQREE(LUNIT)
  11.       STOP
  12.       END
  13.       SUBROUTINE SQREE(LUNIT)
  14.       REAL R(10,10),RX(11,10),RHR(10,10),Z(10,4),S(10)
  15.       REAL C(10)
  16.       LOGICAL NOTWRT
  17.       REAL SMACH,SCALE
  18.       COMMON SCALE,NOTWRT
  19.       INTEGER P
  20.       NOTWRT = .TRUE.
  21.       SCALE = 1.0E0
  22.       LDR = 10
  23.       LDX = 11
  24.       LDZ = 10
  25.       P = 10
  26.       WRITE (LUNIT,280) P
  27.       CALL SGETEX(RX,LDX,Z,LDZ,P,S,LUNIT)
  28.       CALL SRHXR(RX,LDX,P,RHR,LDR)
  29.       CALL SRHXZ(RX,LDX,P,Z(1,1),Z(1,3))
  30.       DO 80 I = 1, 3
  31.          GO TO (10,20,30), I
  32.    10    CONTINUE
  33.             K = 1
  34.             L = 10
  35.          GO TO 40
  36.    20    CONTINUE
  37.             K = 5
  38.             L = 6
  39.          GO TO 40
  40.    30    CONTINUE
  41.             K = 3
  42.             L = 7
  43.    40    CONTINUE
  44.          DO 70 JOB = 1, 2
  45.             DO 60 J1 = 1, P
  46.                Z(J1,2) = Z(J1,1)
  47.                DO 50 I1 = 1, P
  48.                   R(I1,J1) = RX(I1,J1)
  49.    50          CONTINUE
  50.    60       CONTINUE
  51.             CALL SCHEX(R,LDR,P,K,L,Z(1,2),LDZ,1,C,S,JOB)
  52.             CALL SCMPEX(R,LDR,RHR,LDR,Z(1,2),LDZ,P,JOB,K,L,LUNIT)
  53.    70    CONTINUE
  54.    80 CONTINUE
  55.       P = 2
  56.       WRITE (LUNIT,280) P
  57.       CALL SGETEX(RX,LDX,Z,LDZ,P,S,LUNIT)
  58.       CALL SRHXR(RX,LDX,P,RHR,LDR)
  59.       CALL SRHXZ(RX,LDX,P,Z(1,1),Z(1,3))
  60.       K = 1
  61.       L = 2
  62.       DO 110 JOB = 1, 2
  63.          DO 100 J1 = 1, P
  64.             Z(J1,2) = Z(J1,1)
  65.             DO 90 I1 = 1, P
  66.                R(I1,J1) = RX(I1,J1)
  67.    90       CONTINUE
  68.   100    CONTINUE
  69.          CALL SCHEX(R,LDR,P,K,L,Z(1,2),LDZ,1,C,S,JOB)
  70.          CALL SCMPEX(R,LDR,RHR,LDR,Z(1,2),LDZ,P,JOB,K,L,LUNIT)
  71.   110 CONTINUE
  72.       SCALE = SMACH(3)
  73.       P = 10
  74.       WRITE (LUNIT,280) P
  75.       WRITE (LUNIT,290)
  76.       CALL SGETEX(RX,LDX,Z,LDZ,P,S,LUNIT)
  77.       CALL SRHXR(RX,LDX,P,RHR,LDR)
  78.       CALL SRHXZ(RX,LDX,P,Z(1,1),Z(1,3))
  79.       DO 190 I = 1, 3
  80.          GO TO (120,130,140), I
  81.   120    CONTINUE
  82.             K = 1
  83.             L = 10
  84.          GO TO 150
  85.   130    CONTINUE
  86.             K = 5
  87.             L = 6
  88.          GO TO 150
  89.   140    CONTINUE
  90.             K = 3
  91.             L = 7
  92.   150    CONTINUE
  93.          DO 180 JOB = 1, 2
  94.             DO 170 J1 = 1, P
  95.                Z(J1,2) = Z(J1,1)
  96.                DO 160 I1 = 1, P
  97.                   R(I1,J1) = RX(I1,J1)
  98.   160          CONTINUE
  99.   170       CONTINUE
  100.             CALL SCHEX(R,LDR,P,K,L,Z(1,2),LDZ,1,C,S,JOB)
  101.             CALL SCMPEX(R,LDR,RHR,LDR,Z(1,2),LDZ,P,JOB,K,L,LUNIT)
  102.   180    CONTINUE
  103.   190 CONTINUE
  104.       SCALE = SMACH(2)
  105.       P = 10
  106.       WRITE (LUNIT,280) P
  107.       WRITE (LUNIT,300)
  108.       CALL SGETEX(RX,LDX,Z,LDZ,P,S,LUNIT)
  109.       CALL SRHXR(RX,LDX,P,RHR,LDR)
  110.       CALL SRHXZ(RX,LDX,P,Z(1,1),Z(1,3))
  111.       DO 270 I = 1, 3
  112.          GO TO (200,210,220), I
  113.   200    CONTINUE
  114.             K = 1
  115.             L = 10
  116.          GO TO 230
  117.   210    CONTINUE
  118.             K = 5
  119.             L = 6
  120.          GO TO 230
  121.   220    CONTINUE
  122.             K = 3
  123.             L = 7
  124.   230    CONTINUE
  125.          DO 260 JOB = 1, 2
  126.             DO 250 J1 = 1, P
  127.                Z(J1,2) = Z(J1,1)
  128.                DO 240 I1 = 1, P
  129.                   R(I1,J1) = RX(I1,J1)
  130.   240          CONTINUE
  131.   250       CONTINUE
  132.             CALL SCHEX(R,LDR,P,K,L,Z(1,2),LDZ,1,C,S,JOB)
  133.             CALL SCMPEX(R,LDR,RHR,LDR,Z(1,2),LDZ,P,JOB,K,L,LUNIT)
  134.   260    CONTINUE
  135.   270 CONTINUE
  136.       RETURN
  137. C
  138. C
  139.   280 FORMAT (1H1, 3X, 22H *** SCHEX ***     P =, I3 ///)
  140.   290 FORMAT (19H      OVERFLOW TEST ///)
  141.   300 FORMAT (20H      UNDERFLOW TEST ///)
  142.       END
  143.       SUBROUTINE SARRAY(A,LDA,M,N,NNL,LUNIT)
  144.       INTEGER LDA,M,N,NNL,LUNIT
  145.       REAL A(LDA,1)
  146. C     X
  147. C          FORTRAN IABS,MIN0
  148. C
  149.       INTEGER I,J,K,KU,NL
  150.       NL = IABS(NNL)
  151.       IF (NNL .LT. 0) GO TO 30
  152.          DO 20 I = 1, M
  153.             WRITE (LUNIT,70)
  154.             DO 10 K = 1, N, NL
  155.                KU = MIN0(K+NL-1,N)
  156.                WRITE (LUNIT,70) (A(I,J), J = K, KU)
  157.    10       CONTINUE
  158.    20    CONTINUE
  159.       GO TO 60
  160.    30 CONTINUE
  161.          DO 50 J = 1, N
  162.             WRITE (LUNIT,70)
  163.             DO 40 K = 1, M, NL
  164.                KU = MIN0(K+NL-1,M)
  165.                WRITE (LUNIT,70) (A(I,J), I = K, KU)
  166.    40       CONTINUE
  167.    50    CONTINUE
  168.    60 CONTINUE
  169.       RETURN
  170. C
  171.    70 FORMAT (1H , 4(2E13.6, 4X))
  172.       END
  173.       SUBROUTINE SCMPEX(R,LDR,RHR,LDRHR,Z,LDZ,P,JOB,K,L,LUNIT)
  174.       INTEGER LDR,LDRHR,LDZ,P,JOB,K,L,LUNIT
  175.       REAL R(LDR,1),RHR(LDRHR,1),Z(LDZ,1)
  176.       INTEGER I,J,JJ,LM1,PM1
  177.       REAL ERRCR,ERRCZ,ERMAX,EZMAX,RMAX,ZMAX
  178.       REAL T
  179.       LOGICAL NOTWRT
  180.       REAL SCALE,SMACH
  181.       COMMON SCALE,NOTWRT
  182.       IF (NOTWRT) GO TO 10
  183.          WRITE (LUNIT,150)
  184.          CALL SARRAY(R,LDR,P,P,P,LUNIT)
  185.          WRITE (LUNIT,160)
  186.          CALL SARRAY(Z,LDZ,P,1,-P,LUNIT)
  187.    10 CONTINUE
  188.       PM1 = P - 1
  189.       LM1 = L - 1
  190.       DO 30 J = 1, PM1
  191.          JJ = J + 1
  192.          DO 20 I = JJ, P
  193.             R(I,J) = 0.0E0
  194.    20    CONTINUE
  195.    30 CONTINUE
  196.       CALL SRHXZ(R,LDR,P,Z(1,1),Z(1,3))
  197.       CALL SRHXR(R,LDR,P,R,LDR)
  198.       DO 50 J = 1, P
  199.          DO 40 I = 1, P
  200.             R(I,J) = R(I,J)/SCALE
  201.    40    CONTINUE
  202.    50 CONTINUE
  203.       IF (JOB .EQ. 2) GO TO 80
  204.          DO 60 J = K, LM1
  205.             CALL SSWAP(P,R(1,J),1,R(1,J+1),1)
  206.             CALL SSWAP(1,Z(J,3),1,Z(J+1,3),1)
  207.    60    CONTINUE
  208.          DO 70 J = K, LM1
  209.             CALL SSWAP(P,R(J,1),LDR,R(J+1,1),LDR)
  210.    70    CONTINUE
  211.       GO TO 110
  212.    80 CONTINUE
  213.          DO 90 I = K, LM1
  214.             J = LM1 + K - I
  215.             CALL SSWAP(P,R(1,J),1,R(1,J+1),1)
  216.             CALL SSWAP(1,Z(J,3),1,Z(J+1,3),1)
  217.    90    CONTINUE
  218.          DO 100 I = K, LM1
  219.             J = LM1 + K - I
  220.             CALL SSWAP(P,R(J,1),LDR,R(J+1,1),LDR)
  221.   100    CONTINUE
  222.   110 CONTINUE
  223.       ERMAX = 0.0E0
  224.       EZMAX = 0.0E0
  225.       RMAX = 0.0E0
  226.       ZMAX = 0.0E0
  227.       DO 130 J = 1, P
  228.          T = Z(J,2)
  229.          ZMAX = AMAX1(ZMAX,AMAX1(ABS(T),0.0E0))
  230.          T = Z(J,2) - Z(J,3)
  231.          EZMAX = AMAX1(EZMAX,AMAX1(ABS(T),0.0E0))
  232.          DO 120 I = 1, P
  233.             T = RHR(I,J)
  234.             RMAX = AMAX1(RMAX,AMAX1(ABS(T),0.0E0))
  235.             T = RHR(I,J) - R(I,J)
  236.             ERMAX = AMAX1(ERMAX,AMAX1(ABS(T),0.0E0))
  237.   120    CONTINUE
  238.   130 CONTINUE
  239.       ERRCR = ERMAX/RMAX/SMACH(1)
  240.       ERRCZ = EZMAX/ZMAX/SMACH(1)
  241.       WRITE (LUNIT,140) JOB,K,L,ERRCR,ERRCZ
  242.       RETURN
  243.   140 FORMAT ( /// 25H     STATISTICS     JOB =, I2, 5H   K=, I2,
  244.      *         5H   L=, I2 / 38H0         RH*R    ................... ,
  245.      *         E10.2 / 38H          RH*Z    ................... ,
  246.      *         E10.2)
  247.   150 FORMAT (6H   REX)
  248.   160 FORMAT (6H   ZEX)
  249.       END
  250.       SUBROUTINE SGETEX(X,LDX,Z,LDZ,P,S,LUNIT)
  251.       INTEGER LDX,LDZ,P,LUNIT,JOB
  252.       LOGICAL NOTWRT
  253.       REAL SCALE
  254.       COMMON SCALE,NOTWRT
  255.       REAL X(LDX,1),Z(LDZ,1),S(1)
  256.       INTEGER I,JJ,J,PP1,PM1
  257.       PP1 = P + 1
  258.       PM1 = P - 1
  259.       CALL SGETRX(X,LDX,PP1,P,S)
  260.       DO 10 I = 1, P
  261.          Z(I,1) = X(PP1,I)*SCALE
  262.    10 CONTINUE
  263.       DO 30 J = 1, PM1
  264.          JJ = J + 1
  265.          DO 20 I = JJ, P
  266.             X(I,J) = 0.0E0
  267.    20    CONTINUE
  268.    30 CONTINUE
  269.       DO 50 J = 1, P
  270.          DO 40 I = 1, P
  271.             X(I,J) = X(I,J)*SCALE
  272.    40    CONTINUE
  273.    50 CONTINUE
  274.       IF (NOTWRT) GO TO 60
  275.          WRITE (LUNIT,80)
  276.          CALL SARRAY(X,LDX,P,P,P,LUNIT)
  277.          WRITE (LUNIT,70)
  278.          CALL SARRAY(Z,LDZ,P,1,-P,LUNIT)
  279.    60 CONTINUE
  280.       RETURN
  281. C
  282.    70 FORMAT ( /// 4H   Z)
  283.    80 FORMAT (4H   R)
  284.       END
  285.       SUBROUTINE SGETRX(X,LDX,N,P,S)
  286.       INTEGER N,P,LDX
  287.       REAL X(LDX,1),S(1)
  288.       INTEGER PD2,PD2D1,I
  289.       PD2 = P/2
  290.       PD2D1 = PD2 + 1
  291.       DO 10 I = 1, PD2
  292.          S(I) = 1.0E0
  293.    10 CONTINUE
  294.       IF (P .LT. PD2D1) GO TO 30
  295.       DO 20 I = PD2D1, P
  296.          S(I) = 0.5E0
  297.    20 CONTINUE
  298.    30 CONTINUE
  299.       CALL SXGEN(X,LDX,N,P,S)
  300.       RETURN
  301.       END
  302.       SUBROUTINE SRHXR(R,LDR,P,RHR,LDRHR)
  303.       INTEGER LDR,LDRHR,P
  304.       REAL R(LDR,1),RHR(LDRHR,1)
  305.       LOGICAL DUMMY
  306.       REAL SCALE
  307.       COMMON SCALE,DUMMY
  308.       REAL SDOT
  309.       DO 20 J = 1, P
  310.          DO 10 I = 1, J
  311.             R(I,J) = R(I,J)/SCALE
  312.    10    CONTINUE
  313.    20 CONTINUE
  314.       DO 40 J = 1, P
  315.          DO 30 I = J, P
  316.             II = P + J - I
  317.             RHR(II,J) = SDOT(J,R(1,II),1,R(1,J),1)
  318.    30    CONTINUE
  319.    40 CONTINUE
  320.       DO 60 J = 2, P
  321.          JJ = J - 1
  322.          DO 50 I = 1, JJ
  323.             RHR(I,J) = RHR(J,I)
  324.    50    CONTINUE
  325.    60 CONTINUE
  326.       DO 80 J = 1, P
  327.          DO 70 I = 1, P
  328.             R(I,J) = R(I,J)*SCALE
  329.    70    CONTINUE
  330.    80 CONTINUE
  331.       RETURN
  332.       END
  333.       SUBROUTINE SRHXZ(R,LDR,P,Z,RHZ)
  334.       INTEGER LDR,P
  335.       REAL R(LDR,1),Z(1),RHZ(1)
  336.       LOGICAL DUMMY
  337.       REAL SCALE
  338.       COMMON SCALE,DUMMY
  339.       REAL SDOT
  340.       DO 20 J = 1, P
  341.          Z(J) = Z(J)/SCALE
  342.          DO 10 I = 1, J
  343.             R(I,J) = R(I,J)/SCALE
  344.    10    CONTINUE
  345.    20 CONTINUE
  346.       DO 30 J = 1, P
  347.          RHZ(J) = SDOT(J,R(1,J),1,Z(1),1)
  348.    30 CONTINUE
  349.       DO 50 J = 1, P
  350.          Z(J) = Z(J)*SCALE
  351.          DO 40 I = 1, J
  352.             R(I,J) = R(I,J)*SCALE
  353.    40    CONTINUE
  354.    50 CONTINUE
  355.       RETURN
  356.       END
  357.       SUBROUTINE SXGEN(X,LDX,N,P,S)
  358.       INTEGER LDX,N,P
  359.       REAL X(LDX,1),S(1)
  360.       INTEGER I,J,M,MP1
  361.       REAL FN,FP
  362.       REAL FAC,RU,T
  363.       FP = FLOAT(P)
  364.       FN = FLOAT(N)
  365.       M = MIN0(N,P)
  366.       RU = 1.0E0
  367.       RU = COS(6.28E0/FLOAT(M+1))
  368.       FAC = 1.0E0/FP
  369.       DO 20 I = 1, M
  370.          FAC = FAC*RU
  371.          DO 10 J = 1, P
  372.             X(I,J) = -2.0E0*S(I)*FAC
  373.    10    CONTINUE
  374.          X(I,I) = X(I,I) + FP*S(I)*FAC
  375.    20 CONTINUE
  376.       IF (M .GE. N) GO TO 50
  377.          MP1 = M + 1
  378.          DO 40 J = 1, P
  379.             DO 30 I = MP1, N
  380.                X(I,J) = 0.0E0
  381.    30       CONTINUE
  382.    40    CONTINUE
  383.    50 CONTINUE
  384.       DO 80 J = 1, P
  385.          T = 0.0E0
  386.          DO 60 I = 1, N
  387.             T = T + X(I,J)
  388.    60    CONTINUE
  389.          DO 70 I = 1, N
  390.             X(I,J) = X(I,J) - 2.0E0*T/FN
  391.    70    CONTINUE
  392.    80 CONTINUE
  393.       RETURN
  394.       END
  395.