home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / linpack / ssv.for < prev    next >
Encoding:
Text File  |  1984-01-06  |  11.9 KB  |  458 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 SSVTS(LUNIT)
  11.       STOP
  12.       END
  13.       SUBROUTINE SSVTS(LUNIT)
  14. C     LUNIT IS THE OUTPUT UNIT NUMBER.
  15. C
  16. C     TESTS
  17. C        SSVDC
  18. C
  19. C     LINPACK. THIS VERSION DATED 08/14/78 .
  20. C     G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
  21. C
  22. C     SUBROUTINES AND FUNCTIONS
  23. C
  24. C     EXTERNAL SMACH,SSVT1,SXGEN
  25. C     FORTRAN FLOAT
  26. C
  27. C     INTERNAL VARIABLES
  28. C
  29.       INTEGER LUNIT
  30.       INTEGER I,J,N,P,LDX,LDU,LDV,CASE,NCASE
  31.       REAL X(25,25),XX(25,25),U(25,25),V(25,25),S(25),E(25),WORK(25)
  32.       REAL SMACH,HUGE,TINY
  33.       LOGICAL NOTWRT
  34.       LDU = 25
  35.       LDV = 25
  36.       LDX = 25
  37.       HUGE = SMACH(3)
  38.       TINY = SMACH(2)
  39.       NOTWRT = .TRUE.
  40.       NCASE = 12
  41.       WRITE (LUNIT,430)
  42.       DO 290 CASE = 1, NCASE
  43.          WRITE (LUNIT,300) CASE
  44.          GO TO (10, 40, 70, 90, 110, 130, 170, 210, 240, 250, 260,
  45.      *          270), CASE
  46. C
  47. C        BIDIAGONAL MATRIX WITH ZERO AT END.
  48. C
  49.    10    CONTINUE
  50.             WRITE (LUNIT,310)
  51.             N = 4
  52.             P = 4
  53.             DO 30 I = 1, 4
  54.                DO 20 J = 1, 4
  55.                   X(I,J) = 0.0E0
  56.    20          CONTINUE
  57.    30       CONTINUE
  58.             X(1,1) = 1.0E0
  59.             X(1,2) = 1.0E0
  60.             X(2,2) = 2.0E0
  61.             X(2,3) = 1.0E0
  62.             X(3,3) = 3.0E0
  63.             X(3,4) = 1.0E0
  64.             CALL SSVT1(X,LDX,N,P,S,E,U,LDU,V,LDV,WORK,XX,CASE,NOTWRT,
  65.      *                 LUNIT,11)
  66.          GO TO 280
  67. C
  68. C        BIDIAGONAL MATRIX WITH ZERO IN THE MIDDLE.
  69. C
  70.    40    CONTINUE
  71.             WRITE (LUNIT,320)
  72.             N = 5
  73.             P = 5
  74.             DO 60 I = 1, 5
  75.                DO 50 J = 1, 5
  76.                   X(I,J) = 0.0E0
  77.    50          CONTINUE
  78.    60       CONTINUE
  79.             X(1,1) = 1.0E0
  80.             X(1,2) = 1.0E0
  81.             X(2,3) = 1.0E0
  82.             X(3,3) = 2.0E0
  83.             X(3,4) = 1.0E0
  84.             X(4,4) = 3.0E0
  85.             X(4,5) = 1.0E0
  86.             X(5,5) = 4.0E0
  87.             CALL SSVT1(X,LDX,N,P,S,E,U,LDU,X,LDX,WORK,XX,CASE,NOTWRT,
  88.      *                 LUNIT,11)
  89.          GO TO 280
  90. C
  91. C        TEST CASE WITH N .GT. P.
  92. C
  93.    70    CONTINUE
  94.             WRITE (LUNIT,330)
  95.             N = 8
  96.             P = 4
  97.             DO 80 I = 1, 4
  98.                S(I) = FLOAT(I)
  99.    80       CONTINUE
  100.             CALL SXGEN(X,LDX,N,P,S)
  101.             CALL SSVT1(X,LDX,N,P,S,E,X,LDX,V,LDV,WORK,XX,CASE,NOTWRT,
  102.      *                 LUNIT,21)
  103.          GO TO 280
  104. C
  105. C        TEST CASE WITH N .LT. P.
  106. C
  107.    90    CONTINUE
  108.             WRITE (LUNIT,340)
  109.             N = 4
  110.             P = 8
  111.             DO 100 I = 1, 8
  112.                S(I) = FLOAT(I)
  113.   100       CONTINUE
  114.             CALL SXGEN(X,LDX,N,P,S)
  115.             CALL SSVT1(X,LDX,N,P,S,E,X,LDX,V,LDV,WORK,XX,CASE,NOTWRT,
  116.      *                 LUNIT,11)
  117.          GO TO 280
  118. C
  119. C        TEST CASE WITH N = P = LDX = LDU = LDV.
  120. C
  121.   110    CONTINUE
  122.             WRITE (LUNIT,350)
  123.             N = 25
  124.             P = 25
  125.             DO 120 I = 1, 25
  126.                S(I) = FLOAT(I)
  127.   120       CONTINUE
  128.             CALL SXGEN(X,LDX,N,P,S)
  129.             CALL SSVT1(X,LDX,N,P,S,E,X,LDX,V,LDV,WORK,XX,CASE,NOTWRT,
  130.      *                 LUNIT,11)
  131.          GO TO 280
  132. C
  133. C        TEST FOR OVERFLOW CONTROL.
  134. C
  135.   130    CONTINUE
  136.             WRITE (LUNIT,360)
  137.             N = 4
  138.             P = 8
  139.             DO 140 I = 1, 8
  140.                S(I) = FLOAT(I)
  141.   140       CONTINUE
  142.             CALL SXGEN(X,LDX,N,P,S)
  143.             DO 160 I = 1, 4
  144.                DO 150 J = 1, 8
  145.                   X(I,J) = HUGE*X(I,J)
  146.   150          CONTINUE
  147.   160       CONTINUE
  148.             CALL SSVT1(X,LDX,N,P,S,E,U,LDU,V,LDV,WORK,XX,CASE,NOTWRT,
  149.      *                 LUNIT,11)
  150.          GO TO 280
  151. C
  152. C        TEST FOR UNDERFLOW CONTROL.
  153. C
  154.   170    CONTINUE
  155.             WRITE (LUNIT,370)
  156.             N = 8
  157.             P = 4
  158.             DO 180 I = 1, 8
  159.                S(I) = FLOAT(I)
  160.   180       CONTINUE
  161.             CALL SXGEN(X,LDX,N,P,S)
  162.             DO 200 I = 1, 8
  163.                DO 190 J = 1, 4
  164.                   X(I,J) = TINY*X(I,J)
  165.   190          CONTINUE
  166.   200       CONTINUE
  167.             CALL SSVT1(X,LDX,N,P,S,E,U,LDU,V,LDV,WORK,XX,CASE,NOTWRT,
  168.      *                 LUNIT,11)
  169.          GO TO 280
  170. C
  171. C        ZERO MATRIX.
  172. C
  173.   210    CONTINUE
  174.             WRITE (LUNIT,380)
  175.             N = 8
  176.             P = 4
  177.             DO 230 I = 1, N
  178.                DO 220 J = 1, P
  179.                   X(I,J) = 0.0E0
  180.   220          CONTINUE
  181.   230       CONTINUE
  182.             CALL SSVT1(X,LDX,N,P,S,E,U,LDU,V,LDV,WORK,XX,CASE,NOTWRT,
  183.      *                 LUNIT,11)
  184.          GO TO 280
  185. C
  186. C        1X1 MATRIX.
  187. C
  188.   240    CONTINUE
  189.             WRITE (LUNIT,390)
  190.             N = 1
  191.             P = 1
  192.             X(1,1) = 2.0E0
  193.             CALL SSVT1(X,LDX,N,P,S,E,U,LDU,V,LDV,WORK,XX,CASE,NOTWRT,
  194.      *                 LUNIT,11)
  195.          GO TO 280
  196. C
  197. C        2X2 MATRIX.
  198. C
  199.   250    CONTINUE
  200.             WRITE (LUNIT,400)
  201.             N = 2
  202.             P = 2
  203.             X(1,1) = 3.0E0
  204.             X(1,2) = 1.0E0
  205.             X(2,1) = 1.0E0
  206.             X(2,2) = 2.0E0
  207.             CALL SSVT1(X,LDX,N,P,S,E,U,LDU,V,LDV,WORK,XX,CASE,NOTWRT,
  208.      *                 LUNIT,11)
  209.          GO TO 280
  210. C
  211. C        COLUMN VECTOR.
  212. C
  213.   260    CONTINUE
  214.             WRITE (LUNIT,410)
  215.             N = 4
  216.             P = 1
  217.             X(1,1) = 1.0E0
  218.             X(2,1) = 0.0E0
  219.             X(3,1) = 0.0E0
  220.             X(4,1) = 2.0E0
  221.             CALL SSVT1(X,LDX,N,P,S,E,U,LDU,V,LDV,WORK,XX,CASE,NOTWRT,
  222.      *                 LUNIT,11)
  223.          GO TO 280
  224. C
  225. C        ROW VECTOR.
  226. C
  227.   270    CONTINUE
  228.             WRITE (LUNIT,420)
  229.             N = 1
  230.             P = 4
  231.             X(1,1) = 0.0E0
  232.             X(1,2) = 1.0E0
  233.             X(1,3) = 2.0E0
  234.             X(1,4) = 3.0E0
  235.             CALL SSVT1(X,LDX,N,P,S,E,U,LDU,V,LDV,WORK,XX,CASE,NOTWRT,
  236.      *                 LUNIT,11)
  237.   280    CONTINUE
  238.   290 CONTINUE
  239.       WRITE (LUNIT,440)
  240.       RETURN
  241.   300 FORMAT ( / 5H1CASE, I3)
  242.   310 FORMAT ( / 35H BIDIAGONAL MATRIX WITH ZERO AT END)
  243.   320 FORMAT ( / 38H BIDIAGONAL MATRIX WITH ZERO IN MIDDLE)
  244.   330 FORMAT ( / 13H 8 X 4 MATRIX)
  245.   340 FORMAT ( / 13H 4 X 8 MATRIX)
  246.   350 FORMAT ( / 15H 25 X 25 MATRIX)
  247.   360 FORMAT ( / 14H OVERFLOW TEST)
  248.   370 FORMAT ( / 15H UNDERFLOW TEST)
  249.   380 FORMAT ( / 12H ZERO MATRIX)
  250.   390 FORMAT ( / 13H 1 X 1 MATRIX)
  251.   400 FORMAT ( / 13H 2 X 2 MATRIX)
  252.   410 FORMAT ( / 14H COLUMN VECTOR)
  253.   420 FORMAT ( / 11H ROW VECTOR)
  254.   430 FORMAT (22H1LINPACK TESTER, SSV** /
  255.      *        29H THIS VERSION DATED 08/14/78.)
  256.   440 FORMAT ( / 27H1END OF SINGULAR VALUE TEST)
  257.       END
  258.       SUBROUTINE SARRAY(A,LDA,M,N,NNL,LUNIT)
  259.       INTEGER LDA,M,N,NNL,LUNIT
  260.       REAL A(LDA,1)
  261. C
  262. C     FORTRAN IABS,MIN0
  263. C
  264.       INTEGER I,J,K,KU,NL
  265.       NL = IABS(NNL)
  266.       IF (NNL .LT. 0) GO TO 30
  267.          DO 20 I = 1, M
  268.             WRITE (6,70)
  269.             DO 10 K = 1, N, NL
  270.                KU = MIN0(K+NL-1,N)
  271.                WRITE (6,70) (A(I,J), J = K, KU)
  272.    10       CONTINUE
  273.    20    CONTINUE
  274.       GO TO 60
  275.    30 CONTINUE
  276.          DO 50 J = 1, N
  277.             WRITE (6,70)
  278.             DO 40 K = 1, M, NL
  279.                KU = MIN0(K+NL-1,M)
  280.                WRITE (6,70) (A(I,J), I = K, KU)
  281.    40       CONTINUE
  282.    50    CONTINUE
  283.    60 CONTINUE
  284.       RETURN
  285.    70 FORMAT (1H , 8E13.6)
  286.       END
  287.       SUBROUTINE CSVBM(XX,LDX,N,P,S,U,LDU,V,LDV,X,XSTAT)
  288.       INTEGER LDX,N,P,LDU,LDV
  289.       REAL XX(LDX,1),S(1),U(LDU,1),V(LDV,1),X(LDX,1)
  290.       REAL XSTAT
  291. C
  292. C     EXTERNAL SMACH
  293. C     FORTRAN AMAX1,ABS,MIN0
  294. C
  295.       INTEGER I,J,K,M
  296.       REAL T(25)
  297.       REAL SMACH,EMAX,XMAX
  298. C
  299.       M = MIN0(N,P)
  300.       DO 20 J = 1, P
  301.          DO 10 I = 1, M
  302.             X(I,J) = S(I)*V(J,I)
  303.    10    CONTINUE
  304.    20 CONTINUE
  305.       IF (N .LE. P) GO TO 50
  306.          M = P + 1
  307.          DO 40 J = 1, P
  308.             DO 30 I = M, N
  309.                X(I,J) = 0.0E0
  310.    30       CONTINUE
  311.    40    CONTINUE
  312.    50 CONTINUE
  313.       M = MIN0(N,P)
  314.       DO 90 J = 1, P
  315.          DO 70 I = 1, N
  316.             T(I) = 0.0E0
  317.             DO 60 K = 1, M
  318.                T(I) = T(I) + U(I,K)*X(K,J)
  319.    60       CONTINUE
  320.    70    CONTINUE
  321.          DO 80 I = 1, N
  322.             X(I,J) = T(I)
  323.    80    CONTINUE
  324.    90 CONTINUE
  325.       EMAX = 0.0E0
  326.       XMAX = 0.0E0
  327.       DO 110 J = 1, P
  328.          DO 100 I = 1, N
  329.             EMAX = AMAX1(EMAX,ABS(X(I,J)-XX(I,J)))
  330.             XMAX = AMAX1(XMAX,ABS(XX(I,J)))
  331.   100    CONTINUE
  332.   110 CONTINUE
  333.       XSTAT = 0.0E0
  334.       IF (EMAX .EQ. 0.0E0) GO TO 140
  335.          IF (XMAX .EQ. 0.0E0) GO TO 120
  336.             XSTAT = (EMAX/XMAX)/SMACH(1)
  337.          GO TO 130
  338.   120    CONTINUE
  339.             XSTAT = 1.0E10
  340.   130    CONTINUE
  341.   140 CONTINUE
  342.       RETURN
  343.       END
  344.       SUBROUTINE SSVOT(X,LDX,N,COL,TEST)
  345.       INTEGER LDX,N,COL
  346.       REAL X(LDX,1)
  347. C
  348. C     FORTRAN AMAX1
  349. C
  350.       INTEGER I,J,K
  351.       REAL ELM
  352.       REAL TEST,EMAX,SMACH
  353.       EMAX = 0.0E0
  354.       DO 30 I = 1, COL
  355.          DO 20 J = 1, COL
  356.             ELM = 0.0E0
  357.             DO 10 K = 1, N
  358.                ELM = ELM + X(K,I)*X(K,J)
  359.    10       CONTINUE
  360.             IF (I .EQ. J) ELM = 1.0E0 - ELM
  361.             EMAX = AMAX1(EMAX,ABS(ELM))
  362.    20    CONTINUE
  363.    30 CONTINUE
  364.       TEST = EMAX/SMACH(1)
  365.       RETURN
  366.       END
  367.       SUBROUTINE SSVT1(X,LDX,N,P,S,E,U,LDU,V,LDV,WORK,XX,CASE,NOTWRT,
  368.      *                 LUNIT,JOB)
  369.       INTEGER LDX,N,P,LDU,LDV,CASE,LUNIT,JOB
  370.       REAL X(LDX,1),S(1),E(1),U(LDU,1),V(LDV,1),WORK(1),XX(LDX,1)
  371.       LOGICAL NOTWRT
  372. C
  373. C     EXTERNAL SARRAY,SSVBM,SSVDC,SSVOT
  374. C     FORTRAN AMAX1,MIN0
  375. C
  376.       INTEGER I,J,M,UCOL,INFO
  377.       REAL USTAT,VSTAT,XSTAT
  378.       REAL XNEW(25,25)
  379.       REAL ERRLVL
  380.       ERRLVL = 100.0E0
  381.       UCOL = N
  382.       IF (JOB/10 .GE. 2) UCOL = P
  383.       DO 20 J = 1, P
  384.          DO 10 I = 1, N
  385.             XX(I,J) = X(I,J)
  386.    10    CONTINUE
  387.    20 CONTINUE
  388.       IF (NOTWRT) GO TO 30
  389.          WRITE (LUNIT,50)
  390.          CALL SARRAY(X,LDX,N,P,5,LUNIT)
  391.    30 CONTINUE
  392.       CALL SSVDC(X,LDX,N,P,S,E,U,LDU,V,LDV,WORK,JOB,INFO)
  393.       M = MIN0(N,P)
  394.       IF (NOTWRT) GO TO 40
  395.          WRITE (LUNIT,60)
  396.          CALL SARRAY(S,P,M,1,-5,LUNIT)
  397.          WRITE (LUNIT,70)
  398.          CALL SARRAY(U,LDU,N,N,5,LUNIT)
  399.          WRITE (LUNIT,80)
  400.          CALL SARRAY(V,LDV,P,P,5,LUNIT)
  401.    40 CONTINUE
  402.       CALL SSVOT(U,LDU,N,UCOL,USTAT)
  403.       CALL SSVOT(V,LDV,P,P,VSTAT)
  404.       CALL CSVBM(XX,LDX,N,P,S,U,LDU,V,LDV,XNEW,XSTAT)
  405.       WRITE (LUNIT,90) XSTAT,USTAT,VSTAT
  406.       IF (AMAX1(XSTAT,USTAT,VSTAT) .GT. ERRLVL) WRITE (LUNIT,100)
  407.       RETURN
  408.    50 FORMAT ( / 2H X)
  409.    60 FORMAT ( / 2H S)
  410.    70 FORMAT ( / 2H U)
  411.    80 FORMAT ( / 2H V)
  412.    90 FORMAT ( / 11H STATISTICS //
  413.      *         38H         U*SIGMA*VH .................., E10.2 /
  414.      *         38H         UHU ........................., E10.2 /
  415.      *         38H         VHV ........................., E10.2)
  416.   100 FORMAT ( / 35H ***** STATISTICS ABOVE ERROR LEVEL)
  417.       END
  418.       SUBROUTINE SXGEN(X,LDX,N,P,S)
  419.       INTEGER P,N,LDX
  420.       REAL X(LDX,1),S(1)
  421. C
  422. C     FORTRAN FLOAT,MIN0
  423. C
  424.       INTEGER I,J,M,MP1
  425.       REAL T,RU,FAC
  426.       REAL FP,FN
  427.       FP = FLOAT(P)
  428.       FN = FLOAT(N)
  429.       M = MIN0(N,P)
  430.       RU = COS(6.28E0/FLOAT(M+1))
  431.       FAC = 1.0E0/FP
  432.       DO 20 I = 1, M
  433.          FAC = FAC*RU
  434.          DO 10 J = 1, P
  435.             X(I,J) = -2.0E0*S(I)*FAC
  436.    10    CONTINUE
  437.          X(I,I) = X(I,I) + FP*S(I)*FAC
  438.    20 CONTINUE
  439.       IF (M .GE. N) GO TO 50
  440.          MP1 = M + 1
  441.          DO 40 J = 1, P
  442.             DO 30 I = MP1, N
  443.                X(I,J) = 0.0E0
  444.    30       CONTINUE
  445.    40    CONTINUE
  446.    50 CONTINUE
  447.       DO 80 J = 1, P
  448.          T = 0.0E0
  449.          DO 60 I = 1, N
  450.             T = T + X(I,J)
  451.    60    CONTINUE
  452.          DO 70 I = 1, N
  453.             X(I,J) = X(I,J) - 2.0E0*T/FN
  454.    70    CONTINUE
  455.    80 CONTINUE
  456.       RETURN
  457.       END
  458.