home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / linpack / sqr.for < prev    next >
Encoding:
Text File  |  1984-01-06  |  19.9 KB  |  665 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 SQRTS(LUNIT)
  11.       STOP
  12.       END
  13.       SUBROUTINE SQRTS(LUNIT)
  14. C     LUNIT IS THE OUTPUT UNIT NUMBER
  15. C
  16. C     TESTS
  17. C        SQRDC,SQRSL
  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     LINPACK SQRDC,SQRSL
  25. C     EXTERNAL SQRBM,SQRFM,SQRRX,SXGEN,SMACH
  26. C     FORTRAN AMAX1,ABS,FLOAT
  27. C
  28. C     INTERNAL VARIABLES
  29. C
  30.       INTEGER LUNIT
  31.       INTEGER CASE,LDX,N,P,PD2,PD2P1,NCASE,JPVT(25)
  32.       INTEGER I,INFO,J,JJ,JJJ,L
  33.       REAL BESTAT,RMAX,RSTAT,XBMAX,XBSTAT,ERRLVL,BSTAT,FSTAT
  34.       REAL BETA(25),QRAUX(25),QY(25),RSD(25),S(25),WORK(25)
  35.       REAL X(25,25),XBETA(25),XX(25,25),Y(25),Y1(25),Y2(25)
  36.       REAL SMACH,TINY,HUGE
  37.       LOGICAL NOTWRT
  38.       LDX = 25
  39.       NCASE = 9
  40.       ERRLVL = 100.0E0
  41.       NOTWRT = .TRUE.
  42.       TINY = SMACH(2)
  43.       HUGE = SMACH(3)
  44.       WRITE (LUNIT,600)
  45.       DO 550 CASE = 1, NCASE
  46.          WRITE (LUNIT,560) CASE
  47.          XBSTAT = 0.0E0
  48.          BESTAT = 0.0E0
  49.          GO TO (10, 100, 170, 240, 300, 330, 380, 430, 470), CASE
  50. C
  51. C        WELL CONDITIONED LEAST SQUARES PROBLEM.
  52. C
  53.    10    CONTINUE
  54.             WRITE (LUNIT,640)
  55.             N = 10
  56.             P = 4
  57. C
  58. C           GENERATE A MATRIX X WITH SINGULAR VALUES 1. AND .5.
  59. C
  60.             PD2 = P/2
  61.             PD2P1 = PD2 + 1
  62.             DO 20 L = 1, PD2
  63.                S(L) = 1.0E0
  64.    20       CONTINUE
  65.             DO 30 L = PD2P1, P
  66.                S(L) = 0.5E0
  67.    30       CONTINUE
  68.             CALL SXGEN(X,LDX,N,P,S)
  69.             IF (NOTWRT) GO TO 40
  70.                WRITE (LUNIT,610)
  71.                CALL SARRAY(X,LDX,N,P,4,LUNIT)
  72.    40       CONTINUE
  73. C           THE OBSERVATION VECTOR IS Y = Y1 + Y2, WHERE Y1 IS
  74. C           THE VECTOR OF ROW SUMS OF X AND Y2 IS ORTHOGONAL TO
  75. C           THE COLUMNS OF X.
  76. C
  77.             DO 60 I = 1, N
  78.                Y1(I) = 0.0E0
  79.                Y2(I) = 2.0E0*TINY
  80.                IF (I .EQ. P + 1) Y2(I) = Y2(I) - FLOAT(N)*TINY
  81.                DO 50 J = 1, P
  82.                   X(I,J) = X(I,J)*TINY
  83.                   Y1(I) = Y1(I) + X(I,J)
  84.                   XX(I,J) = X(I,J)
  85.    50          CONTINUE
  86.                Y(I) = Y1(I) + Y2(I)
  87.    60       CONTINUE
  88. C
  89. C           REDUCE THE MATRIX.
  90. C
  91.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,0)
  92.             IF (NOTWRT) GO TO 70
  93.                WRITE (LUNIT,610)
  94.                CALL SARRAY(X,LDX,N,P,4,LUNIT)
  95.                WRITE (LUNIT,620)
  96.                CALL SARRAY(QRAUX,P,P,1,-4,LUNIT)
  97.    70       CONTINUE
  98. C
  99. C           COMPUTE THE STATISTICS FOR THE FOWARD AND BACK
  100. C           MULTIPLICATIONS.
  101. C
  102.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  103.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  104. C
  105. C           SOLVE THE LEAST SQUARES PROBLEM.
  106. C
  107.             CALL SQRSL(X,LDX,N,P,QRAUX,Y,QY,RSD,BETA,RSD,XBETA,111,
  108.      *                 INFO)
  109. C
  110. C           COMPUTE THE LEAST SQUARES TEST STATISTICS.
  111. C
  112.             RSTAT = 0.0E0
  113.             RMAX = 0.0E0
  114.             XBSTAT = 0.0E0
  115.             XBMAX = 0.0E0
  116.             DO 80 I = 1, N
  117.                RSTAT = AMAX1(RSTAT,ABS(Y2(I)-RSD(I)))
  118.                RMAX = AMAX1(RMAX,ABS(Y2(I)))
  119.                XBSTAT = AMAX1(XBSTAT,ABS(Y1(I)-XBETA(I)))
  120.                XBMAX = AMAX1(XBMAX,ABS(Y2(I)))
  121.    80       CONTINUE
  122.             RSTAT = RSTAT/(SMACH(1)*RMAX)
  123.             XBSTAT = XBSTAT/(SMACH(1)*XBMAX)
  124.             BESTAT = 0.0E0
  125.             DO 90 J = 1, P
  126.                BESTAT = AMAX1(BESTAT,ABS(1.0E0-BETA(J)))
  127.    90       CONTINUE
  128.             BESTAT = BESTAT/SMACH(1)
  129.             WRITE (LUNIT,570) FSTAT,BSTAT,BESTAT,XBSTAT,RSTAT
  130.          GO TO 540
  131. C
  132. C        4 X 10 MATRIX.
  133. C
  134.   100    CONTINUE
  135.             WRITE (LUNIT,650)
  136.             N = 4
  137.             P = 10
  138.             PD2 = P/2
  139.             PD2P1 = PD2 + 1
  140.             DO 110 L = 1, PD2
  141.                S(L) = 1.0E0
  142.   110       CONTINUE
  143.             DO 120 L = PD2P1, P
  144.                S(L) = 0.5E0
  145.   120       CONTINUE
  146.             CALL SXGEN(X,LDX,N,P,S)
  147.             IF (NOTWRT) GO TO 130
  148.                WRITE (LUNIT,610)
  149.                CALL SARRAY(X,LDX,N,P,10,LUNIT)
  150.   130       CONTINUE
  151.             DO 150 I = 1, N
  152.                DO 140 J = 1, P
  153.                   XX(I,J) = X(I,J)
  154.   140          CONTINUE
  155.   150       CONTINUE
  156.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,0)
  157.             IF (NOTWRT) GO TO 160
  158.                WRITE (LUNIT,610)
  159.                CALL SARRAY(X,LDX,N,P,10,LUNIT)
  160.                WRITE (LUNIT,620)
  161.                CALL SARRAY(QRAUX,N,N,1,-4,LUNIT)
  162.   160       CONTINUE
  163.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  164.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  165.             WRITE (LUNIT,570) FSTAT,BSTAT
  166.          GO TO 540
  167. C
  168. C        PIVOTING AND OVERFLOW TEST.  COLUMNS 1, 4, 7 FROZEN.
  169. C
  170.   170    CONTINUE
  171.             WRITE (LUNIT,670)
  172.             N = 10
  173.             P = 3
  174.             S(1) = 1.0E0
  175.             S(2) = 0.5E0
  176.             S(3) = 0.25E0
  177.             CALL SXGEN(X,LDX,N,P,S)
  178.             P = 9
  179.             DO 190 I = 1, N
  180.                DO 180 JJJ = 1, 3
  181.                   J = 4 - JJJ
  182.                   JJ = 3*J
  183.                   X(I,JJ) = HUGE*X(I,J)
  184.                   X(I,JJ-1) = X(I,JJ)/2.0E0
  185.                   X(I,JJ-2) = X(I,JJ-1)/2.0E0
  186.   180          CONTINUE
  187.   190       CONTINUE
  188.             IF (NOTWRT) GO TO 200
  189.                WRITE (LUNIT,610)
  190.                CALL SARRAY(X,LDX,N,P,9,LUNIT)
  191.   200       CONTINUE
  192.             DO 220 J = 1, P
  193.                JPVT(J) = 0
  194.                DO 210 I = 1, N
  195.                   XX(I,J) = X(I,J)
  196.   210          CONTINUE
  197.   220       CONTINUE
  198.             JPVT(1) = -1
  199.             JPVT(4) = -1
  200.             JPVT(7) = -1
  201.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
  202.             IF (NOTWRT) GO TO 230
  203.                WRITE (LUNIT,610)
  204.                CALL SARRAY(X,LDX,N,P,9,LUNIT)
  205.                WRITE (LUNIT,620)
  206.                CALL SARRAY(QRAUX,P,P,1,-9,LUNIT)
  207.   230       CONTINUE
  208.             WRITE (LUNIT,630) (JPVT(J), J = 1, P)
  209.             CALL SQRRX(XX,LDX,N,P,JPVT)
  210.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  211.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  212.             WRITE (LUNIT,570) FSTAT,BSTAT
  213.          GO TO 540
  214. C
  215. C        25 BY 25 MATRIX
  216. C
  217.   240    CONTINUE
  218.             WRITE (LUNIT,680)
  219.             N = 25
  220.             P = 25
  221.             DO 250 I = 1, 25
  222.                S(I) = 2.0E0**(1 - I)
  223.   250       CONTINUE
  224.             CALL SXGEN(X,LDX,N,P,S)
  225.             IF (NOTWRT) GO TO 260
  226.                WRITE (LUNIT,610)
  227.                CALL SARRAY(X,LDX,N,P,10,LUNIT)
  228.   260       CONTINUE
  229.             DO 280 J = 1, P
  230.                JPVT(J) = 0
  231.                DO 270 I = 1, N
  232.                   XX(I,J) = X(I,J)
  233.   270          CONTINUE
  234.   280       CONTINUE
  235.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
  236.             IF (NOTWRT) GO TO 290
  237.                WRITE (LUNIT,610)
  238.                CALL SARRAY(X,LDX,N,P,10,LUNIT)
  239.                WRITE (LUNIT,620)
  240.                CALL SARRAY(QRAUX,P,P,1,-10,LUNIT)
  241.   290       CONTINUE
  242.             WRITE (LUNIT,630) (JPVT(J), J = 1, P)
  243.             CALL SQRRX(XX,LDX,N,P,JPVT)
  244.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  245.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  246.             WRITE (LUNIT,570) FSTAT,BSTAT
  247.          GO TO 540
  248. C
  249. C        MONOELEMENTAL MATRIX.
  250. C
  251.   300    CONTINUE
  252.             WRITE (LUNIT,690)
  253.             N = 1
  254.             P = 1
  255.             X(1,1) = 1.0E0
  256.             IF (NOTWRT) GO TO 310
  257.                WRITE (LUNIT,610)
  258.                CALL SARRAY(X,LDX,N,P,1,LUNIT)
  259.   310       CONTINUE
  260.             XX(1,1) = 1.0E0
  261.             JPVT(1) = 0
  262.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
  263.             IF (NOTWRT) GO TO 320
  264.                WRITE (LUNIT,610)
  265.                CALL SARRAY(X,LDX,N,P,1,LUNIT)
  266.                WRITE (LUNIT,620)
  267.                CALL SARRAY(QRAUX,P,P,1,-1,LUNIT)
  268.   320       CONTINUE
  269.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  270.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  271.             WRITE (LUNIT,570) FSTAT,BSTAT
  272.          GO TO 540
  273. C
  274. C        ZERO MATRIX.
  275. C
  276.   330    CONTINUE
  277.             WRITE (LUNIT,700)
  278.             N = 10
  279.             P = 4
  280.             DO 350 J = 1, P
  281.                JPVT(J) = 0
  282.                DO 340 I = 1, N
  283.                   X(I,J) = 0.0E0
  284.                   XX(I,J) = 0.0E0
  285.   340          CONTINUE
  286.   350       CONTINUE
  287.             IF (NOTWRT) GO TO 360
  288.                WRITE (LUNIT,610)
  289.                CALL SARRAY(X,LDX,N,P,4,LUNIT)
  290.   360       CONTINUE
  291.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
  292.             IF (NOTWRT) GO TO 370
  293.                WRITE (LUNIT,610)
  294.                CALL SARRAY(X,LDX,N,P,4,LUNIT)
  295.                WRITE (LUNIT,620)
  296.                CALL SARRAY(QRAUX,P,P,1,-4,LUNIT)
  297.   370       CONTINUE
  298.             WRITE (LUNIT,630) (JPVT(J), J = 1, P)
  299.             CALL SQRRX(XX,LDX,N,P,JPVT)
  300.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  301.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  302.             WRITE (LUNIT,570) FSTAT,BSTAT
  303.          GO TO 540
  304. C
  305. C        10 X 1 MATRIX WITH LEAST SQUARES PROBLEM.
  306. C
  307.   380    CONTINUE
  308.             WRITE (LUNIT,710)
  309.             N = 10
  310.             P = 1
  311.             S(1) = 1.0E0
  312.             CALL SXGEN(X,LDX,N,P,S)
  313.             IF (NOTWRT) GO TO 390
  314.                WRITE (LUNIT,610)
  315.                CALL SARRAY(X,LDX,N,P,1,LUNIT)
  316.   390       CONTINUE
  317.             DO 400 I = 1, N
  318.                Y1(I) = 0.0E0
  319.                Y2(I) = 2.0E0
  320.                IF (I .EQ. P + 1) Y2(I) = Y2(I) - FLOAT(N)
  321.                Y1(I) = Y1(I) + X(I,1)
  322.                Y(I) = Y1(I) + Y2(I)
  323.                XX(I,1) = X(I,1)
  324.   400       CONTINUE
  325.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,0)
  326.             IF (NOTWRT) GO TO 410
  327.                WRITE (LUNIT,610)
  328.                CALL SARRAY(X,LDX,N,P,1,LUNIT)
  329.                WRITE (LUNIT,620)
  330.                CALL SARRAY(QRAUX,P,P,1,1,LUNIT)
  331.   410       CONTINUE
  332.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  333.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  334.             CALL SQRSL(X,LDX,N,P,QRAUX,Y,QY,RSD,BETA,RSD,XBETA,111,
  335.      *                 INFO)
  336.             RSTAT = 0.0E0
  337.             RMAX = 0.0E0
  338.             XBSTAT = 0.0E0
  339.             XBMAX = 0.0E0
  340.             DO 420 I = 1, N
  341.                RSTAT = AMAX1(RSTAT,ABS(Y2(I)-RSD(I)))
  342.                RMAX = AMAX1(RMAX,ABS(Y2(I)))
  343.                XBMAX = AMAX1(XBMAX,ABS(Y2(I)))
  344.                XBSTAT = AMAX1(XBSTAT,ABS(Y1(I)-XBETA(I)))
  345.   420       CONTINUE
  346.             RSTAT = RSTAT/(SMACH(1)*RMAX)
  347.             XBSTAT = XBSTAT/(SMACH(1)*XBMAX)
  348.             BESTAT = ABS(1.0E0-BETA(1))
  349.             BESTAT = BESTAT/SMACH(1)
  350.             WRITE (LUNIT,570) FSTAT,BSTAT,BESTAT,XBSTAT,RSTAT
  351.          GO TO 540
  352. C
  353. C        1 X 4 MATRIX
  354. C
  355.   430    CONTINUE
  356.             WRITE (LUNIT,720)
  357.             N = 1
  358.             P = 4
  359.             X(1,1) = 1.0E0
  360.             X(1,2) = 2.0E0
  361.             X(1,3) = 0.25E0
  362.             X(1,4) = 0.5E0
  363.             DO 440 I = 1, 4
  364.                JPVT(I) = 0
  365.                XX(1,I) = X(1,I)
  366.   440       CONTINUE
  367.             IF (NOTWRT) GO TO 450
  368.                WRITE (LUNIT,610)
  369.                CALL SARRAY(X,LDX,N,P,4,LUNIT)
  370.   450       CONTINUE
  371.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
  372.             IF (NOTWRT) GO TO 460
  373.                WRITE (LUNIT,610)
  374.                CALL SARRAY(X,LDX,N,P,10,LUNIT)
  375.                WRITE (LUNIT,620)
  376.                CALL SARRAY(QRAUX,N,N,1,-1,LUNIT)
  377.   460       CONTINUE
  378.             WRITE (LUNIT,630) (JPVT(J), J = 1, P)
  379.             CALL SQRRX(XX,LDX,N,P,JPVT)
  380.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  381.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  382.             WRITE (LUNIT,570) FSTAT,BSTAT
  383.          GO TO 540
  384. C
  385. C        PIVOTING TEST.
  386. C
  387.   470    CONTINUE
  388.             WRITE (LUNIT,660)
  389.             N = 10
  390.             P = 3
  391.             S(1) = 1.0E0
  392.             S(2) = 0.5E0
  393.             S(3) = 0.25E0
  394.             CALL SXGEN(X,LDX,N,P,S)
  395.             P = 9
  396.             DO 490 I = 1, N
  397.                DO 480 JJJ = 1, 3
  398.                   J = 4 - JJJ
  399.                   JJ = 3*J
  400.                   X(I,JJ) = X(I,J)
  401.                   X(I,JJ-1) = X(I,JJ)/2.0E0
  402.                   X(I,JJ-2) = X(I,JJ-1)/2.0E0
  403.   480          CONTINUE
  404.   490       CONTINUE
  405.             IF (NOTWRT) GO TO 500
  406.                WRITE (LUNIT,610)
  407.                CALL SARRAY(X,LDX,N,P,9,LUNIT)
  408.   500       CONTINUE
  409.             DO 520 J = 1, P
  410.                JPVT(J) = 0
  411.                DO 510 I = 1, N
  412.                   XX(I,J) = X(I,J)
  413.   510          CONTINUE
  414.   520       CONTINUE
  415.             CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
  416.             IF (NOTWRT) GO TO 530
  417.                WRITE (LUNIT,610)
  418.                CALL SARRAY(X,LDX,N,P,9,LUNIT)
  419.                WRITE (LUNIT,620)
  420.                CALL SARRAY(QRAUX,P,P,1,-9,LUNIT)
  421.   530       CONTINUE
  422.             WRITE (LUNIT,630) (JPVT(J), J = 1, P)
  423.             CALL SQRRX(XX,LDX,N,P,JPVT)
  424.             CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
  425.             CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
  426.             WRITE (LUNIT,570) FSTAT,BSTAT
  427.   540    CONTINUE
  428.          IF (AMAX1(FSTAT,BSTAT,BESTAT,XBSTAT,RSTAT) .GE. ERRLVL)
  429.      *      WRITE (LUNIT,580)
  430.   550 CONTINUE
  431.       WRITE (LUNIT,590)
  432.       RETURN
  433. C
  434.   560 FORMAT ( // '1CASE', I3 )
  435.   570 FORMAT ( / 11H STATISTICS //
  436.      *         35H    FORWARD MULTIPLICATION ........, E10.2 /
  437.      *         35H    BACK MULTIPLICATION ..........., E10.2 /
  438.      *         35H    BETA .........................., E10.2 /
  439.      *         35H    X*BETA ........................, E10.2 /
  440.      *         35H    RESIDUAL ......................, E10.2)
  441.   580 FORMAT ( / 34H *****STATISTICS ABOVE ERROR LEVEL)
  442.   590 FORMAT ( / 15H1END OF QR TEST)
  443.   600 FORMAT ('1LINPACK TESTER, SQR**' /
  444.      *        ' THIS VERSION DATED 08/14/78.')
  445.   610 FORMAT ( / 2H X)
  446.   620 FORMAT ( / 6H QRAUX)
  447.   630 FORMAT ( / 5H JPVT // (1H , 10I5))
  448.   640 FORMAT ( / 39H WELL CONDITIONED LEAST SQUARES PROBLEM /
  449.      *         20H AND UNDERFLOW TEST.)
  450.   650 FORMAT ( / 14H 4 X 10 MATRIX)
  451.   660 FORMAT ( / 14H PIVOTING TEST /
  452.      *         42H ON RETURN THE FIRST THREE ENTRIES OF JPVT /
  453.      *         36H SHOULD BE 3,6,9 BUT NOT NECESSARILY /
  454.      *         15H IN THAT ORDER.)
  455.   670 FORMAT ( / 27H PIVOTING AND OVERFLOW TEST /
  456.      *         26H WITH COLUMNS 1,4,7 FROZEN /
  457.      *         42H ON RETURN THE LAST  THREE ENTRIES OF JPVT /
  458.      *         31H SHOULD BE 1,4,7 IN THAT ORDER.)
  459.   680 FORMAT ( / 15H 25 X 25 MATRIX)
  460.   690 FORMAT ( / 21H MONOELEMENTAL MATRIX)
  461.   700 FORMAT ( / 12H ZERO MATRIX)
  462.   710 FORMAT ( / 41H 10 X 1 MATRIX WITH LEAST SQUARES PROBLEM)
  463.   720 FORMAT ( / 13H 1 X 4 MATRIX)
  464.       END
  465.       SUBROUTINE SARRAY(A,LDA,M,N,NNL,LUNIT)
  466.       INTEGER LDA,M,N,NNL,LUNIT
  467.       REAL A(LDA,1)
  468. C
  469. C     FORTRAN IABS,MIN0
  470. C
  471.       INTEGER I,J,K,KU,NL
  472.       NL = IABS(NNL)
  473.       IF (NNL .LT. 0) GO TO 30
  474.          DO 20 I = 1, M
  475.             WRITE (LUNIT,70)
  476.             DO 10 K = 1, N, NL
  477.                KU = MIN0(K+NL-1,N)
  478.                WRITE (LUNIT,70) (A(I,J), J = K, KU)
  479.    10       CONTINUE
  480.    20    CONTINUE
  481.       GO TO 60
  482.    30 CONTINUE
  483.          DO 50 J = 1, N
  484.             WRITE (LUNIT,70)
  485.             DO 40 K = 1, M, NL
  486.                KU = MIN0(K+NL-1,M)
  487.                WRITE (LUNIT,70) (A(I,J), I = K, KU)
  488.    40       CONTINUE
  489.    50    CONTINUE
  490.    60 CONTINUE
  491.       RETURN
  492.    70 FORMAT (1H , 8E13.6)
  493.       END
  494.       SUBROUTINE SQRBM(X,LDX,N,P,XX,QR,QRAUX,BSTAT)
  495.       INTEGER LDX,N,P
  496.       REAL BSTAT
  497.       REAL X(LDX,1),XX(LDX,1),QR(1),QRAUX(1)
  498. C
  499. C     SQRBM TAKES THE OUTPUT OF SQRDC AND MULTIPLIES THE FACTORS
  500. C     Q AND R TOGETHER.  THE RESULTS ARE COMPARED WITH XX,
  501. C     WHICH PRESUMABLY CONTAINS THE ORIGINAL MATRIX.
  502. C
  503. C     LINPACK SQRSL
  504. C     EXTERNAL SMACH
  505. C     FORTRAN AMAX1,ABS,MIN0
  506. C
  507.       INTEGER IU,JP1,I,INFO,J
  508.       REAL SMACH,EMAX,XMAX
  509.       EMAX = 0.0E0
  510.       XMAX = 0.0E0
  511.       DO 50 J = 1, P
  512.          IU = MIN0(N,J)
  513.          DO 10 I = 1, IU
  514.             QR(I) = X(I,J)
  515.    10    CONTINUE
  516.          JP1 = J + 1
  517.          IF (N .LT. JP1) GO TO 30
  518.          DO 20 I = JP1, N
  519.             QR(I) = 0.0E0
  520.    20    CONTINUE
  521.    30    CONTINUE
  522.          CALL SQRSL(X,LDX,N,P,QRAUX,QR,QR,QR,QR,QR,QR,10000,INFO)
  523.          DO 40 I = 1, N
  524.             EMAX = AMAX1(EMAX,ABS(XX(I,J)-QR(I)))
  525.             XMAX = AMAX1(XMAX,ABS(XX(I,J)))
  526.    40    CONTINUE
  527.    50 CONTINUE
  528.       IF (XMAX .NE. 0.0E0) GO TO 80
  529.          IF (EMAX .NE. 0.0E0) GO TO 60
  530.             BSTAT = 0.0E0
  531.          GO TO 70
  532.    60    CONTINUE
  533.             BSTAT = 1.0E20
  534.    70    CONTINUE
  535.       GO TO 90
  536.    80 CONTINUE
  537.          BSTAT = (EMAX/XMAX)/SMACH(1)
  538.    90 CONTINUE
  539.       RETURN
  540.       END
  541.       SUBROUTINE SQRFM(X,LDX,N,P,XX,QTX,QRAUX,FSTAT)
  542.       INTEGER LDX,N,P
  543.       REAL FSTAT
  544.       REAL X(LDX,1),XX(LDX,1),QTX(1),QRAUX(1)
  545. C
  546. C     SQRFM TAKES THE OUTPUT OF SQRDC AND APPLIES THE
  547. C     TRANSFORMATIONS TO XX, WHICH PRESUMABLY CONTAINS THE
  548. C     ORIGINAL MATRIX.  THE RESULTS ARE COMPARED WITH THE
  549. C     R PART OF THE QR DECOMPOSITION.
  550. C
  551. C     LINPACK SQRSL
  552. C     EXTERNAL SMACH
  553. C     FORTRAN AMAX1,ABS,MIN0
  554. C
  555.       INTEGER IU,JP1,I,INFO,J
  556.       REAL SMACH,EMAX,RMAX
  557.       EMAX = 0.0E0
  558.       RMAX = 0.0E0
  559.       DO 50 J = 1, P
  560.          DO 10 I = 1, N
  561.             QTX(I) = XX(I,J)
  562.    10    CONTINUE
  563.          CALL SQRSL(X,LDX,N,P,QRAUX,QTX,QTX,QTX,QTX,QTX,QTX,1000,INFO)
  564.          IU = MIN0(J,N)
  565.          DO 20 I = 1, IU
  566.             EMAX = AMAX1(EMAX,ABS(X(I,J)-QTX(I)))
  567.             RMAX = AMAX1(RMAX,ABS(X(I,J)))
  568.    20    CONTINUE
  569.          JP1 = J + 1
  570.          IF (N .LT. JP1) GO TO 40
  571.          DO 30 I = JP1, N
  572.             EMAX = AMAX1(EMAX,ABS(QTX(I)))
  573.    30    CONTINUE
  574.    40    CONTINUE
  575.    50 CONTINUE
  576.       IF (RMAX .NE. 0.0E0) GO TO 80
  577.          IF (EMAX .NE. 0.0E0) GO TO 60
  578.             FSTAT = 0.0E0
  579.          GO TO 70
  580.    60    CONTINUE
  581.             FSTAT = 1.0E20
  582.    70    CONTINUE
  583.       GO TO 90
  584.    80 CONTINUE
  585.          FSTAT = (EMAX/RMAX)/SMACH(1)
  586.    90 CONTINUE
  587.       RETURN
  588.       END
  589.       SUBROUTINE SQRRX(XX,LDX,N,P,JP)
  590.       INTEGER N,P,LDX
  591.       INTEGER JP(1)
  592.       REAL XX(LDX,1)
  593. C
  594. C     SQRRX REARRANGES THE COLUMNS OF XX IN THE ORDER SPECIFIED
  595. C     BY JPVT.
  596. C
  597. C     EXTERNAL SSWAP
  598. C
  599.       INTEGER I,IP,I1,J,JPVT(50),K
  600.       DO 10 J = 1, P
  601.          JPVT(J) = JP(J)
  602.    10 CONTINUE
  603.       IP = P - 1
  604.       DO 60 I = 1, IP
  605.          IF (JPVT(I) .EQ. I) GO TO 50
  606.             I1 = I + 1
  607.             DO 20 J = I1, P
  608. C           ...EXIT
  609.                IF (JPVT(I) .EQ. J) GO TO 30
  610.    20       CONTINUE
  611.    30       CONTINUE
  612.             CALL SSWAP(N,XX(1,I),1,XX(1,J),1)
  613.             DO 40 K = I1, P
  614.                IF (JPVT(K) .EQ. I) JPVT(K) = JPVT(I)
  615.    40       CONTINUE
  616.    50    CONTINUE
  617.    60 CONTINUE
  618.       RETURN
  619.       END
  620.       SUBROUTINE SXGEN(X,LDX,N,P,S)
  621.       INTEGER LDX,N,P
  622.       REAL X(LDX,1),S(1)
  623. C
  624. C     SXGEN GENERATES AN N X P MATRIX X HAVING SINGULAR
  625. C     VALUES S(1), S(2), ..., S(M), WHERE M = MIN(N,P)
  626. C
  627. C     FORTRAN FLOAT,MIN0
  628. C
  629.       INTEGER I,J,M,MP1
  630.       REAL FN,FP
  631.       REAL FAC,RU,T
  632.       FP = FLOAT(P)
  633.       FN = FLOAT(N)
  634.       M = MIN0(N,P)
  635.       RU = 1.0E0
  636.       FAC = 1.0E0/FP
  637.       DO 20 I = 1, M
  638.          FAC = FAC*RU
  639.          DO 10 J = 1, P
  640.             X(I,J) = -2.0E0*S(I)*FAC
  641.    10    CONTINUE
  642.          X(I,I) = X(I,I) + FP*S(I)*FAC
  643.    20 CONTINUE
  644.       IF (M .GE. N) GO TO 50
  645.          MP1 = M + 1
  646.          DO 40 J = 1, P
  647.             DO 30 I = MP1, N
  648.                X(I,J) = 0.0E0
  649.    30       CONTINUE
  650.    40    CONTINUE
  651.    50 CONTINUE
  652.       DO 80 J = 1, P
  653.          T = 0.0E0
  654.          DO 60 I = 1, N
  655.             T = T + X(I,J)
  656.    60    CONTINUE
  657.          DO 70 I = 1, N
  658.             X(I,J) = X(I,J) - 2.0E0*T/FN
  659.    70    CONTINUE
  660.    80 CONTINUE
  661.       RETURN
  662.       END
  663.  
  664.  
  665.