SUBROUTINE SLARTG( F, G, CS, SN, R )
*
*  -- LAPACK auxiliary routine (version 3.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     Jan 17, 2001
*
*     .. Scalar Arguments ..
      REAL               CS, F, G, R, SN
*     ..
*
*  Purpose
*  =======
*
*  SLARTG generate a plane rotation so that
*
*     [  CS  SN  ]  .  [ F ]  =  [ R ]   where CS**2 + SN**2 = 1.
*     [ -SN  CS  ]     [ G ]     [ 0 ]
*
*  This is a slower, more accurate version of the BLAS1 routine SROTG,
*  with the following other differences:
*     F and G are unchanged on return.
*     If F=0 and G=0, then CS=1, SN=0, and R=0.
*     If F .ne. 0 and G=0, then CS=1, SN=0, and R=F.
*     If F=0 and G .ne. 0, then CS=0, SN=sign(G), and R=abs(G).
*     If F .ne. 0 and (G .ne. 0), then 
*             CS = abs(F)/sqrt(F**2 + G**2)
*             SN = sign(F)*G/sqrt(F**2 + G**2)
*             R  = sign(F)*sqrt(F**2 + G**2)
*
*  If a NaN occurs in the input, then R and possibly also
*  CS and SN will be a NaN.
*  
*  If an infinity occurs in the input, then R and possibly also
*  CS and SN will be infinite or NaNs.
*  
*  The complex routine CLARTG returns the same 
*  CS and SN on complex inputs (F,0) and (G,0).
*
*  Arguments
*  =========
*
*  F       (input) REAL
*          The first component of vector to be rotated.
*
*  G       (input) REAL
*          The second component of vector to be rotated.
*
*  CS      (output) REAL
*          The cosine of the rotation.
*
*  SN      (output) REAL
*          The sine of the rotation.
*
*  R       (output) REAL
*          The nonzero component of the rotated vector.
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO
      PARAMETER          ( ZERO = 0.0E0 )
      REAL               ONE
      PARAMETER          ( ONE = 1.0E0 )
      REAL               TWO
      PARAMETER          ( TWO = 2.0E0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            FIRST
      INTEGER            COUNT, I, MAXCNT
      REAL               EPS, F1, G1, SAFMIN, SAFMN2, SAFMX2, SCALE
      REAL               SAFMN, SAFMX, SCL, ESFMN2, RR

*     ..
*     .. External Functions ..
      REAL               SLAMCH
      EXTERNAL           SLAMCH
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, INT, LOG, MAX, SQRT, SIGN
*     ..
*     .. Save statement ..
      SAVE               FIRST, EPS, SAFMX2, SAFMIN, SAFMN2, SAFMN
      SAVE               SAFMX
*     ..
*     .. Data statements ..
      DATA               FIRST / .TRUE. /
*     ..
*     .. Executable Statements ..
*
      IF( FIRST ) THEN
*
*        On first call to SLARTG, compute
*        SAFMN2 = sqrt(SAFMIN/EPS) rounded down to the nearest power
*                 of the floating point radix
*        This means that scaling by multiplication by SAFMN2 and its 
*        reciprocal SAFMX2 cause no roundoff error
*
         FIRST = .FALSE.
         SAFMIN = SLAMCH( 'S' )
         EPS = SLAMCH( 'E' )
         ESFMN2 = INT( LOG( SAFMIN / EPS ) / LOG( SLAMCH( 'B' ) ) / TWO)
         SAFMN2 = SLAMCH( 'B' )**ESFMN2
         SAFMN = SAFMN2**2
         SAFMX2 = ONE / SAFMN2
         SAFMX = SAFMX2**2
         MAXCNT = INT( -MAX( SLAMCH('L'), SLAMCH('N')+ONE-SLAMCH('M') )
     +            /ESFMN2 - ONE - ONE )
      END IF
      IF( G.EQ.ZERO ) THEN
*
*        Includes the case F=G=0
*
         CS = ONE
         SN = ZERO
         R = F
      ELSE IF( F.EQ.ZERO ) THEN
*
*        G must be nonzero
*
         CS = ZERO
         SN = SIGN( ONE, G )
         R  = ABS(G)
      ELSE
*
*        Both F and G must be nonzero
*
         F1 = F
         G1 = G
         SCALE = MAX( ABS( F1 ), ABS( G1 ) )
         COUNT = 0
         IF( SCALE.GE.SAFMX2 ) THEN
*
*           Handle case where F1**2 + G1**2 might overflow
*
            SCL = SAFMX2
*
            COUNT = COUNT + 1
            F1 = F1*SAFMN2
            G1 = G1*SAFMN2
            SCALE = SCALE*SAFMN2
            IF( SCALE.LT.SAFMX2 ) GOTO 100
*
            COUNT = COUNT + 1
            F1 = F1*SAFMN2
            G1 = G1*SAFMN2
            SCALE = SCALE*SAFMN2
            IF( SCALE.LT.SAFMX2 ) GOTO 100
*
10          CONTINUE
            COUNT = COUNT + 1
            F1 = F1*SAFMN2
            G1 = G1*SAFMN2
            SCALE = SCALE*SAFMN2
            IF( SCALE.GT.SAFMX2 .AND. COUNT.LE.MAXCNT ) GOTO 10
         ELSE IF( SCALE.LE.SAFMN2 ) THEN
*
*           Handle case where F1**2 + G1**2 might underflow
*
            SCL = SAFMN2
*
            COUNT = COUNT + 1
            F1 = F1*SAFMX2
            G1 = G1*SAFMX2
            SCALE = SCALE*SAFMX2
            IF( SCALE.GT.SAFMN2 ) GOTO 100
*
            COUNT = COUNT + 1
            F1 = F1*SAFMX2
            G1 = G1*SAFMX2
            SCALE = SCALE*SAFMX2
            IF( SCALE.GT.SAFMN2 ) GOTO 100
*
20          CONTINUE
            COUNT = COUNT + 1
            F1 = F1*SAFMX2
            G1 = G1*SAFMX2
            SCALE = SCALE*SAFMX2
            IF( SCALE.LT.SAFMN2 .AND. COUNT.LE.MAXCNT ) GOTO 20
         ENDIF
100      CONTINUE
         R = SQRT( F1**2+G1**2 )
         RR = ONE/R
         CS = ABS(F1) * RR
         SN = G1 * RR
         IF (F .LT. ZERO) THEN
            R = -R
            SN = -SN
         ENDIF
         DO 40 I = 1, COUNT
            R = R*SCL
40       CONTINUE
      ENDIF
      RETURN
*
*     End of SLARTG
*
      END