c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
c RCS version control information:
c $Header: /home/paul/rt/sbdart/RCS/sbdart.f,v 1.4 2000/07/13 23:43:37 paul Exp $
c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

c Call tree:
c
c    SGBCO
c       SASUM
c       SDOT
c       SAXPY
c       SGBFA
c           ISAMAX
c           SAXPY
c           SSCAL
c       SSCAL
c   SGBSL
c       SDOT
c       SAXPY
c   SGECO
c       SASUM
c       SDOT
c       SAXPY
c       SGEFA
c           ISAMAX
c           SAXPY
c           SSCAL
c       SSCAL
c   SGESL
c       SDOT
c       SAXPY
c   SSWAP
c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



      SUBROUTINE SGBCO( ABD, LDA, N, ML, MU, IPVT, RCOND, Z ) 2,17

c         Factors a real band matrix by Gaussian elimination
c         and estimates the condition of the matrix.

c         Revision date:  8/1/82
c         Author:  Moler, C. B. (U. of New Mexico)

c     If  RCOND  is not needed, SGBFA is slightly faster.
c     To solve  A*X = B , follow SBGCO by SGBSL.

c     input:

C        ABD     REAL(LDA, N)
c                contains the matrix in band storage.  The columns
c                of the matrix are stored in the columns of  ABD  and
c                the diagonals of the matrix are stored in rows
c                ML+1 through 2*ML+MU+1 of  ABD .
c                See the comments below for details.

C        LDA     INTEGER
c                the leading dimension of the array  ABD .
c                LDA must be .GE. 2*ML + MU + 1 .

C        N       INTEGER
c                the order of the original matrix.

C        ML      INTEGER
c                number of diagonals below the main diagonal.
c                0 .LE. ML .LT. N .

C        MU      INTEGER
c                number of diagonals above the main diagonal.
c                0 .LE. MU .LT. N .
c                more efficient if  ML .LE. MU .

c     on return

c        ABD     an upper triangular matrix in band storage and
c                the multipliers which were used to obtain it.
c                The factorization can be written  A = L*U  where
c                L  is a product of permutation and unit lower
c                triangular matrices and  U  is upper triangular.

C        IPVT    INTEGER(N)
c                an integer vector of pivot indices.

C        RCOND   REAL
c                an estimate of the reciprocal condition of  A .
c                For the system  A*X = B , relative perturbations
c                in  A  and  B  of size  epsilon  may cause
c                relative perturbations in  X  of size  epsilon/RCOND .
c                If  RCOND  is so small that the logical expression
c                           1.0 + RCOND .EQ. 1.0
c                is true, then  A  may be singular to working
c                precision.  In particular,  RCOND  is zero  if
c                exact singularity is detected or the estimate
c                underflows.

C        Z       REAL(N)
c                a work vector whose contents are usually unimportant.
c                If  A  is close to a singular matrix, then  Z  is
c                an approximate null vector in the sense that
c                norm(a*z) = rcond*norm(a)*norm(z) .

c     Band storage

c           If  A  is a band matrix, the following program segment
c           will set up the input.

c                   ML = (band width below the diagonal)
c                   MU = (band width above the diagonal)
c                   M = ML + MU + 1
c                   DO 20 J = 1, N
c                      I1 = MAX(1, J-MU)
c                      I2 = MIN(N, J+ML)
c                      DO 10 I = I1, I2
c                         K = I - J + M
c                         ABD(K,J) = A(I,J)
c                10    CONTINUE
c                20 CONTINUE

c           This uses rows  ML+1  through  2*ML+MU+1  of  ABD .
c           In addition, the first  ML  rows in  ABD  are used for
c           elements generated during the triangularization.
c           The total number of rows needed in  ABD  is  2*ML+MU+1 .
c           The  ML+MU by ML+MU  upper left triangle and the
c           ML by ML  lower right triangle are not referenced.

c     Example:  if the original matrix is

c           11 12 13  0  0  0
c           21 22 23 24  0  0
c            0 32 33 34 35  0
c            0  0 43 44 45 46
c            0  0  0 54 55 56
c            0  0  0  0 65 66

c      then  N = 6, ML = 1, MU = 2, LDA .GE. 5  and ABD should contain

c            *  *  *  +  +  +  , * = not used
c            *  * 13 24 35 46  , + = used for pivoting
c            * 12 23 34 45 56
c           11 22 33 44 55 66
c           21 32 43 54 65  *

c --------------------------------------------------------------------


c     .. Scalar Arguments ..

      INTEGER   LDA, ML, MU, N
      REAL      RCOND
c     ..
c     .. Array Arguments ..

      INTEGER   IPVT( * )
      REAL      ABD( LDA, * ), Z( * )
c     ..
c     .. Local Scalars ..

      INTEGER   INFO, IS, J, JU, K, KB, KP1, L, LA, LM, LZ, M, MM
      REAL      ANORM, EK, S, SM, T, WK, WKM, YNORM
c     ..
c     .. External Functions ..

      REAL      SASUM, SDOT
      EXTERNAL  SASUM, SDOT
c     ..
c     .. External Subroutines ..

      EXTERNAL  SAXPY, SGBFA, SSCAL
c     ..
c     .. Intrinsic Functions ..

      INTRINSIC ABS, MAX, MIN, SIGN
c     ..


