c
c*******************************************************************************
c
       subroutine GSINS ( Q, ldQ, m, k, aq, r, rr, info )
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      integer  ldQ, m, k, info
c
      real  Q(ldQ,k+1), aq(m), r(k+1), rr(k+1)
c
      external  SCOPY, SNRM2, SSCAL, SDOT, SAXPY
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c  A subroutine for updating a Gram-Schmidt QR factorization
c  
c  Adapted from a subroutine written by G.W. Stewart, September 1981.
c 
c   Given  A = Q R    with  A: m x k , 
c                           Q: m x k , with orthogonal columns, and
c                           R: k x k , upper triangular,
c     and  a = an m-vector,
c
c   Want                        ( R |   )
c         ( A | a ) = ( Q | q ) (---| r )
c                               ( 0 |   )
c
c                     with  ( Q | q ): m x (k+1) , with orthogonal columns, and
c                               r    : (k+1)-vector
c
c-------------------------------------------------------------------------------
c
c  INPUT
c  -----
c 
c  Q       Real matrix ( m x k )
c          Old orthogonal matrix
c
c  ldQ     Row dimension of Fortran array Q in the calling program
c
c  m       Integer
c          No. of rows in matrix Q (and vector aq)
c
c  k       Integer
c          No. of columns in matrix Q
c
c  rr      Real vector ( length = k + 1 )
c          Working vector
c
c-------------------------------------------------------------------------------
c
c  INPUT and OUTPUT
c  ----------------
c
c  aq      Real vector ( length = m )
c          On input: Vector (a) to be orthogonalized against Q
c          On output: Portion (q) of original (a) orthogonal to columns of Q 
c                     (normalized to norm 1) 
c
c-------------------------------------------------------------------------------
c
c  OUTPUT
c  ------
c
c  r       Real vector ( length = k + 1 )
c          Coefficients for expanding (a) as a linear combination of the 
c          the columns of ( Q | q )
c
c  info    Integer
c          = -1, if reorthogonalization failed to converge 
c          =  0, if no trouble
c          =  1, if norm of input aq is zero
c
c-------------------------------------------------------------------------------
c
c  EXAMPLE Of USE
c  --------------
c
c  To update matrices Q and R when a new column A(:,jc) is added to the basis
c
c     call SCOPY ( m, A(1,jc), 1, Q(1,k+1), 1 )
c
c     call GSINS ( Q, ldQ, m, k, Q(1,k+1), R(1,k+1), info )
c
c     k = k + 1
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      integer  i, j, ntry, kp1
c
      real  oldnrm
c
      data ntry /5/
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
       kp1 = k + 1
c
      call SCOPY ( kp1, 0.0, 0, r, 1 )
c
       do 60 i = 1, ntry
           r(kp1) = SNRM2 ( m, aq, 1 )
           if ( r(kp1) .gt. 0 )  go to 20 
           info = 1 
           go to 70 
   20  continue
c
       if ( i .eq. 1 )  go to 40
       if ( r(kp1) .lt. 0.5 * oldnrm )  go to 30
       call SSCAL ( m, 1.0/r(kp1), aq, 1 )
       info = 0
       go to 70
c
   30  continue
c
   40  continue
       oldnrm = r(kp1)
c
       do 50 j = 1, k
           rr(j) = SDOT ( m, Q(1,j), 1, aq, 1 )
           call SAXPY ( m, -rr(j), Q(1,j), 1, aq, 1 )
   50  continue
c
      do 55 j=1,k
        r(j) = r(j) + rr(j)
   55 continue
c
   60  continue
c
       info = -1
c
   70  continue
c
c     --------
       return
c     --------
c
       end
c
c
c*******************************************************************************
c
      subroutine GSERR ( info, jA, jQ, ipr, mou )
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      logical ipr
c
      integer  info, jA, jQ, mou
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c  This subroutine prints error messages from GSINS.  It is called by LS if
c  subroutine GSINS fails when trying to add column jA of the matrix A to the
c  basis set, 
c
c-------------------------------------------------------------------------------
c
c  INPUT
c  -----
c
c  info    Integer
c          A integer code (set in subroutine GSINS) specifying the error
c
c  jA      Integer
c          No. of the column in matrix A which caused the failure
c
c  jQ      Integer
c          No. of the column in matrix Q to which the orthogonalized column
c          of A would have been written.
c
c  ipr     Logical
c          Flag for output to a file as well as to standard output unit.
c          = .true.   means write duplicate output to a file for later printing,
c          = .false.  means write output only to standard output unit
c
c  mou     Integer
c          Unit number for duplicate output file (not used if ipr = .false.)
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      if ( info .ne. -1 ) go to 50
      print 10, jA, jQ
   10 format (' '/' Reorthog. failed to converge',
     1 ' in processing col.',i3,' into col.',i3,' of orthog. matrix') 
      if (ipr) write (mou,10) jA, jQ
