      SUBROUTINE DCHUR1(UPLO,TRANS,N,M,R,LDR,X,Z,LDZ,Y,RHO,C,S,INFO)
*     .. Scalar Arguments ..
      CHARACTER UPLO,TRANS
      INTEGER N,M,LDR,LDZ,INFO
*     .. Array Arguments ..
      DOUBLE PRECISION R(LDR,*)
      DOUBLE PRECISION X(*)
      DOUBLE PRECISION Z(LDZ,*)
      DOUBLE PRECISION Y(*)
      DOUBLE PRECISION RHO(*)
      DOUBLE PRECISION C(*)
      DOUBLE PRECISION S(*)
*
*  Purpose
*  =======
*
*  DCHUR1 performs a rank 1 update of a Cholesky factor R and 
*  optionnally, of associated vectors Z
*
*  If UPLO = 'U', this solves the problems
*      R2'*R2 = R1'*R1 - x*x'
*      R2'*op(Z2) = R1'*op(Z1) - x*y'
*  Or if UPLO='L', the transposed versions
*      R2*R2'=R1*R1' - x*x'
*      R2*op(Z2)=R1*op(Z1) - x*y'
*  where R1 is R on entry, R2 is R on exit, ' means transpose,
*  Z1 is Z on entry, Z2 is Z on exit, op is transpose or identity.
*
*  The problem can also be written if UPLO='U' as finding U such that
*
*                              (R  op(Z))     (R2  op(Z2))
*                         U  * (        )  =  (          ) ,
*                              (X      Y)     ( 0     Y2 )
*
*  Or its transposed form if UPLO='L'
*
*                         (R      X)          (R2      0 )
*                         (        ) * U'  =  (          ) ,
*                         (op(Z)' Y)          (op(Z2)' Y2)
*
*  Where U is a product of Givens plane rotations U=U(N)*...*U(1)
*  And U(I) is a rotation in the plane (I,N+1)
*
*                        (     c(i)      s(i) )
*                  U(I)= (                    ) .
*                        (    -s(i)      c(i) )
*
*  When op(Z) are the factors of a least square problems,
*  and RHO is the norm of residuals of least square solution,
*  this is equivalent to adding a new observation (X,Y),
*  and update the factored solution R,Z,RHO of the problem.
*
*  Arguments
*  =========
*
*  UPLO    (input) CHARACTER*1
*          = 'U':  R is a upper triangular factor of R'*R
*          = 'L':  R is a lower triangular factor of R*R'.
*
*  TRANS    (input) CHARACTER*1
*          = 'N':  Z is a N-BY-M Matrix and op=identity.
*          = 'T':  Z is a M-BY-N Matrix and op=transpose
*
*  N       (input) INTEGER
*          The order of the matrix R.  N >= 0.
*
*  M       (input) INTEGER
*          The number of columns the matrix Z.  M >= 0.
*          In the case M = 0, the array Z,Y,RHO are not referenced
*
*  R       (input/output) DOUBLE PRECISION array, dimension (LDR,N)
*          On entry, the triangular matrix R.  If UPLO = 'U', the leading
*          N-by-N upper triangular part of R contains the upper
*          triangular part of the matrix R, and the strictly lower
*          triangular part of R is not referenced.  If UPLO = 'L', the
*          leading N-by-N lower triangular part of R contains the lower
*          triangular part of the matrix R, and the strictly upper
*          triangular part of R is not referenced.
*
*          On exit, if INFO = 0, R contains the update Cholesky factor
*
*  LDR     (input) INTEGER
*          The leading dimension of the array R, LDR >= max(1,N).
*
*  X       (input/output) DOUBLE PRECISION array, dimension (N)
*          On entry, X contains the vector X used to update R.
*          On exit, X is destructed.
*
*  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,k)
*          On entry, Z contains a list of M vectors of dimension N
*          to be updated with R.
*          If TRANS = 'N', Z is considered to be N-by-M (column vectors)
*          If TRANS = 'T', Z is considered to be M-by-N (row vectors)
*          On exit, if INFO = 0, Z contains the updated matrix Z
*
*  LDZ     (input) INTEGER
*          The leading dimension of the array Z.
*          If TRANS = 'N' LDZ >= max(1,N).
*          If TRANS = 'T' LDZ >= max(1,M).
*
*  Y       (input/output) DOUBLE PRECISION array, dimension (M)
*          On entry, Y contains the scalar quantities used to update Z.
*          on exit, Y contains the updated values Y2.
*
*  RHO     (input/output) DOUBLE PRECISION array, dimension (M)
*          on entry, RHO contains the norm of residual vectors.
*          on exit, RHO contains the updated norm of residual vectors.
*
*  C       (output) DOUBLE PRECISION array, dimension (N)
*          C contains the cosines of performed Givens rotations.
*
*  S       (output) DOUBLE PRECISION array, dimension (N)
*          S contains the sines of performed Givens rotations.
*
*  INFO    (output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value
*          > 0:  if INFO = i, the ith residual RHO was negative.
*
*
*  Note:
*  This is a variant of LINPACKs' dchud [1] & [2]
*  but with the following differences taking [3] into account
*  - uses BLAS drot and drotg for Given's rotations
*  - forces positive elements on the diagonal of R
*  - accepts upper or lower Cholesky factor R
*
*  Author: Nicolas Cellier
*
*  References:
*  [1] Linpack dchud http://www.netlib.org/linpack/dchud.f
*  [2] LINPACK User Guide http://www.uploading.com/files/F89HJTGH/linguide.zip.html
*  [3] Low Rank Updates for the Cholesky Decomposition - Matthias Seeger
*      http://infoscience.epfl.ch/record/161468/files/cholupdate.pdf
*     
*     .. Local Scalars ..
      LOGICAL UPPER,DOTRAN
      INTEGER I
      DOUBLE PRECISION SC,W
*     .. External Functions ..
      LOGICAL LSAME
      EXTERNAL LSAME
*     .. External Subroutines ..
      EXTERNAL DROT,DROTG
*     .. Intrinsic Functions ..
      INTRINSIC DSQRT
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      INFO=0
      UPPER=LSAME(UPLO,'U')
      DOTRAN=LSAME(TRANS,'T')
      IF( .NOT. UPPER .AND. .NOT. LSAME(UPLO,'L') ) THEN
        INFO=-1
      ELSE IF( .NOT.DOTRAN .AND. .NOT. LSAME(TRANS,'N') ) THEN
        INFO=-2
      ELSE IF( N.LT.0 ) THEN
        INFO=-3
      ELSE IF( M.LT.0 ) THEN
        INFO=-4
      ELSE IF( LDR.LT.N ) THEN
        INFO=-6
      ELSE IF (M.GT.0) THEN
        IF(DOTRAN) THEN
          IF( LDZ.LT.M ) INFO=-9
        ELSE
          IF( LDZ.LT.N ) INFO=-9
        ENDIF
      ENDIF
      IF(INFO.NE.0) RETURN
*
*     Quick return if possible
*
      IF(N.EQ.0) RETURN
*
* ------ FIRST PART: UPDATE R ------
*
      IF(UPPER) THEN 
*       Compute givens rotations to annihilate x
        DO 10 I=1,N
          CALL DROTG(R(I,I), X(I), C(I), S(I))
*         Force positive diagonal values
          IF(R(I,I).LT.0.0D0) THEN
            R(I,I)=-R(I,I)
            C(I)=-C(I)
            S(I)=-S(I)
          ENDIF
*         Apply rotation on each row
          IF(I.LT.N) THEN
            CALL DROT(N-I,R(I,I+1),LDR,X(I+1),1,C(I),S(I))
          ENDIF
 10     CONTINUE
      ELSE
*       Compute givens rotations to annihilate x
        DO 20 I=1,N
          CALL DROTG(R(I,I), X(I), C(I), S(I))
*         Force positive diagonal values
          IF(R(I,I).LT.0.0D0) THEN
            R(I,I)=-R(I,I)
            C(I)=-C(I)
            S(I)=-S(I)
          ENDIF
*         Apply rotation on each column
          IF(I.LT.N) THEN
            CALL DROT(N-I,R(I+1,I),1,X(I+1),1,C(I),S(I))
          ENDIF
 20     CONTINUE
      ENDIF
*
* ------  SECOND PART: UPDATE Z AND RHO IF REQUESTED ------
*
      IF(M.GT.0) THEN
        IF(DOTRAN) THEN 
*         Update colums of z
          DO 30 I=1,N
            CALL DROT(M,Z(1,I),1,Y(1),1,C(I),S(I));
 30       CONTINUE
        ELSE
*         Update rows of z
          DO 40 I=1,N
            CALL DROT(M,Z(I,1),LDZ,Y(1),1,C(I),S(I));
 40       CONTINUE
        ENDIF
*       Update rho
        DO 50 I=1,M
          W=DABS(Y(I))
          IF(RHO(I).GT.0.0D0) THEN
            SC=RHO(I)+W
            RHO(I)=SC*DSQRT((RHO(I)/SC)**2+(W/SC)**2)
          ELSE
            INFO=I
            GOTO 100
          ENDIF
 50     CONTINUE
      ENDIF
 100  CONTINUE
      RETURN
      END
