SUBROUTINE SQRTS(LUNIT)
C     LUNIT IS THE OUTPUT UNIT NUMBER
C
C     TESTS
C        SQRDC,SQRSL
C
C     LINPACK. THIS VERSION DATED 08/14/78 .
C     G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
C
C     SUBROUTINES AND FUNCTIONS
C
C     LINPACK SQRDC,SQRSL
C     EXTERNAL SQRBM,SQRFM,SQRRX,SXGEN,SMACH
C     FORTRAN AMAX1,ABS,FLOAT
C
C     INTERNAL VARIABLES
C
      INTEGER LUNIT
      INTEGER CASE,LDX,N,P,PD2,PD2P1,NCASE,JPVT(25)
      INTEGER I,INFO,J,JJ,JJJ,L
      REAL BESTAT,RMAX,RSTAT,XBMAX,XBSTAT,ERRLVL,BSTAT,FSTAT
      REAL BETA(25),QRAUX(25),QY(25),RSD(25),S(25),WORK(25)
      REAL X(25,25),XBETA(25),XX(25,25),Y(25),Y1(25),Y2(25)
      REAL SMACH,TINY,HUGE
      LOGICAL NOTWRT
      LDX = 25
      NCASE = 9
      ERRLVL = 100.0E0
      NOTWRT = .TRUE.
      TINY = SMACH(2)
      HUGE = SMACH(3)
      WRITE (LUNIT,600)
      DO 550 CASE = NCASE, NCASE
         WRITE (LUNIT,560) CASE
         XBSTAT = 0.0E0
         BESTAT = 0.0E0
         GO TO (10, 100, 170, 240, 300, 330, 380, 430, 470), CASE
C
C        WELL CONDITIONED LEAST SQUARES PROBLEM.
C
   10    CONTINUE
            WRITE (LUNIT,640)
            N = 10
            P = 4
C
C           GENERATE A MATRIX X WITH SINGULAR VALUES 1. AND .5.
C
            PD2 = P/2
            PD2P1 = PD2 + 1
            DO 20 L = 1, PD2
               S(L) = 1.0E0
   20       CONTINUE
            DO 30 L = PD2P1, P
               S(L) = 0.5E0
   30       CONTINUE
            CALL SXGEN(X,LDX,N,P,S)
            IF (NOTWRT) GO TO 40
               WRITE (LUNIT,610)
               CALL SARRAY(X,LDX,N,P,4,LUNIT)
   40       CONTINUE
C           THE OBSERVATION VECTOR IS Y = Y1 + Y2, WHERE Y1 IS
C           THE VECTOR OF ROW SUMS OF X AND Y2 IS ORTHOGONAL TO
C           THE COLUMNS OF X.
C
            DO 60 I = 1, N
               Y1(I) = 0.0E0
               Y2(I) = 2.0E0*TINY
               IF (I .EQ. P + 1) Y2(I) = Y2(I) - FLOAT(N)*TINY
               DO 50 J = 1, P
                  X(I,J) = X(I,J)*TINY
                  Y1(I) = Y1(I) + X(I,J)
                  XX(I,J) = X(I,J)
   50          CONTINUE
               Y(I) = Y1(I) + Y2(I)
   60       CONTINUE
C
C           REDUCE THE MATRIX.
C
            CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,0)
            IF (NOTWRT) GO TO 70
               WRITE (LUNIT,610)
               CALL SARRAY(X,LDX,N,P,4,LUNIT)
               WRITE (LUNIT,620)
               CALL SARRAY(QRAUX,P,P,1,-4,LUNIT)
   70       CONTINUE
C
C           COMPUTE THE STATISTICS FOR THE FOWARD AND BACK
C           MULTIPLICATIONS.
C
            CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
            CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
C
C           SOLVE THE LEAST SQUARES PROBLEM.
C
            CALL SQRSL(X,LDX,N,P,QRAUX,Y,QY,RSD,BETA,RSD,XBETA,111,
     *                 INFO)
C
C           COMPUTE THE LEAST SQUARES TEST STATISTICS.
C
            RSTAT = 0.0E0
            RMAX = 0.0E0
            XBSTAT = 0.0E0
            XBMAX = 0.0E0
            DO 80 I = 1, N
               RSTAT = AMAX1(RSTAT,ABS(Y2(I)-RSD(I)))
               RMAX = AMAX1(RMAX,ABS(Y2(I)))
               XBSTAT = AMAX1(XBSTAT,ABS(Y1(I)-XBETA(I)))
               XBMAX = AMAX1(XBMAX,ABS(Y2(I)))
   80       CONTINUE
            RSTAT = RSTAT/(SMACH(1)*RMAX)
            XBSTAT = XBSTAT/(SMACH(1)*XBMAX)
            BESTAT = 0.0E0
            DO 90 J = 1, P
               BESTAT = AMAX1(BESTAT,ABS(1.0E0-BETA(J)))
   90       CONTINUE
            BESTAT = BESTAT/SMACH(1)
            WRITE (LUNIT,570) FSTAT,BSTAT,BESTAT,XBSTAT,RSTAT
         GO TO 540
