C
C     THIS DRIVER TESTS  EISPACK  FOR THE CLASS OF REAL
C     BAND MATRICES EXHIBITING THE USE OF  EISPACK  TO FIND THE
C     SOLUTION TO THE EQUATION  A*X = B .  THIS DRIVER SUMMARIZES
C     THE FIGURES OF MERIT FOR THE PATH.
C
C     THIS DRIVER IS CATALOGUED AS EISPDRV4(RBLSUMAR).
C
C     THE DIMENSION OF  A  AND  AHOLD  SHOULD BE  NM  BY  MBB.
C     THE DIMENSION OF  X,B,  AND  BHOLD  SHOULD BE  NM  BY  NM.
C     THE DIMENSION OF  RESDUL  SHOULD BE  NM.
C     THE DIMENSION OF  RV  SHOULD BE  NM*MBB.
C     THE DIMENSION OF  RV1  AND  RV6  SHOULD BE  NM.
C     HERE  NM = 20  AND  MBB = 39.
C
C
      DOUBLE PRECISION A( 20, 39),X( 20, 20),B( 20, 20),RV( 780),
     X        RESDUL( 20),RV1( 20),RV6( 20),TCRIT( 1),EPSLON,
     X        RESMAX,E21,DFL
      DOUBLE PRECISION AHOLD( 20, 39),BHOLD( 20, 20)
      INTEGER IERR( 1),ERROR
      DATA IREAD1/1/,IREADC/5/,IWRITE/6/
C
      OPEN(UNIT=IREAD1,FILE='FILE49')
      OPEN(UNIT=IREADC,FILE='FILE50')
      REWIND IREAD1
      REWIND IREADC
C
      NM = 20
      MBB = 39
      LCOUNT = 0
      WRITE(IWRITE,1)
    1 FORMAT(1H1,19X,57H EXPLANATION OF COLUMN ENTRIES FOR THE SUMMARY S
     XTATISTICS//1H ,95(1H-)/22H ORDER B.W. IP   BANDV /
     X1H ,95(1H-)//
     X49H UNDER 'ORDER' IS THE ORDER OF EACH TEST MATRIX. //
     X52H UNDER 'B.W.' IS THE FULL BAND WIDTH OF THE MATRIX. //
     X59H UNDER 'IP' IS THE COLUMN DIMENSION OF THE CONSTANT MATRIX. //
     X62H UNDER 'BANDV' ARE TWO NUMBERS.  THE FIRST NUMBER, AN INTEGER,,
     X29H IS THE ABSOLUTE VALUE OF THE  /
     X92H ERROR FLAG RETURNED FROM  BANDV.  THE SECOND NUMBER IS THE MEA
     XSURE OF PERFORMANCE BASED     /
     X46H UPON THE RESIDUAL COMPUTED FROM THE SOLUTION. /)
      WRITE(IWRITE,15)
   15 FORMAT(1X,21HD.P. VERSION 04/15/83 )
    5 FORMAT( 80H1       TABULATION OF THE ERROR FLAG  ERROR  AND THE ME
     XASURE OF PERFORMANCE FOR /5X,13HTHE  EISPACK ,
     X73H CODES.  THIS RUN DISPLAYS THESE STATISTICS FOR REAL BAND LINEA
     XR SYSTEMS. //
     X22H ORDER B.W. IP   BANDV )
   10 CALL  RMATIN(NM,MBB,N,MBW,A,AHOLD,BHOLD,0,0)
      CALL  RMATIN(NM,NM,N,IP,X,AHOLD,BHOLD,1,0)
      READ(IREADC,20) IE21
   20 FORMAT(I6)
      E21 = IE21
C
C     RBL  USING  BANDV
C
      ICT = 1
      DO 905 I = 1,N
         RV1(I) = 0.0D0
  905 CONTINUE
      CALL  BANDV(NM,N,MBW,A,E21,IP,RV1,X,ERROR,780,RV,RV6)
      IERR(ICT) = IABS(ERROR)
      CALL  RMATIN(NM,MBB,N,MBW,A,AHOLD,BHOLD,0,1)
      CALL  RMATIN(NM,NM,N,IP,B,AHOLD,BHOLD,1,1)
      CALL  RBLRES(NM,N,MBW,IP,A,B,X,RESDUL,E21)
      RESMAX = 0.0D0
      DO 935 I = 1,IP
         IF( RESDUL(I) .GT. RESMAX ) RESMAX = RESDUL(I)
  935 CONTINUE
      DFL = N
      TCRIT(ICT) = RESMAX/EPSLON(DFL)
