SUBROUTINE DAXPY( N, DA, DX, INCX, DY, INCY )
*
*     constant times a vector plus a vector.
*     uses unrolled loops for increments equal to one.
*     jack dongarra, linpack, 3/11/78.
*
*     .. Scalar Arguments ..
      INTEGER           INCX, INCY, N
      DOUBLE PRECISION  DA
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION  DX( 1 ), DY( 1 )
*     ..
*     .. Local Scalars ..
      INTEGER           I, IX, IY, M, MP1
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC         MOD
*     ..
*     .. Executable Statements ..
*
      IF( N.LE.0 )
     $   RETURN
      IF( DA.EQ.0.0D0 )
     $   RETURN
      IF( INCX.EQ.1 .AND. INCY.EQ.1 )
     $   GO TO 20
*
*        code for unequal increments or equal increments
*          not equal to 1
*
      IX = 1
      IY = 1
      IF( INCX.LT.0 )
     $   IX = ( -N+1 )*INCX + 1
      IF( INCY.LT.0 )
     $   IY = ( -N+1 )*INCY + 1
      DO 10 I = 1, N
         DY( IY ) = DY( IY ) + DA*DX( IX )
         IX = IX + INCX
         IY = IY + INCY
   10 CONTINUE
      RETURN
*
*        code for both increments equal to 1
*
*
*        clean-up loop
*
   20 M = MOD( N, 4 )
      IF( M.EQ.0 )
     $   GO TO 40
      DO 30 I = 1, M
         DY( I ) = DY( I ) + DA*DX( I )
   30 CONTINUE
      IF( N.LT.4 )
     $   RETURN
   40 MP1 = M + 1
      DO 50 I = MP1, N, 4
         DY( I ) = DY( I ) + DA*DX( I )
         DY( I+1 ) = DY( I+1 ) + DA*DX( I+1 )
         DY( I+2 ) = DY( I+2 ) + DA*DX( I+2 )
         DY( I+3 ) = DY( I+3 ) + DA*DX( I+3 )
   50 CONTINUE
      RETURN
      END
      SUBROUTINE DCOPY( N, DX, INCX, DY, INCY )
*
*     copies a vector, x, to a vector, y.
*     uses unrolled loops for increments equal to one.
*     jack dongarra, linpack, 3/11/78.
*
*     .. Scalar Arguments ..
      INTEGER           INCX, INCY, N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION  DX( 1 ), DY( 1 )
*     ..
*     .. Local Scalars ..
      INTEGER           I, IX, IY, M, MP1
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC         MOD
*     ..
*     .. Executable Statements ..
*
      IF( N.LE.0 )
     $   RETURN
      IF( INCX.EQ.1 .AND. INCY.EQ.1 )
     $   GO TO 20
*
*        code for unequal increments or equal increments
*          not equal to 1
*
      IX = 1
      IY = 1
      IF( INCX.LT.0 )
     $   IX = ( -N+1 )*INCX + 1
      IF( INCY.LT.0 )
     $   IY = ( -N+1 )*INCY + 1
      DO 10 I = 1, N
         DY( IY ) = DX( IX )
         IX = IX + INCX
         IY = IY + INCY
   10 CONTINUE
      RETURN
*
*        code for both increments equal to 1
*
*
*        clean-up loop
*
   20 M = MOD( N, 7 )
      IF( M.EQ.0 )
     $   GO TO 40
      DO 30 I = 1, M
         DY( I ) = DX( I )
   30 CONTINUE
      IF( N.LT.7 )
     $   RETURN
   40 MP1 = M + 1
      DO 50 I = MP1, N, 7
         DY( I ) = DX( I )
         DY( I+1 ) = DX( I+1 )
         DY( I+2 ) = DX( I+2 )
         DY( I+3 ) = DX( I+3 )
         DY( I+4 ) = DX( I+4 )
         DY( I+5 ) = DX( I+5 )
         DY( I+6 ) = DX( I+6 )
   50 CONTINUE
      RETURN
      END
      DOUBLE PRECISION FUNCTION DDOT( N, DX, INCX, DY, INCY )
*
*     forms the dot product of two vectors.
*     uses unrolled loops for increments equal to one.
*     jack dongarra, linpack, 3/11/78.
*
*     .. Scalar Arguments ..
      INTEGER                         INCX, INCY, N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION                DX( 1 ), DY( 1 )
*     ..
*     .. Local Scalars ..
      INTEGER                         I, IX, IY, M, MP1
      DOUBLE PRECISION                DTEMP
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC                       MOD
*     ..
*     .. Executable Statements ..
*
      DDOT = 0.0D0
      DTEMP = 0.0D0
      IF( N.LE.0 )
     $   RETURN
      IF( INCX.EQ.1 .AND. INCY.EQ.1 )
     $   GO TO 20
