C$TEST LGEB
C TO RUN AS A MAIN PROGRAM REMOVE NEXT LINE
      SUBROUTINE LGEB
C***********************************************************************
C
C  EXAMPLE OF USE OF THE PORT PROGRAM GECE
C
C***********************************************************************
       INTEGER N, IA, IB, NB, INTER(5), IREAD, I1MACH
       INTEGER I, J, IWRITE, ITER, IEND
       REAL A(5, 5), SAVEA(5, 5), B(5), SAVEB(5), R(5)
       REAL COND, BNORM, R1MACH, ABS, RNORM
       DOUBLE PRECISION DSDOT
C
       N=5
       IA=5
       IB=5
       NB=1
       IREAD=I1MACH(1)
C
       DO 10 I=1,N
   10     READ(IREAD,11) (A(I,J),J=1,N)
   11     FORMAT(1X,5F8.0)
       DO 20 I=1,IB
   20     READ(IREAD,21) B(I)
   21     FORMAT(F8.0)
C
C SAVE THE MATRIX AND RIGHT-HAND SIDE (WHICH WILL BE OVERWRITTEN)
C
       DO 40 I=1,N
         SAVEB(I)=B(I)
         DO 30 J=1,N
   30       SAVEA(I,J)=A(I,J)
   40  CONTINUE
C
C SOLVE AX = B USING SEPARATE CALLS TO GECE, GEFS, GEBS
C
       CALL GECE(N,A,IA,INTER,COND)
       IWRITE=I1MACH(2)
       IF (COND.GE.1.0/R1MACH(4)) WRITE(IWRITE,41)
   41  FORMAT(49H CONDITION NUMBER HIGH,ACCURATE SOLUTION UNLIKELY)
C
       CALL GEFS(N,A,IA,B,IB,NB,INTER)
C
       CALL GEBS(N,A,IA,B,IB,NB)
       WRITE(IWRITE,42)
   42  FORMAT(44H ESTIMATED CONDITION NUMBER OF THE MATRIX A,)
       WRITE(IWRITE,43) COND
   43  FORMAT(27H  USING ONE CALL TO GECE = ,E15.7)
       BNORM=0.0
       WRITE(IWRITE,44)
   44  FORMAT(/22H THE FIRST SOLUTION X,)
       WRITE(IWRITE,45)
   45  FORMAT(41H (USING CALLS TO GECE, GEFS, AND GEBS) = )
C
C COMPUTE NORM OF SOLUTION
C
       DO 50 I=1,N
          BNORM=BNORM + ABS(B(I))
   50     WRITE(IWRITE,51) B(I)
   51     FORMAT(1X, 5F20.7)
C
C REFINE THE SOLUTION DEPENDING ON THE LENGTH OF THE MANTISSA
C
              IEND=I1MACH(11)*IFIX(R1MACH(5)/ALOG10(2.0) + 1.0)
       DO 90 ITER=1,IEND
C  COMPUTE RESIDUAL R = B - AX, IN DOUBLE PRECISION
C
          WRITE(IWRITE,52)
   52     FORMAT(/27H THE RESIDUAL R = B - AX = )
          DO 70 I=1,IA
                DSDOT=0.0
                DO 60 J=1,N
   60              DSDOT = DSDOT  +  DBLE(SAVEA(I,J))*B(J)
                R(I) = SAVEB(I) - DSDOT
   70           WRITE(IWRITE,51) R(I)
C
C  SOLVE LU*(DELTA X) = R USING SEPARATE CALLS TO GEFS AND GEBS
C
          CALL GEFS(N,A,IA,R,IB,NB,INTER)
          CALL GEBS(N,A,IA,R,IB,NB)
C
C  THE NEW SOLUTION X = X  +  DELTA X
C
          WRITE(IWRITE,71)
   71     FORMAT(/38H THE NEW SOLUTION X = X  +  DELTA X = )
C
C DETERMINE NORM OF CORRECTION AND ADD IN CORRECTION
C
          RNORM=0.0
          DO 80 I=1,N
                B(I) = B(I)  +  R(I)
                RNORM=RNORM + ABS(R(I))
   80           WRITE(IWRITE,51) B(I)
C
C TEST FOR CONVERGENCE
C
          IF(RNORM.LT.R1MACH(4)*BNORM) GO TO 100
   90  CONTINUE
       WRITE(IWRITE,91)
   91  FORMAT(/29H ITERATIVE IMPROVEMENT FAILED)
  100  CONTINUE
       STOP
       END
C
C DATA FOR THE EXAMPLE IN THE PORT SHEET...  (REMOVE THE C
C IN COLUMN 1 BEFORE FEEDING THIS DATA TO THE PROGRAM ABOVE.)
C$DATA
C     1.     -2.      3.      7.     -9.
C    -2.      8.     -6.      2.     50.
C     3.     -6.     18.    -15.    -18.
C     7.      2.    -15.    273.    174.
C    -9.     50.    -18.    173.   1667.
C    78.
C  -320.
C   -81.
C   215.
C-10856.
