home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / linpack / sqrts.for < prev    next >
Text File  |  1984-01-06  |  15KB  |  454 lines
  1.       SUBROUTINE SQRTS(LUNIT)
  2. C     LUNIT IS THE OUTPUT UNIT NUMBER
  3. C
  4. C     TESTS
  5. C        SQRDC,SQRSL
  6. C
  7. C     LINPACK. THIS VERSION DATED 08/14/78 .
  8. C     G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
  9. C
  10. C     SUBROUTINES AND FUNCTIONS
  11. C
  12. C     LINPACK SQRDC,SQRSL
  13. C     EXTERNAL SQRBM,SQRFM,SQRRX,SXGEN,SMACH
  14. C     FORTRAN AMAX1,ABS,FLOAT
  15. C
  16. C     INTERNAL VARIABLES
  17. C
  18.       INTEGER LUNIT
  19.       INTEGER CASE,LDX,N,P,PD2,PD2P1,NCASE,JPVT(25)
  20.       INTEGER I,INFO,J,JJ,JJJ,L
  21.       REAL BESTAT,RMAX,RSTAT,XBMAX,XBSTAT,ERRLVL,BSTAT,FSTAT
  22.       REAL BETA(25),QRAUX(25),QY(25),RSD(25),S(25),WORK(25)
  23.       REAL X(25,25),XBETA(25),XX(25,25),Y(25),Y1(25),Y2(25)
  24.       REAL SMACH,TINY,HUGE
  25.       LOGICAL NOTWRT
  26.       LDX = 25
  27.       NCASE = 9
  28.       ERRLVL = 100.0E0
  29.       NOTWRT = .TRUE.
  30.       TINY = SMACH(2)
  31.       HUGE = SMACH(3)
  32.       WRITE (LUNIT,600)
  33.       DO 550 CASE = NCASE, NCASE
  34.          WRITE (LUNIT,560) CASE
  35.          XBSTAT = 0.0E0
  36.          BESTAT = 0.0E0
  37.          GO TO (10, 100, 170, 240, 300, 330, 380, 430, 470), CASE
  38. C
  39. C        WELL CONDITIONED LEAST SQUARES PROBLEM.
  40. C
  41.    10    CONTINUE
  42.             WRITE (LUNIT,640)
  43.             N = 10
  44.             P = 4
  45. C
  46. C           GENERATE A MATRIX X WITH SINGULAR VALUES 1. AND .5.
  47. C
  48.             PD2 = P/2
  49.             PD2P1 = PD2 + 1
  50.             DO 20 L = 1, PD2
  51.                S(L) = 1.0E0
  52.    20       CONTINUE
  53.             DO 30 L = PD2P1, P
  54.                S(L) = 0.5E0
  55.    30       CONTINUE
  56.             CALL SXGEN(X,LDX,N,P,S)
  57.             IF (NOTWRT) GO TO 40
  58.                WRITE (LUNIT,610)
  59.                CALL SARRAY(X,LDX,N,P,4,LUNIT)
  60.    40       CONTINUE
  61. C           THE OBSERVATION VECTOR IS Y = Y1 + Y2, WHERE Y1 IS
  62. C           THE VECTOR OF ROW SUMS OF X AND Y2 IS ORTHOGONAL TO
  63. C           THE COLUMNS OF X.
  64. C
  65.             DO 60 I = 1, N
  66.                Y1(I) = 0.0E0
  67.                Y2(I) = 2.0E0*TINY
  68.                IF (I .EQ. P + 1) Y2(I) = Y2(I) - FLOAT(N)*TINY
  69.                DO 50 J = 1, P
  70.                   X(I,J) = X(I,J)*TINY
  71.                   Y1(I) = Y1(I) + X(I,J)
  72.                   XX(I,J) = X(I,J)
  73.    50          CONTINUE
  74.                Y(I) = Y1(I) + Y2(I)
  75.    60       CONTINUE
  76. C
  77. C           REDUCE THE MATRIX.
  78. C
  79.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,0)
  80.             IF (NOTWRT) GO TO 70
  81.                WRITE (LUNIT,610)
  82.                CALL SARRAY(X,LDX,N,P,4,LUNIT)
  83.                WRITE (LUNIT,620)
  84.                CALL SARRAY(QRAUX,P,P,1,-4,LUNIT)
  85.    70       CONTINUE
  86. C
  87. C           COMPUTE THE STATISTICS FOR THE FOWARD AND BACK
  88. C           MULTIPLICATIONS.
  89. C
  90.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  91.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  92. C
  93. C           SOLVE THE LEAST SQUARES PROBLEM.
  94. C
  95.             CALL SQRSL(X,LDX,N,P,QRAUX,Y,QY,RSD,BETA,RSD,XBETA,111,
  96.      *                 INFO)
  97. C
  98. C           COMPUTE THE LEAST SQUARES TEST STATISTICS.
  99. C
  100.             RSTAT = 0.0E0
  101.             RMAX = 0.0E0
  102.             XBSTAT = 0.0E0
  103.             XBMAX = 0.0E0
  104.             DO 80 I = 1, N
  105.                RSTAT = AMAX1(RSTAT,ABS(Y2(I)-RSD(I)))
  106.                RMAX = AMAX1(RMAX,ABS(Y2(I)))
  107.                XBSTAT = AMAX1(XBSTAT,ABS(Y1(I)-XBETA(I)))
  108.                XBMAX = AMAX1(XBMAX,ABS(Y2(I)))
  109.    80       CONTINUE
  110.             RSTAT = RSTAT/(SMACH(1)*RMAX)
  111.             XBSTAT = XBSTAT/(SMACH(1)*XBMAX)
  112.             BESTAT = 0.0E0
  113.             DO 90 J = 1, P
  114.                BESTAT = AMAX1(BESTAT,ABS(1.0E0-BETA(J)))
  115.    90       CONTINUE
  116.             BESTAT = BESTAT/SMACH(1)
  117.             WRITE (LUNIT,570) FSTAT,BSTAT,BESTAT,XBSTAT,RSTAT
  118.          GO TO 540
  119. C
  120. C        4 X 10 MATRIX.
  121. C
  122.   100    CONTINUE
  123.             WRITE (LUNIT,650)
  124.             N = 4
  125.             P = 10
  126.             PD2 = P/2
  127.             PD2P1 = PD2 + 1
  128.             DO 110 L = 1, PD2
  129.                S(L) = 1.0E0
  130.   110       CONTINUE
  131.             DO 120 L = PD2P1, P
  132.                S(L) = 0.5E0
  133.   120       CONTINUE
  134.             CALL SXGEN(X,LDX,N,P,S)
  135.             IF (NOTWRT) GO TO 130
  136.                WRITE (LUNIT,610)
  137.                CALL SARRAY(X,LDX,N,P,10,LUNIT)
  138.   130       CONTINUE
  139.             DO 150 I = 1, N
  140.                DO 140 J = 1, P
  141.                   XX(I,J) = X(I,J)
  142.   140          CONTINUE
  143.   150       CONTINUE
  144.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,0)
  145.             IF (NOTWRT) GO TO 160
  146.                WRITE (LUNIT,610)
  147.                CALL SARRAY(X,LDX,N,P,10,LUNIT)
  148.                WRITE (LUNIT,620)
  149.                CALL SARRAY(QRAUX,N,N,1,-4,LUNIT)
  150.   160       CONTINUE
  151.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  152.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  153.             WRITE (LUNIT,570) FSTAT,BSTAT
  154.          GO TO 540
  155. C
  156. C        PIVOTING AND OVERFLOW TEST.  COLUMNS 1, 4, 7 FROZEN.
  157. C
  158.   170    CONTINUE
  159.             WRITE (LUNIT,670)
  160.             N = 10
  161.             P = 3
  162.             S(1) = 1.0E0
  163.             S(2) = 0.5E0
  164.             S(3) = 0.25E0
  165.             CALL SXGEN(X,LDX,N,P,S)
  166.             P = 9
  167.             DO 190 I = 1, N
  168.                DO 180 JJJ = 1, 3
  169.                   J = 4 - JJJ
  170.                   JJ = 3*J
  171.                   X(I,JJ) = HUGE*X(I,J)
  172.                   X(I,JJ-1) = X(I,JJ)/2.0E0
  173.                   X(I,JJ-2) = X(I,JJ-1)/2.0E0
  174.   180          CONTINUE
  175.   190       CONTINUE
  176.             IF (NOTWRT) GO TO 200
  177.                WRITE (LUNIT,610)
  178.                CALL SARRAY(X,LDX,N,P,9,LUNIT)
  179.   200       CONTINUE
  180.             DO 220 J = 1, P
  181.                JPVT(J) = 0
  182.                DO 210 I = 1, N
  183.                   XX(I,J) = X(I,J)
  184.   210          CONTINUE
  185.   220       CONTINUE
  186.             JPVT(1) = -1
  187.             JPVT(4) = -1
  188.             JPVT(7) = -1
  189.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
  190.             IF (NOTWRT) GO TO 230
  191.                WRITE (LUNIT,610)
  192.                CALL SARRAY(X,LDX,N,P,9,LUNIT)
  193.                WRITE (LUNIT,620)
  194.                CALL SARRAY(QRAUX,P,P,1,-9,LUNIT)
  195.   230       CONTINUE
  196.             WRITE (LUNIT,630) (JPVT(J), J = 1, P)
  197.             CALL SQRRX(XX,LDX,N,P,JPVT)
  198.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  199.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  200.             WRITE (LUNIT,570) FSTAT,BSTAT
  201.          GO TO 540
  202. C
  203. C        25 BY 25 MATRIX
  204. C
  205.   240    CONTINUE
  206.             WRITE (LUNIT,680)
  207.             N = 25
  208.             P = 25
  209.             DO 250 I = 1, 25
  210.                S(I) = 2.0E0**(1 - I)
  211.   250       CONTINUE
  212.             CALL SXGEN(X,LDX,N,P,S)
  213.             IF (NOTWRT) GO TO 260
  214.                WRITE (LUNIT,610)
  215.                CALL SARRAY(X,LDX,N,P,10,LUNIT)
  216.   260       CONTINUE
  217.             DO 280 J = 1, P
  218.                JPVT(J) = 0
  219.                DO 270 I = 1, N
  220.                   XX(I,J) = X(I,J)
  221.   270          CONTINUE
  222.   280       CONTINUE
  223.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
  224.             IF (NOTWRT) GO TO 290
  225.                WRITE (LUNIT,610)
  226.                CALL SARRAY(X,LDX,N,P,10,LUNIT)
  227.                WRITE (LUNIT,620)
  228.                CALL SARRAY(QRAUX,P,P,1,-10,LUNIT)
  229.   290       CONTINUE
  230.             WRITE (LUNIT,630) (JPVT(J), J = 1, P)
  231.             CALL SQRRX(XX,LDX,N,P,JPVT)
  232.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  233.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  234.             WRITE (LUNIT,570) FSTAT,BSTAT
  235.          GO TO 540
  236. C
  237. C        MONOELEMENTAL MATRIX.
  238. C
  239.   300    CONTINUE
  240.             WRITE (LUNIT,690)
  241.             N = 1
  242.             P = 1
  243.             X(1,1) = 1.0E0
  244.             IF (NOTWRT) GO TO 310
  245.                WRITE (LUNIT,610)
  246.                CALL SARRAY(X,LDX,N,P,1,LUNIT)
  247.   310       CONTINUE
  248.             XX(1,1) = 1.0E0
  249.             JPVT(1) = 0
  250.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
  251.             IF (NOTWRT) GO TO 320
  252.                WRITE (LUNIT,610)
  253.                CALL SARRAY(X,LDX,N,P,1,LUNIT)
  254.                WRITE (LUNIT,620)
  255.                CALL SARRAY(QRAUX,P,P,1,-1,LUNIT)
  256.   320       CONTINUE
  257.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  258.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  259.             WRITE (LUNIT,570) FSTAT,BSTAT
  260.          GO TO 540
  261. C
  262. C        ZERO MATRIX.
  263. C
  264.   330    CONTINUE
  265.             WRITE (LUNIT,700)
  266.             N = 10
  267.             P = 4
  268.             DO 350 J = 1, P
  269.                JPVT(J) = 0
  270.                DO 340 I = 1, N
  271.                   X(I,J) = 0.0E0
  272.                   XX(I,J) = 0.0E0
  273.   340          CONTINUE
  274.   350       CONTINUE
  275.             IF (NOTWRT) GO TO 360
  276.                WRITE (LUNIT,610)
  277.                CALL SARRAY(X,LDX,N,P,4,LUNIT)
  278.   360       CONTINUE
  279.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
  280.             IF (NOTWRT) GO TO 370
  281.                WRITE (LUNIT,610)
  282.                CALL SARRAY(X,LDX,N,P,4,LUNIT)
  283.                WRITE (LUNIT,620)
  284.                CALL SARRAY(QRAUX,P,P,1,-4,LUNIT)
  285.   370       CONTINUE
  286.             WRITE (LUNIT,630) (JPVT(J), J = 1, P)
  287.             CALL SQRRX(XX,LDX,N,P,JPVT)
  288.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  289.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  290.             WRITE (LUNIT,570) FSTAT,BSTAT
  291.          GO TO 540
  292. C
  293. C        10 X 1 MATRIX WITH LEAST SQUARES PROBLEM.
  294. C
  295.   380    CONTINUE
  296.             WRITE (LUNIT,710)
  297.             N = 10
  298.             P = 1
  299.             S(1) = 1.0E0
  300.             CALL SXGEN(X,LDX,N,P,S)
  301.             IF (NOTWRT) GO TO 390
  302.                WRITE (LUNIT,610)
  303.                CALL SARRAY(X,LDX,N,P,1,LUNIT)
  304.   390       CONTINUE
  305.             DO 400 I = 1, N
  306.                Y1(I) = 0.0E0
  307.                Y2(I) = 2.0E0
  308.                IF (I .EQ. P + 1) Y2(I) = Y2(I) - FLOAT(N)
  309.                Y1(I) = Y1(I) + X(I,1)
  310.                Y(I) = Y1(I) + Y2(I)
  311.                XX(I,1) = X(I,1)
  312.   400       CONTINUE
  313.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,0)
  314.             IF (NOTWRT) GO TO 410
  315.                WRITE (LUNIT,610)
  316.                CALL SARRAY(X,LDX,N,P,1,LUNIT)
  317.                WRITE (LUNIT,620)
  318.                CALL SARRAY(QRAUX,P,P,1,1,LUNIT)
  319.   410       CONTINUE
  320.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  321.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  322.             CALL SQRSL(X,LDX,N,P,QRAUX,Y,QY,RSD,BETA,RSD,XBETA,111,
  323.      *                 INFO)
  324.             RSTAT = 0.0E0
  325.             RMAX = 0.0E0
  326.             XBSTAT = 0.0E0
  327.             XBMAX = 0.0E0
  328.             DO 420 I = 1, N
  329.                RSTAT = AMAX1(RSTAT,ABS(Y2(I)-RSD(I)))
  330.                RMAX = AMAX1(RMAX,ABS(Y2(I)))
  331.                XBMAX = AMAX1(XBMAX,ABS(Y2(I)))
  332.                XBSTAT = AMAX1(XBSTAT,ABS(Y1(I)-XBETA(I)))
  333.   420       CONTINUE
  334.             RSTAT = RSTAT/(SMACH(1)*RMAX)
  335.             XBSTAT = XBSTAT/(SMACH(1)*XBMAX)
  336.             BESTAT = ABS(1.0E0-BETA(1))
  337.             BESTAT = BESTAT/SMACH(1)
  338.             WRITE (LUNIT,570) FSTAT,BSTAT,BESTAT,XBSTAT,RSTAT
  339.          GO TO 540
  340. C
  341. C        1 X 4 MATRIX
  342. C
  343.   430    CONTINUE
  344.             WRITE (LUNIT,720)
  345.             N = 1
  346.             P = 4
  347.             X(1,1) = 1.0E0
  348.             X(1,2) = 2.0E0
  349.             X(1,3) = 0.25E0
  350.             X(1,4) = 0.5E0
  351.             DO 440 I = 1, 4
  352.                JPVT(I) = 0
  353.                XX(1,I) = X(1,I)
  354.   440       CONTINUE
  355.             IF (NOTWRT) GO TO 450
  356.                WRITE (LUNIT,610)
  357.                CALL SARRAY(X,LDX,N,P,4,LUNIT)
  358.   450       CONTINUE
  359.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
  360.             IF (NOTWRT) GO TO 460
  361.                WRITE (LUNIT,610)
  362.                CALL SARRAY(X,LDX,N,P,10,LUNIT)
  363.                WRITE (LUNIT,620)
  364.                CALL SARRAY(QRAUX,N,N,1,-1,LUNIT)
  365.   460       CONTINUE
  366.             WRITE (LUNIT,630) (JPVT(J), J = 1, P)
  367.             CALL SQRRX(XX,LDX,N,P,JPVT)
  368.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  369.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  370.             WRITE (LUNIT,570) FSTAT,BSTAT
  371.          GO TO 540
  372. C
  373. C        PIVOTING TEST.
  374. C
  375.   470    CONTINUE
  376.             WRITE (LUNIT,660)
  377.             N = 10
  378.             P = 3
  379.             S(1) = 1.0E0
  380.             S(2) = 0.5E0
  381.             S(3) = 0.25E0
  382.             CALL SXGEN(X,LDX,N,P,S)
  383.             P = 9
  384.             DO 490 I = 1, N
  385.                DO 480 JJJ = 1, 3
  386.                   J = 4 - JJJ
  387.                   JJ = 3*J
  388.                   X(I,JJ) = X(I,J)
  389.                   X(I,JJ-1) = X(I,JJ)/2.0E0
  390.                   X(I,JJ-2) = X(I,JJ-1)/2.0E0
  391.   480          CONTINUE
  392.   490       CONTINUE
  393.             IF (NOTWRT) GO TO 500
  394.                WRITE (LUNIT,610)
  395.                CALL SARRAY(X,LDX,N,P,9,LUNIT)
  396.   500       CONTINUE
  397.             DO 520 J = 1, P
  398.                JPVT(J) = 0
  399.                DO 510 I = 1, N
  400.                   XX(I,J) = X(I,J)
  401.   510          CONTINUE
  402.   520       CONTINUE
  403.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
  404.             IF (NOTWRT) GO TO 530
  405.                WRITE (LUNIT,610)
  406.                CALL SARRAY(X,LDX,N,P,9,LUNIT)
  407.                WRITE (LUNIT,620)
  408.                CALL SARRAY(QRAUX,P,P,1,-9,LUNIT)
  409.   530       CONTINUE
  410.             WRITE (LUNIT,630) (JPVT(J), J = 1, P)
  411.             CALL SQRRX(XX,LDX,N,P,JPVT)
  412.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  413.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  414.             WRITE (LUNIT,570) FSTAT,BSTAT
  415.   540    CONTINUE
  416.          IF (AMAX1(FSTAT,BSTAT,BESTAT,XBSTAT,RSTAT) .GE. ERRLVL)
  417.      *      WRITE (LUNIT,580)
  418.   550 CONTINUE
  419.       WRITE (LUNIT,590)
  420.       RETURN
  421. C
  422.   560 FORMAT ( // '1CASE', I3 )
  423.   570 FORMAT ( / 11H STATISTICS //
  424.      *         35H    FORWARD MULTIPLICATION ........, E10.2 /
  425.      *         35H    BACK MULTIPLICATION ..........., E10.2 /
  426.      *         35H    BETA .........................., E10.2 /
  427.      *         35H    X*BETA ........................, E10.2 /
  428.      *         35H    RESIDUAL ......................, E10.2)
  429.   580 FORMAT ( / 34H *****STATISTICS ABOVE ERROR LEVEL)
  430.   590 FORMAT ( / 15H1END OF QR TEST)
  431.   600 FORMAT ('1LINPACK TESTER, SQR**' /
  432.      *        ' THIS VERSION DATED 08/14/78.')
  433.   610 FORMAT ( / 2H X)
  434.   620 FORMAT ( / 6H QRAUX)
  435.   630 FORMAT ( / 5H JPVT // (1H , 10I5))
  436.   640 FORMAT ( / 39H WELL CONDITIONED LEAST SQUARES PROBLEM /
  437.      *         20H AND UNDERFLOW TEST.)
  438.   650 FORMAT ( / 14H 4 X 10 MATRIX)
  439.   660 FORMAT ( / 14H PIVOTING TEST /
  440.      *         42H ON RETURN THE FIRST THREE ENTRIES OF JPVT /
  441.      *         36H SHOULD BE 3,6,9 BUT NOT NECESSARILY /
  442.      *         15H IN THAT ORDER.)
  443.   670 FORMAT ( / 27H PIVOTING AND OVERFLOW TEST /
  444.      *         26H WITH COLUMNS 1,4,7 FROZEN /
  445.      *         42H ON RETURN THE LAST  THREE ENTRIES OF JPVT /
  446.      *         31H SHOULD BE 1,4,7 IN THAT ORDER.)
  447.   680 FORMAT ( / 15H 25 X 25 MATRIX)
  448.   690 FORMAT ( / 21H MONOELEMENTAL MATRIX)
  449.   700 FORMAT ( / 12H ZERO MATRIX)
  450.   710 FORMAT ( / 41H 10 X 1 MATRIX WITH LEAST SQUARES PROBLEM)
  451.   720 FORMAT ( / 13H 1 X 4 MATRIX)
  452.       END
  453.  
  454.