*
*        code for unequal increments or equal increments
*          not equal to 1
*
      IX = 1
      IY = 1
      IF( INCX.LT.0 )
     $   IX = ( -N+1 )*INCX + 1
      IF( INCY.LT.0 )
     $   IY = ( -N+1 )*INCY + 1
      DO 10 I = 1, N
         DTEMP = DTEMP + DX( IX )*DY( IY )
         IX = IX + INCX
         IY = IY + INCY
   10 CONTINUE
      DDOT = DTEMP
      RETURN
*
*        code for both increments equal to 1
*
*
*        clean-up loop
*
   20 M = MOD( N, 5 )
      IF( M.EQ.0 )
     $   GO TO 40
      DO 30 I = 1, M
         DTEMP = DTEMP + DX( I )*DY( I )
   30 CONTINUE
      IF( N.LT.5 )
     $   GO TO 60
   40 MP1 = M + 1
      DO 50 I = MP1, N, 5
         DTEMP = DTEMP + DX( I )*DY( I ) + DX( I+1 )*DY( I+1 ) +
     $           DX( I+2 )*DY( I+2 ) + DX( I+3 )*DY( I+3 ) +
     $           DX( I+4 )*DY( I+4 )
   50 CONTINUE
   60 DDOT = DTEMP
      RETURN
      END
      DOUBLE PRECISION FUNCTION DNRM2( N, DX, INCX )
*
*     euclidean norm of the n-vector stored in dx() with storage
*     increment incx .
*     if    n .le. 0 return with result = 0.
*     if n .ge. 1 then incx must be .ge. 1
*
*           c.l.lawson, 1978 jan 08
*
*     four phase method     using two built-in constants that are
*     hopefully applicable to all machines.
*         cutlo = maximum of  dsqrt(u/eps)  over all known machines.
*         cuthi = minimum of  dsqrt(v)      over all known machines.
*     where
*         eps = smallest no. such that eps + 1. .gt. 1.
*         u   = smallest positive no.   (underflow limit)
*         v   = largest  no.            (overflow  limit)
*
*     brief outline of algorithm..
*
*     phase 1    scans zero components.
*     move to phase 2 when a component is nonzero and .le. cutlo
*     move to phase 3 when a component is .gt. cutlo
*     move to phase 4 when a component is .ge. cuthi/m
*     where m = n for x() real and m = 2*n for complex.
*
*     values for cutlo and cuthi..
*     from the environmental parameters listed in the imsl converter
*     document the limiting values are as follows..
*     cutlo, s.p.   u/eps = 2**(-102) for  honeywell.  close seconds are
*                   univac and dec at 2**(-103)
*                   thus cutlo = 2**(-51) = 4.44089e-16
*     cuthi, s.p.   v = 2**127 for univac, honeywell, and dec.
*                   thus cuthi = 2**(63.5) = 1.30438e19
*     cutlo, d.p.   u/eps = 2**(-67) for honeywell and dec.
*                   thus cutlo = 2**(-33.5) = 8.23181d-11
*     cuthi, d.p.   same as s.p.  cuthi = 1.30438d19
*     data cutlo, cuthi / 8.232d-11,  1.304d19 /
*     data cutlo, cuthi / 4.441e-16,  1.304e19 /
*
*     .. Scalar Arguments ..
      INTEGER                          INCX, N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION                 DX( 1 )
*     ..
*     .. Local Scalars ..
      INTEGER                          I, IX, J, NEXT, NN
      DOUBLE PRECISION                 CUTHI, CUTLO, HITEST, ONE, SUM,
     $                                 XMAX, ZERO
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC                        DABS, DSQRT, FLOAT
*     ..
*     .. Data statements ..
      DATA                             ZERO, ONE / 0.0D0, 1.0D0 /
      DATA                             CUTLO, CUTHI / 8.232D-11,
     $                                 1.304D19 /
*     ..
*     .. Executable Statements ..
*
      IF( N.GT.0 )
     $   GO TO 10
      DNRM2 = ZERO
      GO TO 140
*
   10 ASSIGN 30 TO NEXT
      SUM = ZERO
*
*        begin main loop
*
      IX = 1
      IF( INCX.LT.0 )
     $   IX = 1 - ( N-1 )*INCX
      I = IX
      NN = IX + ( N-1 )*INCX
   20 GO TO NEXT( 30, 40, 70, 80 )
   30 IF( DABS( DX( I ) ).GT.CUTLO )
     $   GO TO 110
      ASSIGN 40 TO NEXT
      XMAX = ZERO
*
*        phase 1.  sum is zero
*
   40 IF( DX( I ).EQ.ZERO )
     $   GO TO 130
      IF( DABS( DX( I ) ).GT.CUTLO )
     $   GO TO 110
*
*        prepare for phase 2.
*
      ASSIGN 70 TO NEXT
      GO TO 60
