SUBROUTINE SCHUD(R,LDR,P,X,Z,LDZ,NZ,Y,RHO,C,S)
      INTEGER LDR,P,LDZ,NZ
      REAL RHO(1),C(1)
      REAL R(LDR,1),X(1),Z(LDZ,1),Y(1),S(1)
C
C     SCHUD UPDATES AN AUGMENTED CHOLESKY DECOMPOSITION OF THE
C     TRIANGULAR PART OF AN AUGMENTED QR DECOMPOSITION.  SPECIFICALLY,
C     GIVEN AN UPPER TRIANGULAR MATRIX R OF ORDER P, A ROW VECTOR
C     X, A COLUMN VECTOR Z, AND A SCALAR Y, SCHUD DETERMINES A
C     UNTIARY MATRIX U AND A SCALAR ZETA SUCH THAT
C
C
C                              (R  Z)     (RR   ZZ )
C                         U  * (    )  =  (        ) ,
C                              (X  Y)     ( 0  ZETA)
C
C     WHERE RR IS UPPER TRIANGULAR.  IF R AND Z HAVE BEEN
C     OBTAINED FROM THE FACTORIZATION OF A LEAST SQUARES
C     PROBLEM, THEN RR AND ZZ ARE THE FACTORS CORRESPONDING TO
C     THE PROBLEM WITH THE OBSERVATION (X,Y) APPENDED.  IN THIS
C     CASE, IF RHO IS THE NORM OF THE RESIDUAL VECTOR, THEN THE
C     NORM OF THE RESIDUAL VECTOR OF THE UPDATED PROBLEM IS
C     SQRT(RHO**2 + ZETA**2).  SCHUD WILL SIMULTANEOUSLY UPDATE
C     SEVERAL TRIPLETS (Z,Y,RHO).
C     FOR A LESS TERSE DESCRIPTION OF WHAT SCHUD DOES AND HOW
C     IT MAY BE APPLIED, SEE THE LINPACK GUIDE.
C
C     THE MATRIX U IS DETERMINED AS THE PRODUCT U(P)*...*U(1),
C     WHERE U(I) IS A ROTATION IN THE (I,P+1) PLANE OF THE
C     FORM
C
C                       (     C(I)      S(I) )
C                       (                    ) .
C                       (    -S(I)      C(I) )
C
C     THE ROTATIONS ARE CHOSEN SO THAT C(I) IS REAL.
C
C     ON ENTRY
C
C         R      REAL(LDR,P), WHERE LDR .GE. P.
C                R CONTAINS THE UPPER TRIANGULAR MATRIX
C                THAT IS TO BE UPDATED.  THE PART OF R
C                BELOW THE DIAGONAL IS NOT REFERENCED.
C
C         LDR    INTEGER.
C                LDR IS THE LEADING DIMENSION OF THE ARRAY R.
C
C         P      INTEGER.
C                P IS THE ORDER OF THE MATRIX R.
C
C         X      REAL(P).
C                X CONTAINS THE ROW TO BE ADDED TO R.  X IS
C                NOT ALTERED BY SCHUD.
C
C         Z      REAL(LDZ,NZ), WHERE LDZ .GE. P.
C                Z IS AN ARRAY CONTAINING NZ P-VECTORS TO
C                BE UPDATED WITH R.
C
C         LDZ    INTEGER.
C                LDZ IS THE LEADING DIMENSION OF THE ARRAY Z.
C
C         NZ     INTEGER.
C                NZ IS THE NUMBER OF VECTORS TO BE UPDATED
C                NZ MAY BE ZERO, IN WHICH CASE Z, Y, AND RHO
C                ARE NOT REFERENCED.
C
C         Y      REAL(NZ).
C                Y CONTAINS THE SCALARS FOR UPDATING THE VECTORS
C                Z.  Y IS NOT ALTERED BY SCHUD.
C
C         RHO    REAL(NZ).
C                RHO CONTAINS THE NORMS OF THE RESIDUAL
C                VECTORS THAT ARE TO BE UPDATED.  IF RHO(J)
C                IS NEGATIVE, IT IS LEFT UNALTERED.
C
C     ON RETURN
C
C         RC
C         RHO    CONTAIN THE UPDATED QUANTITIES.
C         Z
C
C         C      REAL(P).
C                C CONTAINS THE COSINES OF THE TRANSFORMING
C                ROTATIONS.
C
C         S      REAL(P).
C                S CONTAINS THE SINES OF THE TRANSFORMING
C                ROTATIONS.
C
C     LINPACK. THIS VERSION DATED 08/14/78 .
C     G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
C
C     SCHUD USES THE FOLLOWING FUNCTIONS AND SUBROUTINES.
C
C     EXTENDED BLAS SROTG
C     FORTRAN SQRT
C
      INTEGER I,J,JM1
      REAL AZETA,SCALE
      REAL T,XJ,ZETA
C
C     UPDATE R.
C
      DO 30 J = 1, P
         XJ = X(J)
C
C        APPLY THE PREVIOUS ROTATIONS.
C
         JM1 = J - 1
         IF (JM1 .LT. 1) GO TO 20
         DO 10 I = 1, JM1
            T = C(I)*R(I,J) + S(I)*XJ
            XJ = C(I)*XJ - S(I)*R(I,J)
            R(I,J) = T
   10    CONTINUE
   20    CONTINUE
C
C        COMPUTE THE NEXT ROTATION.
C
         CALL SROTG(R(J,J),XJ,C(J),S(J))
   30 CONTINUE
C
C     IF REQUIRED, UPDATE Z AND RHO.
C
      IF (NZ .LT. 1) GO TO 70
      DO 60 J = 1, NZ
         ZETA = Y(J)
         DO 40 I = 1, P
            T = C(I)*Z(I,J) + S(I)*ZETA
            ZETA = C(I)*ZETA - S(I)*Z(I,J)
            Z(I,J) = T
   40    CONTINUE
         AZETA = ABS(ZETA)
         IF (AZETA .EQ. 0.0E0 .OR. RHO(J) .LT. 0.0E0) GO TO 50
            SCALE = AZETA + RHO(J)
            RHO(J) = SCALE*SQRT((AZETA/SCALE)**2+(RHO(J)/SCALE)**2)
   50    CONTINUE
   60 CONTINUE
   70 CONTINUE
      RETURN
      END