C
      IF( MOD(LCOUNT,35) .EQ. 0 ) WRITE(IWRITE,5)
      LCOUNT = LCOUNT + 1
      MBF = MBW
      IF( E21 .EQ. 1.0D0 ) MBF = MBW*2 - 1
      WRITE(IWRITE,520) N,MBF,IP,IERR(1),TCRIT(1)
  520 FORMAT(1X,I3,2I5,I4,F6.2)
      GO TO 10
      END
      SUBROUTINE RBLRES(NM,N,MBW,IP,A,B,X,NORM,E21)
C
      DOUBLE PRECISION A(NM,MBW),B(NM,IP),X(NM,IP),NORM(N)
      DOUBLE PRECISION NORMA,NORMX,NORMR,SUM,SUMA,E21
C
C     THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX
C     A*X-B  WHERE  A  IS A REAL BAND MATRIX, WHICH MAY BE SYMMETRIC,
C     X  IS A REAL  N  BY  IP  MATRIX, AND   B IS A REAL
C     N  BY  IP  MATRIX.
C
C     THIS SUBROUTINE IS CATALOGUED AS EISPDRV4(RBLRES).
C
C     INPUT.
C        NM IS THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS
C           AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT;
C
C        N IS THE ROW DIMENSION OF THE MATRIX  A;
C
C        MBW IS THE NUMBER OF COLUMNS OF THE ARRAY A USED TO STORE THE
C           BAND MATRIX.  IF THE MATRIX IS SYMMETRIC, MBW IS ITS (HALF)
C           BAND WIDTH, DENOTED MB AND DEFINED AS THE NUMBER OF ADJACENT
C           DIAGONALS, INCLUDING THE PRINCIPAL DIAGONAL, REQUIRED TO
C           SPECIFY THE NON-ZERO PORTION OF THE LOWER TRIANGLE OF THE
C           MATRIX.  IF THE MATRIX IS NON-SYMMETRIC IT MUST HAVE THE SAM
C           NUMBER OF ADJACENT DIAGONALS ABOVE THE MAIN DIAGONAL AS
C           BELOW.  IN THIS CASE  MBW = 2*MB-1;
C
C        IP  IS THE COLUMN DIMENSION OF  B  AND  X;
C
C        A(N,MBW) IS AN ARRAY WHICH CONTAINS IN ITS COLUMNS THE
C           SUBDIAGONALS AND DIAGONAL OF THE SYMMETRIC BAND
C           MATRIX. IF  A  IS NON-SYMMETRIC IT ALSO CONTAINS
C           THE SUPER-DIAGONALS;
C
C        B  IS A REAL FULL  N  BY  IP  MATRIX;
C
C        X  IS A REAL FULL  N  BY  IP  MATRIX;
C
C        E21 TELLS IF THE MATRIX IS SYMMETRIC.  IF E21 EQUALS ONE THE
C           MATRIX IS SYMMETRIC AND IF E21 EQUALS MINUS ONE THE
C           MATRIX IS NON-SYMMETRIC.
C
C     OUTPUT.
C
C        NORM(N)  IS AN ARRAY SUCH THAT FOR EACH K,
C               NORM(K) = !!A*X(K)-B(K)!!/(!!A!!*!!X(K)!!) .
C
C     ------------------------------------------------------------------
C
      IF( E21 .EQ. 1.0D0 ) GO TO 75
      MBWT = MBW
      MB = MBW/2 + 1
      NORMA = 0.0D0
C
      DO 35 I = 1,N
         L = MAX0(MB - I + 1,1)
         M = MIN0(MB - I + N,MBWT)
         SUMA = 0.0D0
C
         DO 34 J = L,M
            SUMA = SUMA + DABS(A(I,J))
   34    CONTINUE
C
         NORMA = DMAX1(NORMA,SUMA)
   35 CONTINUE
C
      DO 70 IPP = 1,IP
         NORMR = 0.0D0
         NORMX = 0.0D0
         DO 50 I = 1,N
            L = MAX0(MB - I + 1,1)
            M = MIN0(MB - I + N,MBWT)
            K = MAX0(-MB + 1 + I,1)
            SUM = -B(I,IPP)
            SUMA = 0.0D0
C
            DO 40 J = L,M
               SUM = SUM + A(I,J)*X(K,IPP)
               K = K + 1
               SUMA = SUMA + DABS(A(I,J))
   40       CONTINUE
C
            NORMR = NORMR + DABS(SUM)
            NORMX = NORMX + DABS(X(I,IPP))
   50    CONTINUE