C
C        4 X 10 MATRIX.
C
  100    CONTINUE
            WRITE (LUNIT,650)
            N = 4
            P = 10
            PD2 = P/2
            PD2P1 = PD2 + 1
            DO 110 L = 1, PD2
               S(L) = 1.0E0
  110       CONTINUE
            DO 120 L = PD2P1, P
               S(L) = 0.5E0
  120       CONTINUE
            CALL SXGEN(X,LDX,N,P,S)
            IF (NOTWRT) GO TO 130
               WRITE (LUNIT,610)
               CALL SARRAY(X,LDX,N,P,10,LUNIT)
  130       CONTINUE
            DO 150 I = 1, N
               DO 140 J = 1, P
                  XX(I,J) = X(I,J)
  140          CONTINUE
  150       CONTINUE
            CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,0)
            IF (NOTWRT) GO TO 160
               WRITE (LUNIT,610)
               CALL SARRAY(X,LDX,N,P,10,LUNIT)
               WRITE (LUNIT,620)
               CALL SARRAY(QRAUX,N,N,1,-4,LUNIT)
  160       CONTINUE
            CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
            CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
            WRITE (LUNIT,570) FSTAT,BSTAT
         GO TO 540
C
C        PIVOTING AND OVERFLOW TEST.  COLUMNS 1, 4, 7 FROZEN.
C
  170    CONTINUE
            WRITE (LUNIT,670)
            N = 10
            P = 3
            S(1) = 1.0E0
            S(2) = 0.5E0
            S(3) = 0.25E0
            CALL SXGEN(X,LDX,N,P,S)
            P = 9
            DO 190 I = 1, N
               DO 180 JJJ = 1, 3
                  J = 4 - JJJ
                  JJ = 3*J
                  X(I,JJ) = HUGE*X(I,J)
                  X(I,JJ-1) = X(I,JJ)/2.0E0
                  X(I,JJ-2) = X(I,JJ-1)/2.0E0
  180          CONTINUE
  190       CONTINUE
            IF (NOTWRT) GO TO 200
               WRITE (LUNIT,610)
               CALL SARRAY(X,LDX,N,P,9,LUNIT)
  200       CONTINUE
            DO 220 J = 1, P
               JPVT(J) = 0
               DO 210 I = 1, N
                  XX(I,J) = X(I,J)
  210          CONTINUE
  220       CONTINUE
            JPVT(1) = -1
            JPVT(4) = -1
            JPVT(7) = -1
            CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
            IF (NOTWRT) GO TO 230
               WRITE (LUNIT,610)
               CALL SARRAY(X,LDX,N,P,9,LUNIT)
               WRITE (LUNIT,620)
               CALL SARRAY(QRAUX,P,P,1,-9,LUNIT)
  230       CONTINUE
            WRITE (LUNIT,630) (JPVT(J), J = 1, P)
            CALL SQRRX(XX,LDX,N,P,JPVT)
            CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
            CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
            WRITE (LUNIT,570) FSTAT,BSTAT
         GO TO 540
C
C        25 BY 25 MATRIX
C
  240    CONTINUE
            WRITE (LUNIT,680)
            N = 25
            P = 25
            DO 250 I = 1, 25
               S(I) = 2.0E0**(1 - I)
  250       CONTINUE
            CALL SXGEN(X,LDX,N,P,S)
            IF (NOTWRT) GO TO 260
               WRITE (LUNIT,610)
               CALL SARRAY(X,LDX,N,P,10,LUNIT)
  260       CONTINUE
            DO 280 J = 1, P
               JPVT(J) = 0
               DO 270 I = 1, N
                  XX(I,J) = X(I,J)
  270          CONTINUE
  280       CONTINUE
            CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
            IF (NOTWRT) GO TO 290
               WRITE (LUNIT,610)
               CALL SARRAY(X,LDX,N,P,10,LUNIT)
               WRITE (LUNIT,620)
               CALL SARRAY(QRAUX,P,P,1,-10,LUNIT)
  290       CONTINUE
            WRITE (LUNIT,630) (JPVT(J), J = 1, P)
            CALL SQRRX(XX,LDX,N,P,JPVT)
            CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
            CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
            WRITE (LUNIT,570) FSTAT,BSTAT
         GO TO 540