*
*        prepare for phase 4.
*
   50 I = J
      ASSIGN 80 TO NEXT
      SUM = ( SUM / DX( I ) ) / DX( I )
   60 XMAX = DABS( DX( I ) )
      GO TO 90
*
*        phase 2.  sum is small.
*                  scale to avoid destructive underflow.
*
   70 IF( DABS( DX( I ) ).GT.CUTLO )
     $   GO TO 100
*
*        common code for phases 2 and 4.
*        in phase 4 sum is large.  scale to avoid overflow.
*
   80 IF( DABS( DX( I ) ).LE.XMAX )
     $   GO TO 90
      SUM = ONE + SUM*( XMAX / DX( I ) )**2
      XMAX = DABS( DX( I ) )
      GO TO 130
*
   90 SUM = SUM + ( DX( I ) / XMAX )**2
      GO TO 130
*
*        prepare for phase 3.
*
  100 SUM = ( SUM*XMAX )*XMAX
*
*        for real or d.p. set hitest = cuthi/n
*        for complex      set hitest = cuthi/(2*n)
*
  110 HITEST = CUTHI / FLOAT( N )
*
*        phase 3.  sum is mid-range.  no scaling.
*
      DO 120 J = I, NN, INCX
         IF( DABS( DX( J ) ).GE.HITEST )
     $      GO TO 50
         SUM = SUM + DX( J )**2
  120 CONTINUE
      DNRM2 = DSQRT( SUM )
      GO TO 140
*
  130 CONTINUE
      IF( I.NE.NN ) THEN
         I = I + INCX
         GO TO 20
      ENDIF
*
*        end of main loop.
*
*        compute square root and adjust for scaling.
*
      DNRM2 = XMAX*DSQRT( SUM )
  140 CONTINUE
      RETURN
      END
      SUBROUTINE DROT( N, DX, INCX, DY, INCY, C, S )
*
*     applies a plane rotation.
*     jack dongarra, linpack, 3/11/78.
*
*     .. Scalar Arguments ..
      INTEGER          INCX, INCY, N
      DOUBLE PRECISION C, S
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION DX( 1 ), DY( 1 )
*     ..
*     .. Local Scalars ..
      INTEGER          I, IX, IY
      DOUBLE PRECISION DTEMP
*     ..
*     .. Executable Statements ..
*
      IF( N.LE.0 )
     $   RETURN
      IF( INCX.EQ.1 .AND. INCY.EQ.1 )
     $   GO TO 20
*
*        code for unequal increments or equal increments
*          not equal to 1
*
      IX = 1
      IY = 1
      IF( INCX.LT.0 )
     $   IX = ( -N+1 )*INCX + 1
      IF( INCY.LT.0 )
     $   IY = ( -N+1 )*INCY + 1
      DO 10 I = 1, N
         DTEMP = C*DX( IX ) + S*DY( IY )
         DY( IY ) = C*DY( IY ) - S*DX( IX )
         DX( IX ) = DTEMP
         IX = IX + INCX
         IY = IY + INCY
   10 CONTINUE
      RETURN
*
*        code for both increments equal to 1
*
   20 DO 30 I = 1, N
         DTEMP = C*DX( I ) + S*DY( I )
         DY( I ) = C*DY( I ) - S*DX( I )
         DX( I ) = DTEMP
   30 CONTINUE
      RETURN
      END
      SUBROUTINE DROTG( DA, DB, C, S )
*
*     construct givens plane rotation.
*     jack dongarra, linpack, 3/11/78.
*                    modified 9/27/86.
*
*     .. Scalar Arguments ..
      DOUBLE PRECISION  C, DA, DB, S
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION  R, ROE, SCALE, Z
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC         DABS, DSIGN, DSQRT
*     ..
*     .. Executable Statements ..
*
      ROE = DB
      IF( DABS( DA ).GT.DABS( DB ) )
     $   ROE = DA
      SCALE = DABS( DA ) + DABS( DB )
      IF( SCALE.NE.0.0D0 )
     $   GO TO 10
      C = 1.0D0
      S = 0.0D0
      R = 0.0D0
      GO TO 20
   10 R = SCALE*DSQRT( ( DA / SCALE )**2+( DB / SCALE )**2 )
      R = DSIGN( 1.0D0, ROE )*R
      C = DA / R
      S = DB / R
   20 Z = S
      IF( DABS( C ).GT.0.0D0 .AND. DABS( C ).LE.S )
     $   Z = 1.0D0 / C
      DA = R
      DB = Z
      RETURN
      END
      SUBROUTINE DSCAL( N, DA, DX, INCX )