c                       ** compute 1-norm of A
      ANORM  = 0.0E0
      L  = ML + 1
      IS = L + MU

      DO 10 J = 1, N

         ANORM  = MAX( ANORM, SASUM( L,ABD( IS,J ),1 ) )

         IF( IS.GT.ML + 1 ) IS = IS - 1

         IF( J.LE.MU ) L  = L + 1

         IF( J.GE.N - ML ) L  = L - 1

   10 CONTINUE
c                                               ** factor

      CALL SGBFA( ABD, LDA, N, ML, MU, IPVT, INFO )

c     RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))) .
c     estimate = norm(Z)/norm(Y) where  A*Z = Y  and  trans(A)*Y = E.
c     trans(A) is the transpose of A.  The components of E  are
c     chosen to cause maximum local growth in the elements of W  where
c     trans(U)*W = E.  The vectors are frequently rescaled to avoid
c     overflow.

c                     ** solve trans(U)*W = E
      EK = 1.0E0

      DO 20 J = 1, N
         Z( J ) = 0.0E0
   20 CONTINUE


      M  = ML + MU + 1
      JU = 0

      DO 50 K = 1, N

         IF( Z( K ).NE.0.0E0 ) EK = SIGN( EK, -Z( K ) )

         IF( ABS( EK - Z( K ) ).GT.ABS( ABD( M,K ) ) ) THEN

            S  = ABS( ABD( M,K ) ) / ABS( EK - Z( K ) )

            CALL SSCAL( N, S, Z, 1 )

            EK = S*EK

         END IF

         WK   = EK - Z( K )
         WKM  = -EK - Z( K )
         S    = ABS( WK )
         SM   = ABS( WKM )

         IF( ABD( M,K ).NE.0.0E0 ) THEN

            WK   = WK / ABD( M, K )
            WKM  = WKM / ABD( M, K )

         ELSE

            WK   = 1.0E0
            WKM  = 1.0E0

         END IF

         KP1  = K + 1
         JU   = MIN( MAX( JU,MU + IPVT( K ) ), N )
         MM   = M

         IF( KP1.LE.JU ) THEN

            DO 30 J = KP1, JU
               MM     = MM - 1
               SM     = SM + ABS( Z( J ) + WKM*ABD( MM,J ) )
               Z( J ) = Z( J ) + WK*ABD( MM, J )
               S      = S + ABS( Z( J ) )
   30       CONTINUE

            IF( S.LT.SM ) THEN

               T  = WKM - WK
               WK = WKM
               MM = M

               DO 40 J = KP1, JU
                  MM = MM - 1
                  Z( J ) = Z( J ) + T*ABD( MM, J )
   40          CONTINUE

            END IF

         END IF

         Z( K ) = WK

   50 CONTINUE


      S  = 1.0E0 / SASUM( N, Z, 1 )

      CALL SSCAL( N, S, Z, 1 )

c                         ** solve trans(L)*Y = W
      DO 60 KB = 1, N
         K  = N + 1 - KB
         LM = MIN( ML, N - K )

         IF( K.LT.N )
     &       Z( K ) = Z( K ) + SDOT( LM, ABD( M+1, K ), 1, Z( K+1 ), 1 )

         IF( ABS( Z( K ) ).GT.1.0E0 ) THEN

            S  = 1.0E0 / ABS( Z( K ) )

            CALL SSCAL( N, S, Z, 1 )

         END IF

         L      = IPVT( K )
         T      = Z( L )
         Z( L ) = Z( K )
         Z( K ) = T

   60 CONTINUE


      S  = 1.0E0 / SASUM( N, Z, 1 )

      CALL SSCAL( N, S, Z, 1 )

      YNORM  = 1.0E0
c                         ** solve L*V = Y
      DO 70 K = 1, N

         L      = IPVT( K )
         T      = Z( L )
         Z( L ) = Z( K )
         Z( K ) = T
         LM     = MIN( ML, N - K )

         IF( K.LT.N )
     &       CALL SAXPY( LM, T, ABD( M+1, K ), 1, Z( K+1 ), 1 )

         IF( ABS( Z(K) ).GT.1.0E0 ) THEN

            S  = 1.0E0 / ABS( Z(K) )

            CALL SSCAL( N, S, Z, 1 )

            YNORM  = S*YNORM

         END IF

   70 CONTINUE


      S  = 1.0E0 / SASUM( N, Z, 1 )

      CALL SSCAL( N, S, Z, 1 )

      YNORM  = S*YNORM

c                           ** solve  U*Z = W
      DO 80 KB = 1, N

         K  = N + 1 - KB

         IF( ABS( Z( K ) ).GT.ABS( ABD( M,K ) ) ) THEN

            S  = ABS( ABD( M,K ) ) / ABS( Z( K ) )

            CALL SSCAL( N, S, Z, 1 )

            YNORM  = S*YNORM

         END IF

         IF( ABD( M,K ).NE.0.0E0 ) Z( K ) = Z( K ) / ABD( M, K )
         IF( ABD( M,K ).EQ.0.0E0 ) Z( K ) = 1.0E0

         LM = MIN( K, M ) - 1
         LA = M - LM
         LZ = K - LM
         T  = -Z( K )

         CALL SAXPY( LM, T, ABD( LA,K ), 1, Z( LZ ), 1 )

   80 CONTINUE
c                              ** make znorm = 1.0

      S  = 1.0E0 / SASUM( N, Z, 1 )

      CALL SSCAL( N, S, Z, 1 )

      YNORM  = S*YNORM
      IF( ANORM.NE.0.0E0 ) RCOND  = YNORM / ANORM
      IF( ANORM.EQ.0.0E0 ) RCOND  = 0.0E0

      END