C
C        MONOELEMENTAL MATRIX.
C
  300    CONTINUE
            WRITE (LUNIT,690)
            N = 1
            P = 1
            X(1,1) = 1.0E0
            IF (NOTWRT) GO TO 310
               WRITE (LUNIT,610)
               CALL SARRAY(X,LDX,N,P,1,LUNIT)
  310       CONTINUE
            XX(1,1) = 1.0E0
            JPVT(1) = 0
            CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
            IF (NOTWRT) GO TO 320
               WRITE (LUNIT,610)
               CALL SARRAY(X,LDX,N,P,1,LUNIT)
               WRITE (LUNIT,620)
               CALL SARRAY(QRAUX,P,P,1,-1,LUNIT)
  320       CONTINUE
            CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
            CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
            WRITE (LUNIT,570) FSTAT,BSTAT
         GO TO 540
C
C        ZERO MATRIX.
C
  330    CONTINUE
            WRITE (LUNIT,700)
            N = 10
            P = 4
            DO 350 J = 1, P
               JPVT(J) = 0
               DO 340 I = 1, N
                  X(I,J) = 0.0E0
                  XX(I,J) = 0.0E0
  340          CONTINUE
  350       CONTINUE
            IF (NOTWRT) GO TO 360
               WRITE (LUNIT,610)
               CALL SARRAY(X,LDX,N,P,4,LUNIT)
  360       CONTINUE
            CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
            IF (NOTWRT) GO TO 370
               WRITE (LUNIT,610)
               CALL SARRAY(X,LDX,N,P,4,LUNIT)
               WRITE (LUNIT,620)
               CALL SARRAY(QRAUX,P,P,1,-4,LUNIT)
  370       CONTINUE
            WRITE (LUNIT,630) (JPVT(J), J = 1, P)
            CALL SQRRX(XX,LDX,N,P,JPVT)
            CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
            CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
            WRITE (LUNIT,570) FSTAT,BSTAT
         GO TO 540
C
C        10 X 1 MATRIX WITH LEAST SQUARES PROBLEM.
C
  380    CONTINUE
            WRITE (LUNIT,710)
            N = 10
            P = 1
            S(1) = 1.0E0
            CALL SXGEN(X,LDX,N,P,S)
            IF (NOTWRT) GO TO 390
               WRITE (LUNIT,610)
               CALL SARRAY(X,LDX,N,P,1,LUNIT)
  390       CONTINUE
            DO 400 I = 1, N
               Y1(I) = 0.0E0
               Y2(I) = 2.0E0
               IF (I .EQ. P + 1) Y2(I) = Y2(I) - FLOAT(N)
               Y1(I) = Y1(I) + X(I,1)
               Y(I) = Y1(I) + Y2(I)
               XX(I,1) = X(I,1)
  400       CONTINUE
            CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,0)
            IF (NOTWRT) GO TO 410
               WRITE (LUNIT,610)
               CALL SARRAY(X,LDX,N,P,1,LUNIT)
               WRITE (LUNIT,620)
               CALL SARRAY(QRAUX,P,P,1,1,LUNIT)
  410       CONTINUE
            CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
            CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
            CALL SQRSL(X,LDX,N,P,QRAUX,Y,QY,RSD,BETA,RSD,XBETA,111,
     *                 INFO)
            RSTAT = 0.0E0
            RMAX = 0.0E0
            XBSTAT = 0.0E0
            XBMAX = 0.0E0
            DO 420 I = 1, N
               RSTAT = AMAX1(RSTAT,ABS(Y2(I)-RSD(I)))
               RMAX = AMAX1(RMAX,ABS(Y2(I)))
               XBMAX = AMAX1(XBMAX,ABS(Y2(I)))
               XBSTAT = AMAX1(XBSTAT,ABS(Y1(I)-XBETA(I)))
  420       CONTINUE
            RSTAT = RSTAT/(SMACH(1)*RMAX)
            XBSTAT = XBSTAT/(SMACH(1)*XBMAX)
            BESTAT = ABS(1.0E0-BETA(1))
            BESTAT = BESTAT/SMACH(1)
            WRITE (LUNIT,570) FSTAT,BSTAT,BESTAT,XBSTAT,RSTAT
         GO TO 540