c
c    --------
      return
c    --------
c
   50 continue
      if ( info .ne. +1 ) go to 100
      print 60, jA, jQ
   60 format (' '/' Norm of col.',i3,' is zero so it cannot be',
     1 ' transformed into col.',i3,' of the orthog. matrix')
      if (ipr) write(mou,60) jA, jQ
c
c    --------
      return
c    --------
c
  100 print 110, info
  110 format (' '/' In subr. GSINS  info =',i3,'  which is not',
     1 ' a standard error return.')
      if (ipr) write (mou,110) info
c
c    --------
      return
c    --------
c
      end 
c
c
c*******************************************************************************
c
       subroutine GSDEL ( Q, ldQ, m, k, R, ldR, jdr, info )
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      integer  ldQ, m, k, ldR, jdr, info
c
      real  Q(ldQ,k), R(ldR,k)
c
      external  SCOPY, SROTG, SROT
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c  A subroutine for downdating (deleting a column from) a QR factorization
c
c  Written by Dianne P. O'leary, September 1981.
c
c  Given  A = Q R    with  A: m x k , 
c                          Q: m x k , with orthogonal columns, and
c                          R: k x k , upper triangular,
c    and with
c               A = ( A_1 | a_jdr | A_2 )
c
c               R = ( R_1 | r_jdr | R_2 )
c
c    where:     a_jdr is the column of A to be dropped from the basis, and
c               r_jdr is the corresponding column of R
c
c  Want  ( A_1 | A_2 ) = Q' R'
c
c                     with  Q': m x (k-1) , with orthogonal columns, and
c                           R': (k-1) x (k-1) , upper triangular
c
c-------------------------------------------------------------------------------
c
c  INPUT
c  -----
c
c  Q       Real matrix ( m x k )
c          The old orthogonal matrix
c
c  ldQ     Integer
c          Row dimension of Fortran array Q in the calling program
c
c  m       Integer
c          No. rows in matrix Q 
c
c  k       Integer
c          No. columns in matrix and Q and no. rows and columns in R
c
c  R       Real matrix ( k x k )
c          The old upper triangular matrix
c
c  ldR     Integer
c          Row dimension of Fortran array R in the calling program
c
c  jdr     Integer
c          Index of column of  R  to be dropped
c
c-------------------------------------------------------------------------------
c
c  OUTPUT
c  ------
c  
c  Q       Real matrix ( m x k-1 )
c          Downdated orthogonal matrix
c
c  R       Real matrix ( k-1 x k-1 )
c          Downdated upper triangular matrix 
c
c  info    Integer
c           = 0  if no trouble
c           = 1  if trouble (jdr > k)
c
c-------------------------------------------------------------------------------
c
c  USAGE
c  -----
c 
c     Note:  The value of k is not changed by this subroutine.
c 
c     To delete a column, call GSDEL and then set  k = k - 1
c
c-------------------------------------------------------------------------------
c
c  LOCAL VARIABLES
c  ----- ---------
c
      integer  j
c
      real  s, c
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
       info = 0
       if ( jdr .gt. k )  info = 1
       if ( jdr .gt. k )  return
c
c  Shift the  R  matrix to delete column  jdr
c
       do 40 j = jdr+1,k
           call SCOPY ( j, R(1,j), 1, R(1,j-1), 1 )
   40  continue
c
c  Apply rotations to convert the resulting upper Hessenberg  R
c                    to upper triangular
c
       do 80 j = jdr,k-1
           call SROTG ( R(j,j), R(j+1,j), c, s )
           R(j+1,j) = 0.
           call SROT ( k-j, R(j,j+1), ldR, R(j+1,j+1), ldR, c, s )
           call SROT ( m, Q(1,j), 1, Q(1,j+1), 1, c, s )
   80  continue
c
c     --------
       return
c     --------
c
       end
c
c
c*******************************************************************************
c
      subroutine GSDERR ( jd, k, ipr, mou )
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      logical ipr
c
      integer  jd, k, mou
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c   This subroutine prints error messages from GSDEL.  It is called by LS if
c   subroutine GSDEL fails trying to drop column jd from the factorization.
c
c-------------------------------------------------------------------------------
c
c  INPUT
c  -----
c
c  jd      Integer
c          No. of the column that GSDEL was attempting to drop
c
c  k       Integer
c          No. of columns in the Q and R matrices for the QR factorization
c 
c  ipr     Logical
c          Flag for output to a file as well as to standard output unit.
c          = .true.   means write duplicate output to a file for later printing,
c          = .false.  means write output only to standard output unit
c
c  mou     Integer
c          Unit number for duplicate output file (not used if ipr = .false.)
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      print 10, jd,k
   10 format (' '//' Col.',i3,' cannot be eliminated from the',
     1 ' matrix because it has only ',i3,' columns'/)
      if (ipr) write (mou,10) jd,k
c
c    --------
      return
c    --------
c
      end 