*
*     scales a vector by a constant.
*     uses unrolled loops for increment equal to one.
*     jack dongarra, linpack, 3/11/78.
*
*     .. Scalar Arguments ..
      INTEGER           INCX, N
      DOUBLE PRECISION  DA
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION  DX( 1 )
*     ..
*     .. Local Scalars ..
      INTEGER           I, IX, M, MP1, NINCX
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC         MOD
*     ..
*     .. Executable Statements ..
*
      IF( N.LE.0 )
     $   RETURN
      IF( INCX.EQ.1 )
     $   GO TO 20
*
*        code for increment not equal to 1
*
      IX = 1
      IF( INCX.LT.0 )
     $   IX = 1 - ( N-1 )*INCX
      NINCX = IX + ( N-1 )*INCX
      DO 10 I = IX, NINCX, INCX
         DX( I ) = DA*DX( I )
   10 CONTINUE
      RETURN
*
*        code for increment equal to 1
*
*
*        clean-up loop
*
   20 M = MOD( N, 5 )
      IF( M.EQ.0 )
     $   GO TO 40
      DO 30 I = 1, M
         DX( I ) = DA*DX( I )
   30 CONTINUE
      IF( N.LT.5 )
     $   RETURN
   40 MP1 = M + 1
      DO 50 I = MP1, N, 5
         DX( I ) = DA*DX( I )
         DX( I+1 ) = DA*DX( I+1 )
         DX( I+2 ) = DA*DX( I+2 )
         DX( I+3 ) = DA*DX( I+3 )
         DX( I+4 ) = DA*DX( I+4 )
   50 CONTINUE
      RETURN
      END
      SUBROUTINE DSWAP( N, DX, INCX, DY, INCY )
*
*     interchanges two vectors.
*     uses unrolled loops for increments equal one.
*     jack dongarra, linpack, 3/11/78.
*
*     .. Scalar Arguments ..
      INTEGER           INCX, INCY, N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION  DX( 1 ), DY( 1 )
*     ..
*     .. Local Scalars ..
      INTEGER           I, IX, IY, M, MP1
      DOUBLE PRECISION  DTEMP
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC         MOD
*     ..
*     .. Executable Statements ..
*
      IF( N.LE.0 )
     $   RETURN
      IF( INCX.EQ.1 .AND. INCY.EQ.1 )
     $   GO TO 20
*
*        code for unequal increments or equal increments
*          not equal to 1
*
      IX = 1
      IY = 1
      IF( INCX.LT.0 )
     $   IX = ( -N+1 )*INCX + 1
      IF( INCY.LT.0 )
     $   IY = ( -N+1 )*INCY + 1
      DO 10 I = 1, N
         DTEMP = DX( IX )
         DX( IX ) = DY( IY )
         DY( IY ) = DTEMP
         IX = IX + INCX
         IY = IY + INCY
   10 CONTINUE
      RETURN
*
*        code for both increments equal to 1
*
*
*        clean-up loop
*
   20 M = MOD( N, 3 )
      IF( M.EQ.0 )
     $   GO TO 40
      DO 30 I = 1, M
         DTEMP = DX( I )
         DX( I ) = DY( I )
         DY( I ) = DTEMP
   30 CONTINUE
      IF( N.LT.3 )
     $   RETURN
   40 MP1 = M + 1
      DO 50 I = MP1, N, 3
         DTEMP = DX( I )
         DX( I ) = DY( I )
         DY( I ) = DTEMP
         DTEMP = DX( I+1 )
         DX( I+1 ) = DY( I+1 )
         DY( I+1 ) = DTEMP
         DTEMP = DX( I+2 )
         DX( I+2 ) = DY( I+2 )
         DY( I+2 ) = DTEMP
   50 CONTINUE
      RETURN
      END
      DOUBLE PRECISION FUNCTION DASUM(N,DX,INCX)
C
C     TAKES THE SUM OF THE ABSOLUTE VALUES.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DX(1),DTEMP
      INTEGER I,INCX,M,MP1,N,NINCX
C
      DASUM = 0.0D0
      DTEMP = 0.0D0
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1)GO TO 20
C
C        CODE FOR INCREMENT NOT EQUAL TO 1
C
      NINCX = N*INCX
      DO 10 I = 1,NINCX,INCX
        DTEMP = DTEMP + DABS(DX(I))
   10 CONTINUE
      DASUM = DTEMP
      RETURN
C
C        CODE FOR INCREMENT EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,6)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        DTEMP = DTEMP + DABS(DX(I))
   30 CONTINUE
      IF( N .LT. 6 ) GO TO 60
   40 MP1 = M + 1
      DO 50 I = MP1,N,6
        DTEMP = DTEMP + DABS(DX(I)) + DABS(DX(I + 1)) + DABS(DX(I + 2))
     *  + DABS(DX(I + 3)) + DABS(DX(I + 4)) + DABS(DX(I + 5))
   50 CONTINUE
   60 DASUM = DTEMP
      RETURN
      END