home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / CLIPPER / MISC / LINPKDRV.ZIP / SUD.FOR < prev    next >
Encoding:
Text File  |  1984-04-30  |  12.8 KB  |  440 lines
  1. C     MAIN PROGRAM
  2.       INTEGER LUNIT
  3. C     ALLOW 5000 UNDERFLOWS.
  4. C     CALL TRAPS(0,0,5001,0,0)  This is doesn't work on PC's!
  5. C
  6. C     OUTPUT UNIT NUMBER
  7. C
  8.       LUNIT = 6
  9. C
  10.       CALL SQRXX(LUNIT)
  11.       STOP
  12.       END
  13.       SUBROUTINE SQRXX(LUNIT)
  14.       INTEGER LUNIT
  15.       REAL RX(20,10),R(10,10),XROW(10),Z(10,4),YROW(10),S(10)
  16.       REAL SCALE,TINY,HUGE,RHO(2),C(10),SMACH
  17.       INTEGER N,P,LDX,LDR,LDZ,CASE
  18.       LOGICAL NOTWRT,OFLOW,UFLOW
  19.       COMMON SCALE,NOTWRT,OFLOW,UFLOW
  20.       LDZ = 10
  21.       LDX = 20
  22.       LDR = 10
  23.       NOTWRT = .TRUE.
  24.       OFLOW = .FALSE.
  25.       UFLOW = .FALSE.
  26.       TINY = SMACH(2)
  27.       HUGE = SMACH(3)
  28.       SCALE = 1.0E0
  29.       DO 60 CASE = 1, 3
  30.          GO TO (10, 20, 30), CASE
  31.    10    CONTINUE
  32.             N = 20
  33.             P = 10
  34.          GO TO 40
  35.    20    CONTINUE
  36.             N = 10
  37.             P = 4
  38.          GO TO 40
  39.    30    CONTINUE
  40.             N = 10
  41.             P = 1
  42.    40    CONTINUE
  43.          CALL SGETRX(RX,LDX,N,P,S)
  44.          WRITE (LUNIT,130) CASE,N,P
  45.          IF (NOTWRT) GO TO 50
  46.             WRITE (LUNIT,160)
  47.             CALL SARRAY(RX,LDX,N,P,P,LUNIT)
  48.    50    CONTINUE
  49.          CALL SDRUD1(R,LDR,RX,LDX,N,P,XROW,C,S,LUNIT)
  50.          CALL SDRUD2(R,LDR,RX,LDX,P,XROW,YROW,Z,LDZ,RHO,C,S,LUNIT)
  51.          CALL SDRDD(R,LDR,RX,LDX,P,Z,LDZ,XROW,YROW,RHO,C,S,LUNIT)
  52.    60 CONTINUE
  53.       CASE = 4
  54.       N = 10
  55.       P = 4
  56.       WRITE (LUNIT,130) CASE,N,P
  57.       WRITE (LUNIT,140)
  58.       CALL SGETRX(RX,LDX,N,P,S)
  59.       DO 80 J = 1, P
  60.          DO 70 I = 1, N
  61.             RX(I,J) = HUGE*RX(I,J)
  62.    70    CONTINUE
  63.    80 CONTINUE
  64.       SCALE = HUGE
  65.       OFLOW = .TRUE.
  66.       IF (NOTWRT) GO TO 90
  67.          WRITE (LUNIT,160)
  68.          CALL SARRAY(RX,LDX,N,P,P,LUNIT)
  69.    90 CONTINUE
  70.       CALL SDRUD1(R,LDR,RX,LDX,N,P,XROW,C,S,LUNIT)
  71.       CALL SDRUD2(R,LDR,RX,LDX,P,XROW,YROW,Z,LDZ,RHO,C,S,LUNIT)
  72.       CALL SDRDD(R,LDR,RX,LDX,P,Z,LDZ,XROW,YROW,RHO,C,S,LUNIT)
  73.       OFLOW = .FALSE.
  74.       CASE = 5
  75.       N = 10
  76.       P = 4
  77.       WRITE (LUNIT,130) CASE,N,P
  78.       WRITE (LUNIT,150)
  79.       CALL SGETRX(RX,LDX,N,P,S)
  80.       DO 110 J = 1, P
  81.          DO 100 I = 1, N
  82.             RX(I,J) = TINY*RX(I,J)
  83.   100    CONTINUE
  84.   110 CONTINUE
  85.       SCALE = TINY
  86.       UFLOW = .TRUE.
  87.       IF (NOTWRT) GO TO 120
  88.          WRITE (LUNIT,160)
  89.          CALL SARRAY(RX,LDX,N,P,P,LUNIT)
  90.   120 CONTINUE
  91.       CALL SDRUD1(R,LDR,RX,LDX,N,P,XROW,C,S,LUNIT)
  92.       CALL SDRUD2(R,LDR,RX,LDX,P,XROW,YROW,Z,LDZ,RHO,C,S,LUNIT)
  93.       CALL SDRDD(R,LDR,RX,LDX,P,Z,LDZ,XROW,YROW,RHO,C,S,LUNIT)
  94.       RETURN
  95. C
  96.   130 FORMAT (11H1    CASE =, I2, 5X, 3HN =, I2, 5X, 3HP =, I2 /////)
  97.   140 FORMAT (22H         OVERFLOW TEST /////)
  98.   150 FORMAT (24H          UNDERFLOW TEST /////)
  99.   160 FORMAT (5H   RX)
  100.       END
  101.       SUBROUTINE SARRAY(A,LDA,M,N,NNL,LUNIT)
  102.       INTEGER LDA,M,N,NNL,LUNIT
  103.       REAL A(LDA,1)
  104. C     X
  105. C          FORTRAN IABS,MIN0
  106. C
  107.       INTEGER I,J,K,KU,NL
  108.       NL = IABS(NNL)
  109.       IF (NNL .LT. 0) GO TO 30
  110.          DO 20 I = 1, M
  111.             WRITE (LUNIT,70)
  112.             DO 10 K = 1, N, NL
  113.                KU = MIN0(K+NL-1,N)
  114.                WRITE (LUNIT,70) (A(I,J), J = K, KU)
  115.    10       CONTINUE
  116.    20    CONTINUE
  117.       GO TO 60
  118.    30 CONTINUE
  119.          DO 50 J = 1, N
  120.             WRITE (LUNIT,70)
  121.             DO 40 K = 1, M, NL
  122.                KU = MIN0(K+NL-1,M)
  123.                WRITE (LUNIT,70) (A(I,J), I = K, KU)
  124.    40       CONTINUE
  125.    50    CONTINUE
  126.    60 CONTINUE
  127.       RETURN
  128. C
  129.    70 FORMAT (1H , 4(2E13.6, 4X))
  130.       END
  131.       SUBROUTINE SDRDD(R,LDR,RX,LDX,P,Z,LDZ,XROW,YROW,RHO,C,S,LUNIT)
  132.       INTEGER LDR,LDX,LDZ,P,LUNIT
  133.       REAL R(LDR,1),RX(LDX,1),Z(LDZ,1),XROW(1),YROW(1),S(1)
  134.       REAL RHO(1),SCALE,C(1)
  135.       LOGICAL NOTWRT,OFLOW,UFLOW
  136.       COMMON SCALE,NOTWRT,OFLOW,UFLOW
  137.       INTEGER I,J,JP1,PM1
  138.       WRITE (LUNIT,50)
  139.       CALL SCHDD(R,LDR,P,XROW,Z,LDZ,2,YROW,RHO,C,S,INFO)
  140.       PM1 = P - 1
  141.       IF (PM1 .LT. 1) GO TO 30
  142.       DO 20 J = 1, PM1
  143.          JP1 = J + 1
  144.          DO 10 I = JP1, P
  145.             R(I,J) = 0.0E0
  146.    10    CONTINUE
  147.    20 CONTINUE
  148.    30 CONTINUE
  149.       IF (NOTWRT) GO TO 40
  150.          WRITE (LUNIT,80)
  151.          CALL SARRAY(R,LDR,P,P,P,LUNIT)
  152.          WRITE (LUNIT,60)
  153.          CALL SARRAY(Z,LDZ,P,2,2,LUNIT)
  154.          WRITE (LUNIT,70) (RHO(I), I = 1, 2)
  155.    40 CONTINUE
  156.       CALL SMPDD(R,LDR,P,RX,LDX,Z,Z(1,3),LDZ,RHO,LUNIT)
  157.       RETURN
  158.    50 FORMAT ( /////
  159.      *         46H     STEP THREE    DOWNDATING XROW,YROW AND Z, ///)
  160.    60 FORMAT ( /// 4H   Z)
  161.    70 FORMAT ( /// 6H   RHO // 1H , 2E13.6)
  162.    80 FORMAT (4H   R)
  163.       END
  164.       SUBROUTINE SDRUD1(R,LDR,RX,LDX,N,P,XROW,C,S,LUNIT)
  165.       INTEGER N,P,LDR,LDX,LUNIT
  166.       REAL R(LDR,1),RX(LDX,1),XROW(1),S(1)
  167.       REAL C(1)
  168.       LOGICAL NOTWRT,OFLOW
  169.       REAL RELM,XELM,Y,Z
  170.       REAL XMAX,RMAX,TEST,T,SCALE,SMACH
  171.       COMMON SCALE,NOTWRT,OFLOW,UFLOW
  172.       DO 20 I = 1, P
  173.          DO 10 J = 1, P
  174.             R(I,J) = 0.0E0
  175.    10    CONTINUE
  176.    20 CONTINUE
  177.       DO 40 I = 1, N
  178.          DO 30 J = 1, P
  179.             XROW(J) = RX(I,J)
  180.    30    CONTINUE
  181.          CALL SCHUD(R,LDR,P,XROW,Z,1,0,Y,Y,C,S)
  182.    40 CONTINUE
  183.       WRITE (LUNIT,100)
  184.       IF (NOTWRT) GO TO 50
  185.          WRITE (LUNIT,120)
  186.          CALL SARRAY(R,LDR,P,P,P,LUNIT)
  187.    50 CONTINUE
  188.       RMAX = 0.0E0
  189.       XMAX = 0.0E0
  190.       DO 90 I = 1, P
  191.          DO 80 J = 1, I
  192.             RELM = 0.0E0
  193.             XELM = 0.0E0
  194.             NIM = MIN0(I,J)
  195.             DO 60 K = 1, NIM
  196.                RELM = RELM + R(K,I)/SCALE*(R(K,J)/SCALE)
  197.    60       CONTINUE
  198.             DO 70 K = 1, N
  199.                XELM = XELM + RX(K,I)/SCALE*(RX(K,J)/SCALE)
  200.    70       CONTINUE
  201.             T = AMAX1(ABS(XELM),0.0E0)
  202.             XMAX = AMAX1(XMAX,T)
  203.             T = AMAX1(ABS(XELM-RELM),0.0E0)
  204.             RMAX = AMAX1(RMAX,T)
  205.    80    CONTINUE
  206.    90 CONTINUE
  207.       TEST = RMAX/XMAX/SMACH(1)
  208.       WRITE (LUNIT,110) TEST
  209.       RETURN
  210. C
  211.   100 FORMAT ( ///// 25H    STEP ONE   UPDATING X ///)
  212.   110 FORMAT ( /// 15H     STATISTICS //
  213.      *         42H      RH*R    ............................, E10.2)
  214.   120 FORMAT (4H   R)
  215.       END
  216.       SUBROUTINE SDRUD2(R,LDR,RX,LDX,P,XROW,YROW,Z,LDZ,RHO,C,S,LUNIT)
  217.       INTEGER LDR,LDX,LDZ,P,LUNIT
  218.       REAL R(LDR,1),RX(LDX,1),XROW(1),YROW(1),Z(LDZ,1),S(1)
  219.       REAL RHO(1),C(1)
  220.       REAL SCALE
  221.       LOGICAL NOTWRT,OFLOW,UFLOW
  222.       COMMON SCALE,NOTWRT,OFLOW,UFLOW
  223.       DO 20 I = 1, P
  224.          DO 10 J = 1, P
  225.             RX(I,J) = R(I,J)
  226.    10    CONTINUE
  227.    20 CONTINUE
  228.       DO 40 I = 1, P
  229.          Z(I,1) = 0.0E0
  230.          DO 30 J = I, P
  231.             Z(I,1) = Z(I,1) + R(I,J)
  232.    30    CONTINUE
  233.          Z(I,2) = Z(I,1) + SCALE
  234.          Z(I,3) = Z(I,1)
  235.          Z(I,4) = Z(I,2)
  236.    40 CONTINUE
  237.       DO 60 I = 1, P
  238.          XROW(I) = 0.0E0
  239.          DO 50 J = 1, I
  240.             XROW(I) = XROW(I) + R(J,I)
  241.    50    CONTINUE
  242.    60 CONTINUE
  243.       YROW(1) = 0.0E0
  244.       DO 70 I = 1, P
  245.          YROW(1) = YROW(1) + XROW(I)
  246.    70 CONTINUE
  247.       YROW(2) = YROW(1) - SCALE
  248.       RHO(1) = SCALE
  249.       RHO(2) = SCALE
  250.       IF (NOTWRT) GO TO 80
  251.          WRITE (LUNIT,100)
  252.          CALL SARRAY(XROW,P,P,1,-P,LUNIT)
  253.          WRITE (LUNIT,110)
  254.          CALL SARRAY(Z,LDZ,P,2,2,LUNIT)
  255.          WRITE (LUNIT,120)
  256.          CALL SARRAY(YROW,2,2,1,-2,LUNIT)
  257.    80 CONTINUE
  258.       CALL SCHUD(R,LDR,P,XROW,Z,LDZ,2,YROW,RHO,C,S)
  259.       WRITE (LUNIT,130)
  260.       IF (NOTWRT) GO TO 90
  261.          WRITE (LUNIT,150)
  262.          CALL SARRAY(R,LDR,P,P,P,LUNIT)
  263.          WRITE (LUNIT,110)
  264.          CALL SARRAY(Z,LDZ,P,2,2,LUNIT)
  265.          WRITE (LUNIT,140) (RHO(I), I = 1, 2)
  266.    90 CONTINUE
  267.       CALL SMPUD(R,LDR,RX,LDX,P,XROW,Z,LDZ,RHO,LUNIT)
  268.       RETURN
  269. C
  270.   100 FORMAT ( ///// 7H   XROW)
  271.   110 FORMAT ( /// 4H   Z)
  272.   120 FORMAT ( /// 7H   YROW)
  273.   130 FORMAT ( ///// 41H     STEP TWO   UPDATING XROW, YROW AND Z ///)
  274.   140 FORMAT ( /// 6H   RHO // 1H , 2E13.6)
  275.   150 FORMAT (4H   R)
  276.       END
  277.       SUBROUTINE SGETRX(X,LDX,N,P,S)
  278.       INTEGER N,P,LDX
  279.       REAL X(LDX,1),S(1)
  280.       INTEGER PD2,PD2D1,I
  281.       PD2 = MAX0(P/2,1)
  282.       PD2D1 = PD2 + 1
  283.       DO 10 I = 1, PD2
  284.          S(I) = 1.0E0
  285.    10 CONTINUE
  286.       IF (P .LT. PD2D1) GO TO 30
  287.       DO 20 I = PD2D1, P
  288.          S(I) = 0.5E0
  289.    20 CONTINUE
  290.    30 CONTINUE
  291.       CALL SXGEN(X,LDX,N,P,S)
  292.       RETURN
  293.       END
  294.       SUBROUTINE SMPDD(R,LDR,P,ROLD,LDRD,Z,ZUP,LDZ,RHO,LUNIT)
  295.       INTEGER LDR,P,LDRD,LDZ,LUNIT
  296.       REAL R(LDR,1),ROLD(LDRD,1),Z(LDZ,1),ZUP(LDZ,1)
  297.       REAL RHO(1),SMACH
  298.       REAL T
  299.       INTEGER I,J
  300.       REAL TR,TZ(2),TRHO(2),MAXOLD,MAXREL,SCALE
  301.       LOGICAL NOTWRT,OFLOW,UFLOW
  302.       COMMON SCALE,NOTWRT,OFLOW,UFLOW
  303.       MAXOLD = 0.0E0
  304.       MAXREL = 0.0E0
  305.       DO 20 J = 1, P
  306.          DO 10 I = 1, J
  307.             T = ROLD(I,J)
  308.             MAXOLD = AMAX1(MAXOLD,AMAX1(ABS(T),0.0E0))
  309.             T = ROLD(I,J) - R(I,J)
  310.             MAXREL = AMAX1(MAXREL,AMAX1(ABS(T),0.0E0))
  311.    10    CONTINUE
  312.    20 CONTINUE
  313.       TR = MAXREL/MAXOLD/SMACH(1)
  314.       DO 40 I = 1, 2
  315.          MAXOLD = 0.0E0
  316.          MAXREL = 0.0E0
  317.          DO 30 J = 1, P
  318.             T = ZUP(J,I)
  319.             MAXOLD = AMAX1(MAXOLD,AMAX1(ABS(T),0.0E0))
  320.             T = ZUP(J,I) - Z(J,I)
  321.             MAXREL = AMAX1(MAXREL,AMAX1(ABS(T),0.0E0))
  322.    30    CONTINUE
  323.          TZ(I) = MAXREL/MAXOLD/SMACH(1)
  324.          TRHO(I) = ABS(RHO(I)/SCALE-1.0E0)/SMACH(1)
  325.    40 CONTINUE
  326.       WRITE (LUNIT,50) TR,TZ,TRHO
  327.       RETURN
  328. C
  329.    50 FORMAT ( /// 25H     STATSTICS STEP THREE //
  330.      *         33H        R   ....................., E10.2 /
  331.      *         33H        Z(1)   .................., E10.2 /
  332.      *         33H        Z(2)   .................., E10.2 /
  333.      *         33H        RHO(1)   ................, E10.2 /
  334.      *         33H        RHO(2)   ................, E10.2)
  335.       END
  336.       SUBROUTINE SMPUD(R,LDR,RX,LDX,P,XROW,Z,LDZ,RHO,LUNIT)
  337.       INTEGER LDR,LDX,P,LDZ
  338.       REAL R(LDR,1),RX(LDX,1),XROW(1),Z(LDZ,1)
  339.       REAL RHO(1)
  340.       REAL RELM,XELM,TMP1(10),TMP
  341.       REAL TEST1,TEST2(2),TEST3,TEST4,MAXR,MAXX,T
  342.       LOGICAL NOTWRT,OFLOW,UFLOW
  343.       REAL SCALE,SMACH
  344.       COMMON SCALE,NOTWRT,OFLOW,UFLOW
  345.       MAXR = 0.0E0
  346.       MAXX = 0.0E0
  347.       DO 30 I = 1, P
  348.          DO 20 J = 1, I
  349.             RELM = 0.0E0
  350.             XELM = 0.0E0
  351.             NIM = MIN0(I,J)
  352.             DO 10 K = 1, NIM
  353.                RELM = RELM + R(K,I)/SCALE*(R(K,J)/SCALE)
  354.                XELM = XELM + RX(K,I)/SCALE*(RX(K,J)/SCALE)
  355.    10       CONTINUE
  356.             XELM = XELM + XROW(I)/SCALE*(XROW(J)/SCALE)
  357.             T = AMAX1(ABS(XELM),0.0E0)
  358.             MAXX = AMAX1(MAXX,T)
  359.             T = AMAX1(ABS(XELM-RELM),0.0E0)
  360.             MAXR = AMAX1(MAXR,T)
  361.    20    CONTINUE
  362.    30 CONTINUE
  363.       TEST1 = MAXR/MAXX/SMACH(1)
  364.       TMP = 0.0E0
  365.       DO 50 I = 1, P
  366.          TMP1(I) = 0.0E0
  367.          DO 40 J = I, P
  368.             TMP1(I) = TMP1(I) + RX(I,J)
  369.    40    CONTINUE
  370.          TMP = TMP + XROW(I)
  371.    50 CONTINUE
  372.       DO 80 K = 1, 2
  373.          MAXR = 0.0E0
  374.          MAXX = 0.0E0
  375.          DO 70 I = 1, P
  376.             RELM = 0.0E0
  377.             XELM = 0.0E0
  378.             DO 60 J = 1, I
  379.                RELM = RELM + R(J,I)/SCALE*(Z(J,K)/SCALE)
  380.                XELM = XELM + RX(J,I)/SCALE*(TMP1(J)/SCALE)
  381.    60       CONTINUE
  382.             XELM = XELM + XROW(I)/SCALE*(TMP/SCALE)
  383.             T = AMAX1(ABS(XELM-RELM),0.0E0)
  384.             MAXR = AMAX1(MAXR,T)
  385.             T = AMAX1(ABS(XELM),0.0E0)
  386.             MAXX = AMAX1(MAXX,T)
  387.    70    CONTINUE
  388.          TEST2(K) = MAXR/MAXX/SMACH(1)
  389.    80 CONTINUE
  390.       TEST3 = ABS(RHO(1)/SCALE-1.0E0)/SMACH(1)
  391.       T = SQRT(2.0E0+FLOAT(P))
  392.       TEST4 = ABS(RHO(2)/SCALE-T)/T/SMACH(1)
  393.       WRITE (LUNIT,90) TEST1,TEST2,TEST3,TEST4
  394.       RETURN
  395.    90 FORMAT ( /// 24H     STATISTICS STEP TWO //
  396.      *         31H        RH*R   ................, E10.2 /
  397.      *         31H        RH*Z(1)   ............., E10.2 /
  398.      *         31H        RH*Z(2)   ............., E10.2 /
  399.      *         31H        RHO(1)   .............., E10.2 /
  400.      *         31H        RHO(2)   .............., E10.2)
  401.       END
  402.       SUBROUTINE SXGEN(X,LDX,N,P,S)
  403.       INTEGER LDX,N,P
  404.       REAL X(LDX,1),S(1)
  405.       INTEGER I,J,M,MP1
  406.       REAL FN,FP
  407.       REAL FAC,RU,T
  408.       FP = FLOAT(P)
  409.       FN = FLOAT(N)
  410.       M = MIN0(N,P)
  411.       RU = 1.0E0
  412.       RU = COS(6.28E0/FLOAT(M+1))
  413.       FAC = 1.0E0/FP
  414.       DO 20 I = 1, M
  415.          FAC = FAC*RU
  416.          DO 10 J = 1, P
  417.             X(I,J) = -2.0E0*S(I)*FAC
  418.    10    CONTINUE
  419.          X(I,I) = X(I,I) + FP*S(I)*FAC
  420.    20 CONTINUE
  421.       IF (M .GE. N) GO TO 50
  422.          MP1 = M + 1
  423.          DO 40 J = 1, P
  424.             DO 30 I = MP1, N
  425.                X(I,J) = 0.0E0
  426.    30       CONTINUE
  427.    40    CONTINUE
  428.    50 CONTINUE
  429.       DO 80 J = 1, P
  430.          T = 0.0E0
  431.          DO 60 I = 1, N
  432.             T = T + X(I,J)
  433.    60    CONTINUE
  434.          DO 70 I = 1, N
  435.             X(I,J) = X(I,J) - 2.0E0*T/FN
  436.    70    CONTINUE
  437.    80 CONTINUE
  438.       RETURN
  439.       END
  440.