C
         IF( NORMX .EQ. 0.0D0 ) NORMX = 1.0D0
         IF( NORMA .EQ. 0.0D0 ) NORMA = 1.0D0
         NORM(IPP) = NORMR/(NORMA*NORMX)
   70 CONTINUE
C
      RETURN
C
   75 MB = MBW
      NORMA = 0.0D0
      MB1 = MB - 1
C
      DO 110 I = 1,N
         SUMA = 0.0D0
         IF( I .EQ. 1 ) GO TO 90
         LSTART = MAX0(1,MB + 1 - I)
C
         DO 80 L = LSTART,MB1
            SUMA = SUMA + DABS(A(I,L))
   80    CONTINUE
C
   90    LSTOP = MIN0(MB,N + 1 - I)
         J = I
C
         DO 100 L = 1,LSTOP
            L1 = MB + 1 - L
            SUMA = SUMA + DABS(A(J,L1))
            J = J + 1
  100    CONTINUE
C
         NORMA = DMAX1(NORMA,SUMA)
  110 CONTINUE
C
      DO 170 I = 1,IP
         NORMR = 0.0D0
         NORMX = 0.0D0
C
         DO 150 L = 1,N
            SUM = -B(L,I)
            J = MAX0(0,L - MB)
            IF( L .EQ. 1 ) GO TO 130
            KSTART = MAX0(1,MB + 1 - L)
C
            DO 120 K = KSTART,MB1
               J = J + 1
               SUM = SUM + A(L,K)*X(J,I)
  120       CONTINUE
C
  130       KSTOP = MIN0(MB,N + 1 - L)
C
            DO 140 K = 1,KSTOP
               J = J + 1
               K1 = MB + 1 - K
               SUM = SUM + A(J,K1)*X(J,I)
  140       CONTINUE
C
            NORMR = NORMR + DABS(SUM)
  150    CONTINUE
C
         DO 160 K = 1,N
            NORMX = NORMX + DABS(X(K,I))
  160    CONTINUE
C
         IF( NORMX .EQ. 0.0D0 ) NORMX = 1.0D0
         IF( NORMA .EQ. 0.0D0 ) NORMA = 1.0D0
         NORM(I) = NORMR/(NORMA*NORMX)
  170 CONTINUE
      RETURN
C
      END
      SUBROUTINE RMATIN(NM,NN,M,N,A,AHOLD,BHOLD,AORB,INITIL)
C
C     THIS INPUT SUBROUTINE READS A REAL MATRIX, WITH DIMENSION
C     M  BY  N, FROM SYSIN.
C     TO GENERATE THE MATRIX  A (OR  B)  INITIALLY,  INITIL  IS
C     SET TO 0.  TO REGENERATE THE MATRIX  A (OR  B)  FOR THE
C     PURPOSE OF THE RESIDUAL CALCULATION,  INITIL  IS SET TO 1.
C
C     THIS ROUTINE IS CATALOGUED AS  EISPDRV4(RLREADI).
C
      INTEGER IA( 20),NM,NN,M,N,AORB,INITIL
      DOUBLE PRECISION A(NM,NN),AHOLD(NM,NN),BHOLD(NM,NM)
      DATA IREADA/1/,IWRITE/6/
C
      IF( INITIL .EQ. 1 ) GO TO 90
      READ(IREADA,5) M,N
    5 FORMAT(I6,6X,I6)
      IF( M .EQ. 0 ) GO TO 170
      DO 20 I = 1,M
         READ(IREADA,10) (IA(J),J=1,N)
   10    FORMAT(6I12)
         DO 15 J = 1,N
            A(I,J) = IA(J)
   15    CONTINUE
   20 CONTINUE
      IF( AORB .EQ. 1 ) GO TO 50
      DO 40 J = 1,N
         DO 30 I = 1,M
            AHOLD(I,J) = A(I,J)
   30    CONTINUE
   40 CONTINUE
      GO TO 80
   50 DO 70 J = 1,N
         DO 60 I = 1,M
            BHOLD(I,J) = A(I,J)
   60    CONTINUE
   70 CONTINUE
   80 RETURN
C
   90 IF( AORB .EQ. 1 ) GO TO 120
      DO 110 J = 1,N
         DO 100 I = 1,M
            A(I,J) = AHOLD(I,J)
  100    CONTINUE
  110 CONTINUE
      GO TO 150
  120 DO 140 J = 1,N
         DO 130 I = 1,M
            A(I,J) = BHOLD(I,J)
  130    CONTINUE
  140 CONTINUE
  150 RETURN
C
  170 WRITE(IWRITE,175)
  175 FORMAT(45H0END OF DATA FOR SUBROUTINE RMATIN(RLREADI). /1H1)
      STOP
      END