C
C     THIS DRIVER TESTS  EISPACK  FOR THE CLASS OF REAL MATRICES
C     EXHIBITING THE USE OF  EISPACK  TO FIND THE SINGULAR VALUES
C     AND THE SOLUTION TO THE EQUATION  A*X = B .  THIS DRIVER
C     SUMMARIZES THE FIGURES OF MERIT FOR ALL PATHS.
C
C     THIS DRIVER IS CATALOGUED AS EISPDRV4(RLSUMARY).
C
C     THE DIMENSION OF  A,U,  AND  V  SHOULD BE  NM  BY  NM.
C     THE DIMENSION OF  SIGMA  AND  RV1  SHOULD BE  NM.
C     THE DIMENSION OF  AHOLD  SHOULD BE  NM  BY  NM.
C     HERE  NM = 20.
C
      REAL A( 20, 20),U( 20, 20),V( 20, 20),
     X        SIGMA( 20),RV1( 20),AHOLD( 20, 20),TCRIT( 2),
     X        RESDUL,EPSLON,DFL
      INTEGER IERR( 2),ERROR
      DATA IREAD1/1/,IWRITE/6/
C
      OPEN(IREAD1,FILE='FILE48')
      REWIND IREAD1
C
      NM = 20
      LCOUNT = 0
      WRITE(*,1)
    1 FORMAT(1H1,19X,57H EXPLANATION OF COLUMN ENTRIES FOR THE SUMMARY S
     XTATISTICS//1H ,95(1H-)/  31H M   N  IP   SVD         MINFIT /
     X1H ,95(1H-)//
     X53H UNDER 'M' IS THE ROW DIMENSION OF EACH TEST MATRIX. //
     X49H UNDER 'N' IS THE COLUMN DIMENSION OF THE MATRIX. //
     X93H UNDER 'IP' IS THE COLUMN DIMENSION OF THE CONSTANT MATRIX.  (A
     XS PRESENTLY CONSTITUTED,      /
     X58H THE IDENTITY MATRIX IS GENERATED AS THE CONSTANT MATRIX.) //
     X93H UNDER 'SVD' ARE TWO NUMBERS.  THE FIRST NUMBER, AN INTEGER, IS
     X THE ERROR FLAG RETURNED      /
     X93H FROM  SVD.  THE SECOND NUMBER IS THE MEASURE OF PERFORMANCE BA
     XSED UPON A RESIDUAL.          //
     X93H UNDER 'MINFIT' ARE TWO NUMBERS.  THE FIRST NUMBER, AN INTEGER,
     X IS THE ERROR FLAG RETURNED   /
     X93H FROM  MINFIT.  THE SECOND NUMBER IS THE MEASURE OF PERFORMANCE
     X BASED UPON A RESIDUAL.       //
     X62H -1.0  AS THE MEASURE OF PERFORMANCE IS PRINTED IF AN ERROR IN,
     X27H THE CORRESPONDING PATH HAS        /
     X50H PREVENTED THE COMPUTATION OF THE SINGULAR VALUES.  /)
      WRITE(*,15)
   15 FORMAT(1X,21HS.P. VERSION 04/15/83 )
    5 FORMAT( 53H1       TABULATION OF THE ERROR FLAG  ERROR  AND THE ,
     X26HMEASURE OF PERFORMANCE FOR /5X, 21HTHE  EISPACK  CODES. ,
     X60H THIS RUN DISPLAYS THESE STATISTICS FOR REAL LINEAR SYSTEMS.//
     X33H   M   N  IP   SVD         MINFIT  )
   10 CALL  RMATIN(NM,NM,M,N,A,AHOLD,AHOLD,0,0)
C
      DO 230 ICALL = 1,2
         IF( ICALL .EQ. 1 ) GO TO 90
         CALL  RMATIN(NM,NM,M,N,A,AHOLD,AHOLD,0,1)
         GO TO 100
C
C     RL  USING  SVD
C
   90    ICT = 1
         CALL  SVD(NM,M,N,A,SIGMA,.TRUE.,U,.TRUE.,V,ERROR,RV1)
         IERR(ICT) = ERROR
         IF( IERR(ICT) .NE. 0 ) GO TO 225
         GO TO 226
