home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / linpack / sch.for < prev    next >
Text File  |  1985-01-13  |  12KB  |  414 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 SCHTS(LUNIT)
  11.       STOP
  12.       END
  13.       SUBROUTINE SCHTS(LUNIT)
  14. C     LUNIT IS THE OUTPUT UNIT NUMBER
  15. C
  16. C     TESTS
  17. C        SCHDC
  18. C
  19. C     LINPACK. THIS VERSION DATED 08/14/78 .
  20. C     J.J. DONGARRA AND G.W. STEWART, ARGONNE NATIONAL LABORATORY,
  21. C     UNIVERSITY OF MARYLAND
  22. C
  23. C     SUBROUTINES AND FUNCTIONS
  24. C
  25. C     LINPACK SCHDC
  26. C     EXTERNAL SAGEN,SCHXX,SMACH,SARRAY
  27. C     FORTRAN AMAX1,ABS,FLOAT
  28. C
  29. C     INTERNAL VARIABLES
  30. C
  31.       INTEGER LUNIT
  32.       INTEGER CASE,LDA,P,JPVT(25),KPVT(25),KD
  33.       REAL RESDUL
  34.       REAL A(25,25),AA(25,25),ASAVE(25,25),D(25),WORK(25)
  35.       REAL SMACH,TINY,HUGE,EPS
  36.       LOGICAL NOTWRT
  37.       LDA = 25
  38.       NCASE = 7
  39.       ERRLVL = 100.0E0
  40.       NOTWRT = .TRUE.
  41.       EPS = SMACH(1)
  42.       TINY = SMACH(2)*10000.0E0
  43.       HUGE = SMACH(3)
  44.       WRITE (LUNIT,300)
  45.       DO 290 CASE = 1, NCASE
  46.          WRITE (LUNIT,310) CASE
  47.          GO TO (10,50,80,120,160,200,240), CASE
  48.    10    CONTINUE
  49. C
  50. C           5 X 5 CASE NO PIVOTING.
  51. C
  52.             WRITE (LUNIT,320)
  53.             P = 5
  54. C
  55. C           GENERATE MATRIX WITH POSITIVE EIGENVALUES.
  56. C
  57.             DO 20 I = 1, P
  58.                D(I) = FLOAT(I)
  59.                JPVT(I) = 0
  60.    20       CONTINUE
  61.             CALL SAGEN(A,LDA,P,D,ASAVE)
  62.             IF (NOTWRT) GO TO 30
  63.                WRITE (LUNIT,390)
  64.                CALL SARRAY(A,LDA,P,P,P,LUNIT)
  65.    30       CONTINUE
  66. C
  67. C           DECOMPOSE THE MATRIX.
  68. C
  69.             JOB = 0
  70.             WRITE (LUNIT,420) JOB,(JPVT(I), I = 1, P)
  71.             CALL SCHDC(A,LDA,P,WORK,JPVT,JOB,KD)
  72.             WRITE (LUNIT,430) KD
  73.             CALL SCHXX(A,LDA,P,AA,JPVT,KPVT,ASAVE,RESDUL,JOB)
  74.             WRITE (LUNIT,400) (JPVT(I), I = 1, P)
  75.             RESDUL = RESDUL/EPS
  76.             WRITE (LUNIT,410) RESDUL
  77.             IF (NOTWRT) GO TO 40
  78.                WRITE (LUNIT,390)
  79.                CALL SARRAY(A,LDA,P,P,P,LUNIT)
  80.    40       CONTINUE
  81.          GO TO 280
  82.    50    CONTINUE
  83. C
  84. C           1 X 1 CASE WITH PIVOTING.
  85. C
  86.             WRITE (LUNIT,330)
  87.             P = 1
  88. C
  89. C           GENERATE MATRIX WITH POSITIVE EIGENVALUES.
  90. C
  91.             A(1,1) = 1.0E0
  92.             ASAVE(1,1) = A(1,1)
  93.             JPVT(1) = 0
  94.             IF (NOTWRT) GO TO 60
  95.                WRITE (LUNIT,390)
  96.                CALL SARRAY(A,LDA,P,P,P,LUNIT)
  97.    60       CONTINUE
  98. C
  99. C           DECOMPOSE THE MATRIX.
  100. C
  101.             JOB = 1
  102.             WRITE (LUNIT,420) JOB,(JPVT(I), I = 1, P)
  103.             CALL SCHDC(A,LDA,P,WORK,JPVT,JOB,KD)
  104.             WRITE (LUNIT,430) KD
  105.             CALL SCHXX(A,LDA,P,AA,JPVT,KPVT,ASAVE,RESDUL,JOB)
  106.             WRITE (LUNIT,400) (JPVT(I), I = 1, P)
  107.             RESDUL = RESDUL/EPS
  108.             WRITE (LUNIT,410) RESDUL
  109.             IF (NOTWRT) GO TO 70
  110.                WRITE (LUNIT,390)
  111.                CALL SARRAY(A,LDA,P,P,P,LUNIT)
  112.    70       CONTINUE
  113.          GO TO 280
  114.    80    CONTINUE
  115. C
  116. C           8 X 8 CASE PIVOT LOGIC TEST
  117. C
  118.             WRITE (LUNIT,340)
  119.             P = 8
  120. C
  121. C           GENERATE MATRIX WITH POSITIVE EIGENVALUES.
  122. C
  123.             DO 90 I = 1, P
  124.                D(I) = FLOAT(I)
  125.                JPVT(I) = 1 - 2*MOD(I+1,2)
  126.    90       CONTINUE
  127.             CALL SAGEN(A,LDA,P,D,ASAVE)
  128.             IF (NOTWRT) GO TO 100
  129.                WRITE (LUNIT,390)
  130.                CALL SARRAY(A,LDA,P,P,P,LUNIT)
  131.   100       CONTINUE
  132. C
  133. C           DECOMPOSE THE MATRIX.
  134. C
  135.             JOB = 1
  136.             WRITE (LUNIT,420) JOB,(JPVT(I), I = 1, P)
  137.             CALL SCHDC(A,LDA,P,WORK,JPVT,JOB,KD)
  138.             WRITE (LUNIT,430) KD
  139.             CALL SCHXX(A,LDA,P,AA,JPVT,KPVT,ASAVE,RESDUL,JOB)
  140.             WRITE (LUNIT,400) (JPVT(I), I = 1, P)
  141.             RESDUL = RESDUL/EPS
  142.             WRITE (LUNIT,410) RESDUL
  143.             IF (NOTWRT) GO TO 110
  144.                WRITE (LUNIT,390)
  145.                CALL SARRAY(A,LDA,P,P,P,LUNIT)
  146.   110       CONTINUE
  147.          GO TO 280
  148.   120    CONTINUE
  149. C
  150. C           6 X 6 CASE PIVOTING WITH NEGATIVE EIGENVALUE.
  151. C
  152.             WRITE (LUNIT,350)
  153.             P = 6
  154. C
  155. C           GENERATE MATRIX
  156. C
  157.             D(1) = 1.0E0
  158.             D(2) = 1.0E0/2.0E0
  159.             D(3) = 1.0E0/4.0E0
  160.             D(4) = -1.0E0
  161.             D(5) = 1.0E0/2.0E0
  162.             D(6) = 1.0E0/4.0E0
  163.             DO 130 I = 1, P
  164.                JPVT(I) = 0
  165.   130       CONTINUE
  166.             CALL SAGEN(A,LDA,P,D,ASAVE)
  167.             IF (NOTWRT) GO TO 140
  168.                WRITE (LUNIT,390)
  169.                CALL SARRAY(A,LDA,P,P,P,LUNIT)
  170.   140       CONTINUE
  171. C
  172. C           DECOMPOSE THE MATRIX.
  173. C
  174.             JOB = 1
  175.             WRITE (LUNIT,420) JOB,(JPVT(I), I = 1, P)
  176.             CALL SCHDC(A,LDA,P,WORK,JPVT,JOB,KD)
  177.             WRITE (LUNIT,430) KD
  178.             CALL SCHXX(A,LDA,P,AA,JPVT,KPVT,ASAVE,RESDUL,JOB)
  179.             WRITE (LUNIT,400) (JPVT(I), I = 1, P)
  180.             RESDUL = RESDUL/EPS
  181.             WRITE (LUNIT,410) RESDUL
  182.             IF (NOTWRT) GO TO 150
  183.                WRITE (LUNIT,390)
  184.                CALL SARRAY(A,LDA,P,P,P,LUNIT)
  185.   150       CONTINUE
  186.          GO TO 280
  187.   160    CONTINUE
  188. C
  189. C           25 X 25 CASE WITH PIVOTING.
  190. C
  191.             WRITE (LUNIT,360)
  192.             P = 25
  193. C
  194. C           GENERATE MATRIX WITH POSITIVE EIGENVALUES.
  195. C
  196.             DO 170 I = 1, P
  197.                D(I) = FLOAT(I)
  198.                JPVT(I) = 0
  199.   170       CONTINUE
  200.             CALL SAGEN(A,LDA,P,D,ASAVE)
  201.             IF (NOTWRT) GO TO 180
  202.                WRITE (LUNIT,390)
  203.                CALL SARRAY(A,LDA,P,P,P,LUNIT)
  204.   180       CONTINUE
  205. C
  206. C           DECOMPOSE THE MATRIX.
  207. C
  208.             JOB = 1
  209.             WRITE (LUNIT,420) JOB,(JPVT(I), I = 1, P)
  210.             CALL SCHDC(A,LDA,P,WORK,JPVT,JOB,KD)
  211.             WRITE (LUNIT,430) KD
  212.             CALL SCHXX(A,LDA,P,AA,JPVT,KPVT,ASAVE,RESDUL,JOB)
  213.             WRITE (LUNIT,400) (JPVT(I), I = 1, P)
  214.             RESDUL = RESDUL/EPS
  215.             WRITE (LUNIT,410) RESDUL
  216.             IF (NOTWRT) GO TO 190
  217.                WRITE (LUNIT,390)
  218.                CALL SARRAY(A,LDA,P,P,P,LUNIT)
  219.   190       CONTINUE
  220.          GO TO 280
  221.   200    CONTINUE
  222. C
  223. C           5 X 5 CASE WITH PIVOTING AND UNDERFLOW CHECK.
  224. C
  225.             WRITE (LUNIT,370)
  226.             P = 5
  227. C
  228. C           GENERATE MATRIX WITH POSITIVE EIGENVALUES.
  229. C
  230.             DO 210 I = 1, P
  231.                D(I) = FLOAT(I)*TINY
  232.                JPVT(I) = 0
  233.   210       CONTINUE
  234.             CALL SAGEN(A,LDA,P,D,ASAVE)
  235.             IF (NOTWRT) GO TO 220
  236.                WRITE (LUNIT,390)
  237.                CALL SARRAY(A,LDA,P,P,P,LUNIT)
  238.   220       CONTINUE
  239. C
  240. C           DECOMPOSE THE MATRIX.
  241. C
  242.             JOB = 1
  243.             WRITE (LUNIT,420) JOB,(JPVT(I), I = 1, P)
  244.             CALL SCHDC(A,LDA,P,WORK,JPVT,JOB,KD)
  245.             WRITE (LUNIT,430) KD
  246.             CALL SCHXX(A,LDA,P,AA,JPVT,KPVT,ASAVE,RESDUL,JOB)
  247.             WRITE (LUNIT,400) (JPVT(I), I = 1, P)
  248.             RESDUL = RESDUL/EPS
  249.             WRITE (LUNIT,410) RESDUL
  250.             IF (NOTWRT) GO TO 230
  251.                WRITE (LUNIT,390)
  252.                CALL SARRAY(A,LDA,P,P,P,LUNIT)
  253.   230       CONTINUE
  254.          GO TO 280
  255.   240    CONTINUE
  256. C
  257. C           5 X 5 CASE WITH PIVOTING AND OVERFLOW CHECK.
  258. C
  259.             WRITE (LUNIT,380)
  260.             P = 5
  261. C
  262. C           GENERATE MATRIX WITH POSITIVE EIGENVALUES.
  263. C
  264.             DO 250 I = 1, P
  265.                D(I) = FLOAT(I)*HUGE
  266.                JPVT(I) = 0
  267.   250       CONTINUE
  268.             CALL SAGEN(A,LDA,P,D,ASAVE)
  269.             IF (NOTWRT) GO TO 260
  270.                WRITE (LUNIT,390)
  271.                CALL SARRAY(A,LDA,P,P,P,LUNIT)
  272.   260       CONTINUE
  273. C
  274. C           DECOMPOSE THE MATRIX.
  275. C
  276.             JOB = 1
  277.             WRITE (LUNIT,420) JOB,(JPVT(I), I = 1, P)
  278.             CALL SCHDC(A,LDA,P,WORK,JPVT,JOB,KD)
  279.             WRITE (LUNIT,430) KD
  280.             CALL SCHXX(A,LDA,P,AA,JPVT,KPVT,ASAVE,RESDUL,JOB)
  281.             WRITE (LUNIT,400) (JPVT(I), I = 1, P)
  282.             RESDUL = RESDUL/EPS
  283.             WRITE (LUNIT,410) RESDUL
  284.             IF (NOTWRT) GO TO 270
  285.                WRITE (LUNIT,390)
  286.                CALL SARRAY(A,LDA,P,P,P,LUNIT)
  287.   270       CONTINUE
  288.   280    CONTINUE
  289.   290 CONTINUE
  290.       WRITE (LUNIT,440)
  291.       RETURN
  292.   300 FORMAT (22H1LINPACK TESTER, SCH** /
  293.      *        29H THIS VERSION DATED 08/14/78.)
  294.   310 FORMAT ( // 5H1CASE, I3)
  295.   320 FORMAT ( / 19H 5 X 5 NO PIVOTING.)
  296.   330 FORMAT ( / 22H MONOELEMENTAL MATRIX.)
  297.   340 FORMAT ( / 24H 8 X 8 PIVOT LOGIC TEST.)
  298.   350 FORMAT ( / 32H 6 X 6 NEGATIVE EIGENVALUE TEST.)
  299.   360 FORMAT ( / 16H 25 X 25 MATRIX.)
  300.   370 FORMAT ( / 32H 5 X 5 PIVOT AND UNDERFLOW TEST.)
  301.   380 FORMAT ( / 31H 5 X 5 PIVOT AND OVERFLOW TEST.)
  302.   390 FORMAT ( / 2H A)
  303.   400 FORMAT ( / 5H JPVT // (1H , 10I5))
  304.   410 FORMAT ( // 18H A - TRANS(R)*R = , 1PE16.8)
  305.   420 FORMAT ( / 39H JOB AND JPVT BEFORE THE DECOMPOSITION. / I5 /
  306.      *         (10I5))
  307.   430 FORMAT ( / 20H THE VALUE OF  KD  =, I5)
  308.   440 FORMAT ( /// 12H END OF TEST / 1H1)
  309.       END
  310.       SUBROUTINE SCHXX(A,LDA,P,AA,JPVT,KPVT,ASAVE,RESDUL,JOB)
  311. C
  312.       INTEGER LDA,P,JPVT(1),KPVT(1),JOB
  313.       REAL A(LDA,1),AA(LDA,1),ASAVE(LDA,1)
  314.       REAL RESDUL,ANORM,SASUM
  315. C
  316.       REAL TEMP,SDOT
  317. C
  318. C     FORM TRANS(R)*R
  319. C
  320.       DO 20 I = 1, P
  321.          DO 10 J = I, P
  322.             AA(I,J) = SDOT(I,A(1,I),1,A(1,J),1)
  323.             AA(J,I) = AA(I,J)
  324.    10    CONTINUE
  325.          KPVT(I) = JPVT(I)
  326.    20 CONTINUE
  327.       IF (JOB .EQ. 0) GO TO 70
  328. C
  329. C        UNSCRAMBLE TRANS(R)*R
  330. C
  331.          DO 60 J = 1, P
  332.    30       IF (KPVT(J) .EQ. J) GO TO 50
  333.                IK = KPVT(J)
  334.                CALL SSWAP(P,AA(1,J),1,AA(1,IK),1)
  335.                DO 40 I = 1, P
  336.                   TEMP = AA(J,I)
  337.                   AA(J,I) = AA(IK,I)
  338.                   AA(IK,I) = TEMP
  339.    40          CONTINUE
  340.                KPVT(J) = KPVT(IK)
  341.                KPVT(IK) = IK
  342.             GO TO 30
  343.    50       CONTINUE
  344.    60    CONTINUE
  345.    70 CONTINUE
  346.       ANORM = 0.0E0
  347.       RESDUL = 0.0E0
  348.       DO 90 J = 1, P
  349.          ANORM = AMAX1(ANORM,SASUM(P,ASAVE(1,J),1))
  350.          DO 80 I = 1, P
  351.             ASAVE(I,J) = ASAVE(I,J) - AA(I,J)
  352.    80    CONTINUE
  353.          RESDUL = AMAX1(RESDUL,SASUM(P,ASAVE(1,J),1))
  354.    90 CONTINUE
  355.       RESDUL = RESDUL/ANORM
  356.       RETURN
  357.       END
  358.       SUBROUTINE SARRAY(A,LDA,M,N,NNL,LUNIT)
  359.       INTEGER LDA,M,N,NNL,LUNIT
  360.       REAL A(LDA,1)
  361. C
  362. C     FORTRAN IABS,MIN0
  363. C
  364.       INTEGER I,J,K,KU,NL
  365.       NL = IABS(NNL)
  366.       IF (NNL .LT. 0) GO TO 30
  367.          DO 20 I = 1, M
  368.             WRITE (LUNIT,70)
  369.             DO 10 K = 1, N, NL
  370.                KU = MIN0(K+NL-1,N)
  371.                WRITE (LUNIT,70) (A(I,J), J = K, KU)
  372.    10       CONTINUE
  373.    20    CONTINUE
  374.       GO TO 60
  375.    30 CONTINUE
  376.          DO 50 J = 1, N
  377.             WRITE (LUNIT,70)
  378.             DO 40 K = 1, M, NL
  379.                KU = MIN0(K+NL-1,M)
  380.                WRITE (LUNIT,70) (A(I,J), I = K, KU)
  381.    40       CONTINUE
  382.    50    CONTINUE
  383.    60 CONTINUE
  384.       RETURN
  385.    70 FORMAT (1H , 8E13.6)
  386.       END
  387.       SUBROUTINE SAGEN(A,LDA,P,D,ASAVE)
  388.       INTEGER LDA,P
  389.       REAL A(LDA,1),D(1),ASAVE(LDA,1)
  390. C
  391.       REAL EN,ES,FDPTP,S,TDP
  392. C
  393.       S = 0.0E0
  394.       DO 10 I = 1, P
  395.          S = S + D(I)
  396.    10 CONTINUE
  397.       TDP = 2.0E0/FLOAT(P)
  398.       FDPTP = 4.0E0*S/FLOAT(P*P)
  399.       DO 30 I = 1, P
  400.          ES = -1.0E0
  401.          DO 20 J = I, P
  402.             ES = (-1.0E0)*ES
  403.             A(I,J) = -ES*TDP*(D(J) + D(I)) + ES*FDPTP
  404.             A(I,J) = A(I,J)*FLOAT(I+1)*FLOAT(J+1)/FLOAT((J+1)**2+J**2)
  405.             ASAVE(I,J) = A(I,J)
  406.             A(J,I) = A(I,J)
  407.             ASAVE(J,I) = A(J,I)
  408.    20    CONTINUE
  409.          A(I,I) = A(I,I) + D(I)
  410.          ASAVE(I,I) = A(I,I)
  411.    30 CONTINUE
  412.       RETURN
  413.       END
  414.