C
C        1 X 4 MATRIX
C
  430    CONTINUE
            WRITE (LUNIT,720)
            N = 1
            P = 4
            X(1,1) = 1.0E0
            X(1,2) = 2.0E0
            X(1,3) = 0.25E0
            X(1,4) = 0.5E0
            DO 440 I = 1, 4
               JPVT(I) = 0
               XX(1,I) = X(1,I)
  440       CONTINUE
            IF (NOTWRT) GO TO 450
               WRITE (LUNIT,610)
               CALL SARRAY(X,LDX,N,P,4,LUNIT)
  450       CONTINUE
            CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
            IF (NOTWRT) GO TO 460
               WRITE (LUNIT,610)
               CALL SARRAY(X,LDX,N,P,10,LUNIT)
               WRITE (LUNIT,620)
               CALL SARRAY(QRAUX,N,N,1,-1,LUNIT)
  460       CONTINUE
            WRITE (LUNIT,630) (JPVT(J), J = 1, P)
            CALL SQRRX(XX,LDX,N,P,JPVT)
            CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
            CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
            WRITE (LUNIT,570) FSTAT,BSTAT
         GO TO 540
C
C        PIVOTING TEST.
C
  470    CONTINUE
            WRITE (LUNIT,660)
            N = 10
            P = 3
            S(1) = 1.0E0
            S(2) = 0.5E0
            S(3) = 0.25E0
            CALL SXGEN(X,LDX,N,P,S)
            P = 9
            DO 490 I = 1, N
               DO 480 JJJ = 1, 3
                  J = 4 - JJJ
                  JJ = 3*J
                  X(I,JJ) = X(I,J)
                  X(I,JJ-1) = X(I,JJ)/2.0E0
                  X(I,JJ-2) = X(I,JJ-1)/2.0E0
  480          CONTINUE
  490       CONTINUE
            IF (NOTWRT) GO TO 500
               WRITE (LUNIT,610)
               CALL SARRAY(X,LDX,N,P,9,LUNIT)
  500       CONTINUE
            DO 520 J = 1, P
               JPVT(J) = 0
               DO 510 I = 1, N
                  XX(I,J) = X(I,J)
  510          CONTINUE
  520       CONTINUE
            CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
            IF (NOTWRT) GO TO 530
               WRITE (LUNIT,610)
               CALL SARRAY(X,LDX,N,P,9,LUNIT)
               WRITE (LUNIT,620)
               CALL SARRAY(QRAUX,P,P,1,-9,LUNIT)
  530       CONTINUE
            WRITE (LUNIT,630) (JPVT(J), J = 1, P)
            CALL SQRRX(XX,LDX,N,P,JPVT)
            CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
            CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
            WRITE (LUNIT,570) FSTAT,BSTAT
  540    CONTINUE
         IF (AMAX1(FSTAT,BSTAT,BESTAT,XBSTAT,RSTAT) .GE. ERRLVL)
     *      WRITE (LUNIT,580)
  550 CONTINUE
      WRITE (LUNIT,590)
      RETURN
C
  560 FORMAT ( // '1CASE', I3 )
  570 FORMAT ( / 11H STATISTICS //
     *         35H    FORWARD MULTIPLICATION ........, E10.2 /
     *         35H    BACK MULTIPLICATION ..........., E10.2 /
     *         35H    BETA .........................., E10.2 /
     *         35H    X*BETA ........................, E10.2 /
     *         35H    RESIDUAL ......................, E10.2)
  580 FORMAT ( / 34H *****STATISTICS ABOVE ERROR LEVEL)
  590 FORMAT ( / 15H1END OF QR TEST)
  600 FORMAT ('1LINPACK TESTER, SQR**' /
     *        ' THIS VERSION DATED 08/14/78.')
  610 FORMAT ( / 2H X)
  620 FORMAT ( / 6H QRAUX)
  630 FORMAT ( / 5H JPVT // (1H , 10I5))
  640 FORMAT ( / 39H WELL CONDITIONED LEAST SQUARES PROBLEM /
     *         20H AND UNDERFLOW TEST.)
  650 FORMAT ( / 14H 4 X 10 MATRIX)
  660 FORMAT ( / 14H PIVOTING TEST /
     *         42H ON RETURN THE FIRST THREE ENTRIES OF JPVT /
     *         36H SHOULD BE 3,6,9 BUT NOT NECESSARILY /
     *         15H IN THAT ORDER.)
  670 FORMAT ( / 27H PIVOTING AND OVERFLOW TEST /
     *         26H WITH COLUMNS 1,4,7 FROZEN /
     *         42H ON RETURN THE LAST  THREE ENTRIES OF JPVT /
     *         31H SHOULD BE 1,4,7 IN THAT ORDER.)
  680 FORMAT ( / 15H 25 X 25 MATRIX)
  690 FORMAT ( / 21H MONOELEMENTAL MATRIX)
  700 FORMAT ( / 12H ZERO MATRIX)
  710 FORMAT ( / 41H 10 X 1 MATRIX WITH LEAST SQUARES PROBLEM)
  720 FORMAT ( / 13H 1 X 4 MATRIX)
      END