C
C     RL  USING  MINFIT
C
  100    ICT = 2
         IP = M
         DO 102 I = 1,M
            DO 101 J = 1,IP
               V(I,J) = 0.0E0
  101       CONTINUE
            V(I,I) = 1.0E0
  102    CONTINUE
         CALL  MINFIT(NM,M,N,A,SIGMA,IP,V,ERROR,RV1)
         IERR(ICT) = ERROR
         IF( IERR(ICT) .NE. 0 ) GO TO 225
         DO 104 I = 1,M
            DO 103 J = 1,N
               U(I,J) = V(J,I)
  103       CONTINUE
  104    CONTINUE
         DO 106 I = 1,N
            DO 105 J = 1,N
               V(I,J) = A(I,J)
  105       CONTINUE
  106    CONTINUE
         GO TO 226
  225    TCRIT(ICT) = -1.0E0
         GO TO 230
C
  226    CALL  RMATIN(NM,NM,M,N,A,AHOLD,AHOLD,0,1)
         CALL  RLRES(NM,M,N,A,U,SIGMA,V,RV1,RESDUL)
         DFL = 10 * MAX0(M,N)
         TCRIT(ICT) = RESDUL/EPSLON(DFL)
  230 CONTINUE
C
      IF( MOD(LCOUNT,35) .EQ. 0 ) WRITE(*,5)
      LCOUNT = LCOUNT + 1
      WRITE(*,520) M,N,IP,IERR(1),TCRIT(1),IERR(2),TCRIT(2)
  520 FORMAT(3I4,I4,F8.3,I4,F8.3)
      GO TO 10
      END
      SUBROUTINE RLRES(NM,M,N,A,U,S,V,NORM,RESDUL)
C
      REAL A(NM,N),U(NM,N),V(NM,N),S(N),NORM(N)
      REAL SUM1,SUM2,X,RESDUL,NORMA,SUMA,SUMUV
C
C     THIS SUBROUTINE PERFORMS A RESIDUAL TEST FOR THE DECOMPOSITION
C                              T
C     OF A MATRIX  A  INTO  USV .
C
C     THIS SUBROUTINE IS CATALOGUED AS EISPDRV4(RLRES).
C
C     INPUT.
C
C       NM  IS THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS
C         AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT;
C
C       M  IS THE ROW DIMENSION OF THE MATRICES  A  AND  U.
C
C       N  IS THE COLUMN DIMENSION OF THE MATRICES  A  AND U , AND
C         THE ORDER OF  V;
C
C       A  IS A REAL GENERAL MATRIX OF DIMEMSION  M  BY  N;
C
C       U  IS A MATRIX WITH ORTHOGONAL COLUMNS OF DIMENSION  M  BY  N;
C
C       S  IS A DIAGONAL MATRIX STORED AS A VECTOR OF DIMENSION  N;
C
C       V  IS AN ORTHOGONAL MATRIX OF DIMENSION  N  BY  N.
C
C     OUTPUT.
C
C       NORM  IS AN ARRAY SUCH THAT FOR EACH K,
C                                              T
C         NORM(K) = !!A*V(K)-S(K)*U(K)!!  + !!A *U(K)-S(K)*V(K)!!
C                   ---------------------------------------------
C                           !!A!!*(!!U(K)!! + !!V(K)!!)
C
C       RESDUL  IS THE LARGEST ELEMENT OF  NORM.
C
C     ------------------------------------------------------------------
C
      RESDUL = 0.0E0
      NORMA = 0.0E0
C
      DO 20 I = 1,M
         SUMA = 0.0E0
C
         DO 10 J = 1,N
            SUMA = SUMA + ABS(A(I,J))
   10    CONTINUE
C
         NORMA = AMAX1(NORMA,SUMA)
   20 CONTINUE
C
      IF( NORMA .EQ. 0.0E0 ) NORMA = 1.0E0
C
      DO 70 K = 1,N
         SUMUV = 0.0E0
         SUM1 = 0.0E0
C
         DO 40 I = 1,M
            SUMUV = SUMUV + ABS(U(I,K))
            X = -S(K)*U(I,K)
C
            DO 30 J = 1,N
               X = X + A(I,J)*V(J,K)
   30       CONTINUE
C
            SUM1 = SUM1 + ABS(X)
   40    CONTINUE
C
         SUM2 = 0.0E0
C
         DO 60 I = 1,N
            SUMUV = SUMUV + ABS(V(I,K))
            X = -S(K)*V(I,K)
C
            DO 50 J = 1,M
               X = X + A(J,I)*U(J,K)
   50       CONTINUE
C
            SUM2 = SUM2 + ABS(X)
   60    CONTINUE
C
         NORM(K) = (SUM1 + SUM2)/(SUMUV*NORMA)
         RESDUL = AMAX1(NORM(K),RESDUL)
   70 CONTINUE
C
      RETURN
      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
      REAL 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(*,175)
  175 FORMAT(45H0END OF DATA FOR SUBROUTINE RMATIN(RLREADI). /1H1)
      STOP
      END
