Index: ThirdParty/mpqc_open/src/lib/math/scmat/pdsteqr.f
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/scmat/pdsteqr.f	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/scmat/pdsteqr.f	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,3143 @@
+* Modified by Curtis Janssen (cljanss@ca.sandia.gov) to update only a
+* portion of the eigenvector matrix.
+      SUBROUTINE PDSTEQR(N, D, E, Z, LDZ, nz, WORK, INFO )
+*
+*  -- LAPACK routine (version 3.0) --
+*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+*     Courant Institute, Argonne National Lab, and Rice University
+*     September 30, 1994
+*
+*     .. Scalar Arguments ..
+      INTEGER            INFO, LDZ, N
+      integer nz
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   D( * ), E( * ), WORK( * ), Z( LDZ, * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DSTEQR computes all eigenvalues and, optionally, eigenvectors of a
+*  symmetric tridiagonal matrix using the implicit QL or QR method.
+*  The eigenvectors of a full or band symmetric matrix can also be found
+*  if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to
+*  tridiagonal form.
+*
+*  Arguments
+*  =========
+*
+*  COMPZ   (input) CHARACTER*1
+*          = 'N':  Compute eigenvalues only.
+*          = 'V':  Compute eigenvalues and eigenvectors of the original
+*                  symmetric matrix.  On entry, Z must contain the
+*                  orthogonal matrix used to reduce the original matrix
+*                  to tridiagonal form.
+*          = 'I':  Compute eigenvalues and eigenvectors of the
+*                  tridiagonal matrix.  Z is initialized to the identity
+*                  matrix.
+*
+*  N       (input) INTEGER
+*          The order of the matrix.  N >= 0.
+*
+*  D       (input/output) DOUBLE PRECISION array, dimension (N)
+*          On entry, the diagonal elements of the tridiagonal matrix.
+*          On exit, if INFO = 0, the eigenvalues in ascending order.
+*
+*  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
+*          On entry, the (n-1) subdiagonal elements of the tridiagonal
+*          matrix.
+*          On exit, E has been destroyed.
+*
+*  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
+*          On entry, if  COMPZ = 'V', then Z contains the orthogonal
+*          matrix used in the reduction to tridiagonal form.
+*          On exit, if INFO = 0, then if  COMPZ = 'V', Z contains the
+*          orthonormal eigenvectors of the original symmetric matrix,
+*          and if COMPZ = 'I', Z contains the orthonormal eigenvectors
+*          of the symmetric tridiagonal matrix.
+*          If COMPZ = 'N', then Z is not referenced.
+*
+*  LDZ     (input) INTEGER
+*          The leading dimension of the array Z.  LDZ >= 1, and if
+*          eigenvectors are desired, then  LDZ >= max(1,N).
+*
+*  WORK    (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2))
+*          If COMPZ = 'N', then WORK is not referenced.
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value
+*          > 0:  the algorithm has failed to find all the eigenvalues in
+*                a total of 30*N iterations; if INFO = i, then i
+*                elements of E have not converged to zero; on exit, D
+*                and E contain the elements of a symmetric tridiagonal
+*                matrix which is orthogonally similar to the original
+*                matrix.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO, ONE, TWO, THREE
+      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
+     $                   THREE = 3.0D0 )
+      INTEGER            MAXIT
+      PARAMETER          ( MAXIT = 30 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND,
+     $                   LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1,
+     $                   NM1, NMAXIT
+      DOUBLE PRECISION   ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2,
+     $                   S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST
+*     ..
+*     .. External Functions ..
+      LOGICAL            PLSAME
+      DOUBLE PRECISION   PDLAMCH, PDLANST, PDLAPY2
+      EXTERNAL           PLSAME, PDLAMCH, PDLANST, PDLAPY2
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           PDLAE2,PDLAEV2,PDLARTG,PDLASCL,PDLASET,PDLASR,
+     $                   PDLASRT, DSWAP, PXERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, MAX, SIGN, SQRT
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input parameters.
+*
+      INFO = 0
+*
+      ICOMPZ = 1
+      IF( ICOMPZ.LT.0 ) THEN
+         INFO = -1
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1,
+     $         nz ) ) ) THEN
+         INFO = -6
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL PXERBLA( 'DSTEQR', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
+      IF( N.EQ.1 ) THEN
+         RETURN
+      END IF
+*
+*     Determine the unit roundoff and over/underflow thresholds.
+*
+      EPS = PDLAMCH( 'E' )
+      EPS2 = EPS**2
+      SAFMIN = PDLAMCH( 'S' )
+      SAFMAX = ONE / SAFMIN
+      SSFMAX = SQRT( SAFMAX ) / THREE
+      SSFMIN = SQRT( SAFMIN ) / EPS2
+*
+*     Compute the eigenvalues and eigenvectors of the tridiagonal
+*     matrix.
+*
+      IF( ICOMPZ.EQ.2 )
+     $   CALL PDLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
+*
+      NMAXIT = N*MAXIT
+      JTOT = 0
+*
+*     Determine where the matrix splits and choose QL or QR iteration
+*     for each block, according to whether top or bottom diagonal
+*     element is smaller.
+*
+      L1 = 1
+      NM1 = N - 1
+*
+   10 CONTINUE
+      IF( L1.GT.N )
+     $   GO TO 160
+      IF( L1.GT.1 )
+     $   E( L1-1 ) = ZERO
+      IF( L1.LE.NM1 ) THEN
+         DO 20 M = L1, NM1
+            TST = ABS( E( M ) )
+            IF( TST.EQ.ZERO )
+     $         GO TO 30
+            IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
+     $          1 ) ) ) )*EPS ) THEN
+               E( M ) = ZERO
+               GO TO 30
+            END IF
+   20    CONTINUE
+      END IF
+      M = N
+*
+   30 CONTINUE
+      L = L1
+      LSV = L
+      LEND = M
+      LENDSV = LEND
+      L1 = M + 1
+      IF( LEND.EQ.L )
+     $   GO TO 10
+*
+*     Scale submatrix in rows and columns L to LEND
+*
+      ANORM = PDLANST( 'I', LEND-L+1, D( L ), E( L ) )
+      ISCALE = 0
+      IF( ANORM.EQ.ZERO )
+     $   GO TO 10
+      IF( ANORM.GT.SSFMAX ) THEN
+         ISCALE = 1
+         CALL PDLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
+     $                INFO )
+         CALL PDLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
+     $                INFO )
+      ELSE IF( ANORM.LT.SSFMIN ) THEN
+         ISCALE = 2
+         CALL PDLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
+     $                INFO )
+         CALL PDLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
+     $                INFO )
+      END IF
+*
+*     Choose between QL and QR iteration
+*
+      IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
+         LEND = LSV
+         L = LENDSV
+      END IF
+*
+      IF( LEND.GT.L ) THEN
+*
+*        QL Iteration
+*
+*        Look for small subdiagonal element.
+*
+   40    CONTINUE
+         IF( L.NE.LEND ) THEN
+            LENDM1 = LEND - 1
+            DO 50 M = L, LENDM1
+               TST = ABS( E( M ) )**2
+               IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M+1 ) )+
+     $             SAFMIN )GO TO 60
+   50       CONTINUE
+         END IF
+*
+         M = LEND
+*
+   60    CONTINUE
+         IF( M.LT.LEND )
+     $      E( M ) = ZERO
+         P = D( L )
+         IF( M.EQ.L )
+     $      GO TO 80
+*
+*        If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
+*        to compute its eigensystem.
+*
+         IF( M.EQ.L+1 ) THEN
+            IF( ICOMPZ.GT.0 ) THEN
+               CALL PDLAEV2( D( L ), E( L ), D( L+1 ), RT1, RT2, C, S )
+               WORK( L ) = C
+               WORK( N-1+L ) = S
+               CALL PDLASR( 'R', 'V', 'B', nz, 2, WORK( L ),
+     $                     WORK( N-1+L ), Z( 1, L ), nz )
+            ELSE
+               CALL PDLAE2( D( L ), E( L ), D( L+1 ), RT1, RT2 )
+            END IF
+            D( L ) = RT1
+            D( L+1 ) = RT2
+            E( L ) = ZERO
+            L = L + 2
+            IF( L.LE.LEND )
+     $         GO TO 40
+            GO TO 140
+         END IF
+*
+         IF( JTOT.EQ.NMAXIT )
+     $      GO TO 140
+         JTOT = JTOT + 1
+*
+*        Form shift.
+*
+         G = ( D( L+1 )-P ) / ( TWO*E( L ) )
+         R = PDLAPY2( G, ONE )
+         G = D( M ) - P + ( E( L ) / ( G+SIGN( R, G ) ) )
+*
+         S = ONE
+         C = ONE
+         P = ZERO
+*
+*        Inner loop
+*
+         MM1 = M - 1
+         DO 70 I = MM1, L, -1
+            F = S*E( I )
+            B = C*E( I )
+            CALL PDLARTG( G, F, C, S, R )
+            IF( I.NE.M-1 )
+     $         E( I+1 ) = R
+            G = D( I+1 ) - P
+            R = ( D( I )-G )*S + TWO*C*B
+            P = S*R
+            D( I+1 ) = G + P
+            G = C*R - B
+*
+*           If eigenvectors are desired, then save rotations.
+*
+            IF( ICOMPZ.GT.0 ) THEN
+               WORK( I ) = C
+               WORK( N-1+I ) = -S
+            END IF
+*
+   70    CONTINUE
+*
+*        If eigenvectors are desired, then apply saved rotations.
+*
+         IF( ICOMPZ.GT.0 ) THEN
+            MM = M - L + 1
+            CALL PDLASR('R', 'V', 'B', nz, MM, WORK( L ), WORK( N-1+L ),
+     $                  Z( 1, L ), nz )
+         END IF
+*
+         D( L ) = D( L ) - P
+         E( L ) = G
+         GO TO 40
+*
+*        Eigenvalue found.
+*
+   80    CONTINUE
+         D( L ) = P
+*
+         L = L + 1
+         IF( L.LE.LEND )
+     $      GO TO 40
+         GO TO 140
+*
+      ELSE
+*
+*        QR Iteration
+*
+*        Look for small superdiagonal element.
+*
+   90    CONTINUE
+         IF( L.NE.LEND ) THEN
+            LENDP1 = LEND + 1
+            DO 100 M = L, LENDP1, -1
+               TST = ABS( E( M-1 ) )**2
+               IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M-1 ) )+
+     $             SAFMIN )GO TO 110
+  100       CONTINUE
+         END IF
+*
+         M = LEND
+*
+  110    CONTINUE
+         IF( M.GT.LEND )
+     $      E( M-1 ) = ZERO
+         P = D( L )
+         IF( M.EQ.L )
+     $      GO TO 130
+*
+*        If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
+*        to compute its eigensystem.
+*
+         IF( M.EQ.L-1 ) THEN
+            IF( ICOMPZ.GT.0 ) THEN
+               CALL PDLAEV2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C,S)
+               WORK( M ) = C
+               WORK( N-1+M ) = S
+               CALL PDLASR( 'R', 'V', 'F', nz, 2, WORK( M ),
+     $                     WORK( N-1+M ), Z( 1, L-1 ), nz )
+            ELSE
+               CALL PDLAE2( D( L-1 ), E( L-1 ), D( L ), RT1, RT2 )
+            END IF
+            D( L-1 ) = RT1
+            D( L ) = RT2
+            E( L-1 ) = ZERO
+            L = L - 2
+            IF( L.GE.LEND )
+     $         GO TO 90
+            GO TO 140
+         END IF
+*
+         IF( JTOT.EQ.NMAXIT )
+     $      GO TO 140
+         JTOT = JTOT + 1
+*
+*        Form shift.
+*
+         G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) )
+         R = PDLAPY2( G, ONE )
+         G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) )
+*
+         S = ONE
+         C = ONE
+         P = ZERO
+*
+*        Inner loop
+*
+         LM1 = L - 1
+         DO 120 I = M, LM1
+            F = S*E( I )
+            B = C*E( I )
+            CALL PDLARTG( G, F, C, S, R )
+            IF( I.NE.M )
+     $         E( I-1 ) = R
+            G = D( I ) - P
+            R = ( D( I+1 )-G )*S + TWO*C*B
+            P = S*R
+            D( I ) = G + P
+            G = C*R - B
+*
+*           If eigenvectors are desired, then save rotations.
+*
+            IF( ICOMPZ.GT.0 ) THEN
+               WORK( I ) = C
+               WORK( N-1+I ) = S
+            END IF
+*
+  120    CONTINUE
+*
+*        If eigenvectors are desired, then apply saved rotations.
+*
+         IF( ICOMPZ.GT.0 ) THEN
+            MM = L - M + 1
+            CALL PDLASR('R', 'V', 'F', nz, MM, WORK( M ), WORK( N-1+M ),
+     $                  Z( 1, M ), nz )
+         END IF
+*
+         D( L ) = D( L ) - P
+         E( LM1 ) = G
+         GO TO 90
+*
+*        Eigenvalue found.
+*
+  130    CONTINUE
+         D( L ) = P
+*
+         L = L - 1
+         IF( L.GE.LEND )
+     $      GO TO 90
+         GO TO 140
+*
+      END IF
+*
+*     Undo scaling if necessary
+*
+  140 CONTINUE
+      IF( ISCALE.EQ.1 ) THEN
+         CALL PDLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
+     $                D( LSV ), N, INFO )
+         CALL PDLASCL('G',0, 0, SSFMAX, ANORM, LENDSV-LSV, 1, E( LSV ),
+     $                N, INFO )
+      ELSE IF( ISCALE.EQ.2 ) THEN
+         CALL PDLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
+     $                D( LSV ), N, INFO )
+         CALL PDLASCL ('G',0, 0, SSFMIN, ANORM, LENDSV-LSV, 1, E( LSV ),
+     $                N, INFO )
+      END IF
+*
+*     Check for no convergence to an eigenvalue after a total
+*     of N*MAXIT iterations.
+*
+      IF( JTOT.LT.NMAXIT )
+     $   GO TO 10
+      DO 150 I = 1, N - 1
+         IF( E( I ).NE.ZERO )
+     $      INFO = INFO + 1
+  150 CONTINUE
+      GO TO 190
+*
+*     Order eigenvalues and eigenvectors.
+*
+  160 CONTINUE
+      IF( ICOMPZ.EQ.0 ) THEN
+*
+*        Use Quick Sort
+*
+         CALL PDLASRT( 'I', N, D, INFO )
+*
+      ELSE
+*
+*        Use Selection Sort to minimize swaps of eigenvectors
+*
+         DO 180 II = 2, N
+            I = II - 1
+            K = I
+            P = D( I )
+            DO 170 J = II, N
+               IF( D( J ).LT.P ) THEN
+                  K = J
+                  P = D( J )
+               END IF
+  170       CONTINUE
+            IF( K.NE.I ) THEN
+               D( K ) = D( I )
+               D( I ) = P
+               CALL DSWAP( nz, Z( 1, I ), 1, Z( 1, K ), 1 )
+            END IF
+  180    CONTINUE
+      END IF
+*
+  190 CONTINUE
+      RETURN
+*
+*     End of DSTEQR
+*
+      END
+
+* Auxiliary routines are not provided with all commerical lapacks,
+* so a local copy for use by PDSTEQR is provided here.
+      SUBROUTINE PDLAE2( A, B, C, RT1, RT2 )
+*
+*  -- LAPACK auxiliary routine (version 3.0) --
+*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+*     Courant Institute, Argonne National Lab, and Rice University
+*     October 31, 1992
+*
+*     .. Scalar Arguments ..
+      DOUBLE PRECISION   A, B, C, RT1, RT2
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DLAE2  computes the eigenvalues of a 2-by-2 symmetric matrix
+*     [  A   B  ]
+*     [  B   C  ].
+*  On return, RT1 is the eigenvalue of larger absolute value, and RT2
+*  is the eigenvalue of smaller absolute value.
+*
+*  Arguments
+*  =========
+*
+*  A       (input) DOUBLE PRECISION
+*          The (1,1) element of the 2-by-2 matrix.
+*
+*  B       (input) DOUBLE PRECISION
+*          The (1,2) and (2,1) elements of the 2-by-2 matrix.
+*
+*  C       (input) DOUBLE PRECISION
+*          The (2,2) element of the 2-by-2 matrix.
+*
+*  RT1     (output) DOUBLE PRECISION
+*          The eigenvalue of larger absolute value.
+*
+*  RT2     (output) DOUBLE PRECISION
+*          The eigenvalue of smaller absolute value.
+*
+*  Further Details
+*  ===============
+*
+*  RT1 is accurate to a few ulps barring over/underflow.
+*
+*  RT2 may be inaccurate if there is massive cancellation in the
+*  determinant A*C-B*B; higher precision or correctly rounded or
+*  correctly truncated arithmetic would be needed to compute RT2
+*  accurately in all cases.
+*
+*  Overflow is possible only if RT1 is within a factor of 5 of overflow.
+*  Underflow is harmless if the input data is 0 or exceeds
+*     underflow_threshold / macheps.
+*
+* =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ONE
+      PARAMETER          ( ONE = 1.0D0 )
+      DOUBLE PRECISION   TWO
+      PARAMETER          ( TWO = 2.0D0 )
+      DOUBLE PRECISION   ZERO
+      PARAMETER          ( ZERO = 0.0D0 )
+      DOUBLE PRECISION   HALF
+      PARAMETER          ( HALF = 0.5D0 )
+*     ..
+*     .. Local Scalars ..
+      DOUBLE PRECISION   AB, ACMN, ACMX, ADF, DF, RT, SM, TB
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, SQRT
+*     ..
+*     .. Executable Statements ..
+*
+*     Compute the eigenvalues
+*
+      SM = A + C
+      DF = A - C
+      ADF = ABS( DF )
+      TB = B + B
+      AB = ABS( TB )
+      IF( ABS( A ).GT.ABS( C ) ) THEN
+         ACMX = A
+         ACMN = C
+      ELSE
+         ACMX = C
+         ACMN = A
+      END IF
+      IF( ADF.GT.AB ) THEN
+         RT = ADF*SQRT( ONE+( AB / ADF )**2 )
+      ELSE IF( ADF.LT.AB ) THEN
+         RT = AB*SQRT( ONE+( ADF / AB )**2 )
+      ELSE
+*
+*        Includes case AB=ADF=0
+*
+         RT = AB*SQRT( TWO )
+      END IF
+      IF( SM.LT.ZERO ) THEN
+         RT1 = HALF*( SM-RT )
+*
+*        Order of execution important.
+*        To get fully accurate smaller eigenvalue,
+*        next line needs to be executed in higher precision.
+*
+         RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
+      ELSE IF( SM.GT.ZERO ) THEN
+         RT1 = HALF*( SM+RT )
+*
+*        Order of execution important.
+*        To get fully accurate smaller eigenvalue,
+*        next line needs to be executed in higher precision.
+*
+         RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
+      ELSE
+*
+*        Includes case RT1 = RT2 = 0
+*
+         RT1 = HALF*RT
+         RT2 = -HALF*RT
+      END IF
+      RETURN
+*
+*     End of DLAE2
+*
+      END
+      SUBROUTINE PDLAEV2( A, B, C, RT1, RT2, CS1, SN1 )
+*
+*  -- LAPACK auxiliary routine (version 3.0) --
+*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+*     Courant Institute, Argonne National Lab, and Rice University
+*     October 31, 1992
+*
+*     .. Scalar Arguments ..
+      DOUBLE PRECISION   A, B, C, CS1, RT1, RT2, SN1
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix
+*     [  A   B  ]
+*     [  B   C  ].
+*  On return, RT1 is the eigenvalue of larger absolute value, RT2 is the
+*  eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right
+*  eigenvector for RT1, giving the decomposition
+*
+*     [ CS1  SN1 ] [  A   B  ] [ CS1 -SN1 ]  =  [ RT1  0  ]
+*     [-SN1  CS1 ] [  B   C  ] [ SN1  CS1 ]     [  0  RT2 ].
+*
+*  Arguments
+*  =========
+*
+*  A       (input) DOUBLE PRECISION
+*          The (1,1) element of the 2-by-2 matrix.
+*
+*  B       (input) DOUBLE PRECISION
+*          The (1,2) element and the conjugate of the (2,1) element of
+*          the 2-by-2 matrix.
+*
+*  C       (input) DOUBLE PRECISION
+*          The (2,2) element of the 2-by-2 matrix.
+*
+*  RT1     (output) DOUBLE PRECISION
+*          The eigenvalue of larger absolute value.
+*
+*  RT2     (output) DOUBLE PRECISION
+*          The eigenvalue of smaller absolute value.
+*
+*  CS1     (output) DOUBLE PRECISION
+*  SN1     (output) DOUBLE PRECISION
+*          The vector (CS1, SN1) is a unit right eigenvector for RT1.
+*
+*  Further Details
+*  ===============
+*
+*  RT1 is accurate to a few ulps barring over/underflow.
+*
+*  RT2 may be inaccurate if there is massive cancellation in the
+*  determinant A*C-B*B; higher precision or correctly rounded or
+*  correctly truncated arithmetic would be needed to compute RT2
+*  accurately in all cases.
+*
+*  CS1 and SN1 are accurate to a few ulps barring over/underflow.
+*
+*  Overflow is possible only if RT1 is within a factor of 5 of overflow.
+*  Underflow is harmless if the input data is 0 or exceeds
+*     underflow_threshold / macheps.
+*
+* =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ONE
+      PARAMETER          ( ONE = 1.0D0 )
+      DOUBLE PRECISION   TWO
+      PARAMETER          ( TWO = 2.0D0 )
+      DOUBLE PRECISION   ZERO
+      PARAMETER          ( ZERO = 0.0D0 )
+      DOUBLE PRECISION   HALF
+      PARAMETER          ( HALF = 0.5D0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            SGN1, SGN2
+      DOUBLE PRECISION   AB, ACMN, ACMX, ACS, ADF, CS, CT, DF, RT, SM,
+     $                   TB, TN
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, SQRT
+*     ..
+*     .. Executable Statements ..
+*
+*     Compute the eigenvalues
+*
+      SM = A + C
+      DF = A - C
+      ADF = ABS( DF )
+      TB = B + B
+      AB = ABS( TB )
+      IF( ABS( A ).GT.ABS( C ) ) THEN
+         ACMX = A
+         ACMN = C
+      ELSE
+         ACMX = C
+         ACMN = A
+      END IF
+      IF( ADF.GT.AB ) THEN
+         RT = ADF*SQRT( ONE+( AB / ADF )**2 )
+      ELSE IF( ADF.LT.AB ) THEN
+         RT = AB*SQRT( ONE+( ADF / AB )**2 )
+      ELSE
+*
+*        Includes case AB=ADF=0
+*
+         RT = AB*SQRT( TWO )
+      END IF
+      IF( SM.LT.ZERO ) THEN
+         RT1 = HALF*( SM-RT )
+         SGN1 = -1
+*
+*        Order of execution important.
+*        To get fully accurate smaller eigenvalue,
+*        next line needs to be executed in higher precision.
+*
+         RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
+      ELSE IF( SM.GT.ZERO ) THEN
+         RT1 = HALF*( SM+RT )
+         SGN1 = 1
+*
+*        Order of execution important.
+*        To get fully accurate smaller eigenvalue,
+*        next line needs to be executed in higher precision.
+*
+         RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
+      ELSE
+*
+*        Includes case RT1 = RT2 = 0
+*
+         RT1 = HALF*RT
+         RT2 = -HALF*RT
+         SGN1 = 1
+      END IF
+*
+*     Compute the eigenvector
+*
+      IF( DF.GE.ZERO ) THEN
+         CS = DF + RT
+         SGN2 = 1
+      ELSE
+         CS = DF - RT
+         SGN2 = -1
+      END IF
+      ACS = ABS( CS )
+      IF( ACS.GT.AB ) THEN
+         CT = -TB / CS
+         SN1 = ONE / SQRT( ONE+CT*CT )
+         CS1 = CT*SN1
+      ELSE
+         IF( AB.EQ.ZERO ) THEN
+            CS1 = ONE
+            SN1 = ZERO
+         ELSE
+            TN = -CS / TB
+            CS1 = ONE / SQRT( ONE+TN*TN )
+            SN1 = TN*CS1
+         END IF
+      END IF
+      IF( SGN1.EQ.SGN2 ) THEN
+         TN = CS1
+         CS1 = -SN1
+         SN1 = TN
+      END IF
+      RETURN
+*
+*     End of DLAEV2
+*
+      END
+      DOUBLE PRECISION FUNCTION PDLAMCH( CMACH )
+*
+*  -- LAPACK auxiliary routine (version 3.0) --
+*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+*     Courant Institute, Argonne National Lab, and Rice University
+*     October 31, 1992
+*
+*     .. Scalar Arguments ..
+      CHARACTER          CMACH
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DLAMCH determines double precision machine parameters.
+*
+*  Arguments
+*  =========
+*
+*  CMACH   (input) CHARACTER*1
+*          Specifies the value to be returned by DLAMCH:
+*          = 'E' or 'e',   DLAMCH := eps
+*          = 'S' or 's ,   DLAMCH := sfmin
+*          = 'B' or 'b',   DLAMCH := base
+*          = 'P' or 'p',   DLAMCH := eps*base
+*          = 'N' or 'n',   DLAMCH := t
+*          = 'R' or 'r',   DLAMCH := rnd
+*          = 'M' or 'm',   DLAMCH := emin
+*          = 'U' or 'u',   DLAMCH := rmin
+*          = 'L' or 'l',   DLAMCH := emax
+*          = 'O' or 'o',   DLAMCH := rmax
+*
+*          where
+*
+*          eps   = relative machine precision
+*          sfmin = safe minimum, such that 1/sfmin does not overflow
+*          base  = base of the machine
+*          prec  = eps*base
+*          t     = number of (base) digits in the mantissa
+*          rnd   = 1.0 when rounding occurs in addition, 0.0 otherwise
+*          emin  = minimum exponent before (gradual) underflow
+*          rmin  = underflow threshold - base**(emin-1)
+*          emax  = largest exponent before overflow
+*          rmax  = overflow threshold  - (base**emax)*(1-eps)
+*
+* =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ONE, ZERO
+      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            FIRST, LRND
+      INTEGER            BETA, IMAX, IMIN, IT
+      DOUBLE PRECISION   BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN,
+     $                   RND, SFMIN, SMALL, T
+*     ..
+*     .. External Functions ..
+      LOGICAL            PLSAME
+      EXTERNAL           PLSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           PDLAMC2
+*     ..
+*     .. Save statement ..
+      SAVE               FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN,
+     $                   EMAX, RMAX, PREC
+*     ..
+*     .. Data statements ..
+      DATA               FIRST / .TRUE. /
+*     ..
+*     .. Executable Statements ..
+*
+      IF( FIRST ) THEN
+         FIRST = .FALSE.
+         CALL PDLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX )
+         BASE = BETA
+         T = IT
+         IF( LRND ) THEN
+            RND = ONE
+            EPS = ( BASE**( 1-IT ) ) / 2
+         ELSE
+            RND = ZERO
+            EPS = BASE**( 1-IT )
+         END IF
+         PREC = EPS*BASE
+         EMIN = IMIN
+         EMAX = IMAX
+         SFMIN = RMIN
+         SMALL = ONE / RMAX
+         IF( SMALL.GE.SFMIN ) THEN
+*
+*           Use SMALL plus a bit, to avoid the possibility of rounding
+*           causing overflow when computing  1/sfmin.
+*
+            SFMIN = SMALL*( ONE+EPS )
+         END IF
+      END IF
+*
+      IF( PLSAME( CMACH, 'E' ) ) THEN
+         RMACH = EPS
+      ELSE IF( PLSAME( CMACH, 'S' ) ) THEN
+         RMACH = SFMIN
+      ELSE IF( PLSAME( CMACH, 'B' ) ) THEN
+         RMACH = BASE
+      ELSE IF( PLSAME( CMACH, 'P' ) ) THEN
+         RMACH = PREC
+      ELSE IF( PLSAME( CMACH, 'N' ) ) THEN
+         RMACH = T
+      ELSE IF( PLSAME( CMACH, 'R' ) ) THEN
+         RMACH = RND
+      ELSE IF( PLSAME( CMACH, 'M' ) ) THEN
+         RMACH = EMIN
+      ELSE IF( PLSAME( CMACH, 'U' ) ) THEN
+         RMACH = RMIN
+      ELSE IF( PLSAME( CMACH, 'L' ) ) THEN
+         RMACH = EMAX
+      ELSE IF( PLSAME( CMACH, 'O' ) ) THEN
+         RMACH = RMAX
+      END IF
+*
+      PDLAMCH = RMACH
+      RETURN
+*
+*     End of DLAMCH
+*
+      END
+*
+************************************************************************
+*
+      SUBROUTINE PDLAMC1( BETA, T, RND, IEEE1 )
+*
+*  -- LAPACK auxiliary routine (version 3.0) --
+*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+*     Courant Institute, Argonne National Lab, and Rice University
+*     October 31, 1992
+*
+*     .. Scalar Arguments ..
+      LOGICAL            IEEE1, RND
+      INTEGER            BETA, T
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DLAMC1 determines the machine parameters given by BETA, T, RND, and
+*  IEEE1.
+*
+*  Arguments
+*  =========
+*
+*  BETA    (output) INTEGER
+*          The base of the machine.
+*
+*  T       (output) INTEGER
+*          The number of ( BETA ) digits in the mantissa.
+*
+*  RND     (output) LOGICAL
+*          Specifies whether proper rounding  ( RND = .TRUE. )  or
+*          chopping  ( RND = .FALSE. )  occurs in addition. This may not
+*          be a reliable guide to the way in which the machine performs
+*          its arithmetic.
+*
+*  IEEE1   (output) LOGICAL
+*          Specifies whether rounding appears to be done in the IEEE
+*          'round to nearest' style.
+*
+*  Further Details
+*  ===============
+*
+*  The routine is based on the routine  ENVRON  by Malcolm and
+*  incorporates suggestions by Gentleman and Marovich. See
+*
+*     Malcolm M. A. (1972) Algorithms to reveal properties of
+*        floating-point arithmetic. Comms. of the ACM, 15, 949-951.
+*
+*     Gentleman W. M. and Marovich S. B. (1974) More on algorithms
+*        that reveal properties of floating point arithmetic units.
+*        Comms. of the ACM, 17, 276-277.
+*
+* =====================================================================
+*
+*     .. Local Scalars ..
+      LOGICAL            FIRST, LIEEE1, LRND
+      INTEGER            LBETA, LT
+      DOUBLE PRECISION   A, B, C, F, ONE, QTR, SAVEC, T1, T2
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   PDLAMC3
+      EXTERNAL           PDLAMC3
+*     ..
+*     .. Save statement ..
+      SAVE               FIRST, LIEEE1, LBETA, LRND, LT
+*     ..
+*     .. Data statements ..
+      DATA               FIRST / .TRUE. /
+*     ..
+*     .. Executable Statements ..
+*
+      IF( FIRST ) THEN
+         FIRST = .FALSE.
+         ONE = 1
+*
+*        LBETA,  LIEEE1,  LT and  LRND  are the  local values  of  BETA,
+*        IEEE1, T and RND.
+*
+*        Throughout this routine  we use the function  DLAMC3  to ensure
+*        that relevant values are  stored and not held in registers,  or
+*        are not affected by optimizers.
+*
+*        Compute  a = 2.0**m  with the  smallest positive integer m such
+*        that
+*
+*           fl( a + 1.0 ) = a.
+*
+         A = 1
+         C = 1
+*
+*+       WHILE( C.EQ.ONE )LOOP
+   10    CONTINUE
+         IF( C.EQ.ONE ) THEN
+            A = 2*A
+            C = PDLAMC3( A, ONE )
+            C = PDLAMC3( C, -A )
+            GO TO 10
+         END IF
+*+       END WHILE
+*
+*        Now compute  b = 2.0**m  with the smallest positive integer m
+*        such that
+*
+*           fl( a + b ) .gt. a.
+*
+         B = 1
+         C = PDLAMC3( A, B )
+*
+*+       WHILE( C.EQ.A )LOOP
+   20    CONTINUE
+         IF( C.EQ.A ) THEN
+            B = 2*B
+            C = PDLAMC3( A, B )
+            GO TO 20
+         END IF
+*+       END WHILE
+*
+*        Now compute the base.  a and c  are neighbouring floating point
+*        numbers  in the  interval  ( beta**t, beta**( t + 1 ) )  and so
+*        their difference is beta. Adding 0.25 to c is to ensure that it
+*        is truncated to beta and not ( beta - 1 ).
+*
+         QTR = ONE / 4
+         SAVEC = C
+         C = PDLAMC3( C, -A )
+         LBETA = C + QTR
+*
+*        Now determine whether rounding or chopping occurs,  by adding a
+*        bit  less  than  beta/2  and a  bit  more  than  beta/2  to  a.
+*
+         B = LBETA
+         F = PDLAMC3( B / 2, -B / 100 )
+         C = PDLAMC3( F, A )
+         IF( C.EQ.A ) THEN
+            LRND = .TRUE.
+         ELSE
+            LRND = .FALSE.
+         END IF
+         F = PDLAMC3( B / 2, B / 100 )
+         C = PDLAMC3( F, A )
+         IF( ( LRND ) .AND. ( C.EQ.A ) )
+     $      LRND = .FALSE.
+*
+*        Try and decide whether rounding is done in the  IEEE  'round to
+*        nearest' style. B/2 is half a unit in the last place of the two
+*        numbers A and SAVEC. Furthermore, A is even, i.e. has last  bit
+*        zero, and SAVEC is odd. Thus adding B/2 to A should not  change
+*        A, but adding B/2 to SAVEC should change SAVEC.
+*
+         T1 = PDLAMC3( B / 2, A )
+         T2 = PDLAMC3( B / 2, SAVEC )
+         LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND
+*
+*        Now find  the  mantissa, t.  It should  be the  integer part of
+*        log to the base beta of a,  however it is safer to determine  t
+*        by powering.  So we find t as the smallest positive integer for
+*        which
+*
+*           fl( beta**t + 1.0 ) = 1.0.
+*
+         LT = 0
+         A = 1
+         C = 1
+*
+*+       WHILE( C.EQ.ONE )LOOP
+   30    CONTINUE
+         IF( C.EQ.ONE ) THEN
+            LT = LT + 1
+            A = A*LBETA
+            C = PDLAMC3( A, ONE )
+            C = PDLAMC3( C, -A )
+            GO TO 30
+         END IF
+*+       END WHILE
+*
+      END IF
+*
+      BETA = LBETA
+      T = LT
+      RND = LRND
+      IEEE1 = LIEEE1
+      RETURN
+*
+*     End of DLAMC1
+*
+      END
+*
+************************************************************************
+*
+      SUBROUTINE PDLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
+*
+*  -- LAPACK auxiliary routine (version 3.0) --
+*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+*     Courant Institute, Argonne National Lab, and Rice University
+*     October 31, 1992
+*
+*     .. Scalar Arguments ..
+      LOGICAL            RND
+      INTEGER            BETA, EMAX, EMIN, T
+      DOUBLE PRECISION   EPS, RMAX, RMIN
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DLAMC2 determines the machine parameters specified in its argument
+*  list.
+*
+*  Arguments
+*  =========
+*
+*  BETA    (output) INTEGER
+*          The base of the machine.
+*
+*  T       (output) INTEGER
+*          The number of ( BETA ) digits in the mantissa.
+*
+*  RND     (output) LOGICAL
+*          Specifies whether proper rounding  ( RND = .TRUE. )  or
+*          chopping  ( RND = .FALSE. )  occurs in addition. This may not
+*          be a reliable guide to the way in which the machine performs
+*          its arithmetic.
+*
+*  EPS     (output) DOUBLE PRECISION
+*          The smallest positive number such that
+*
+*             fl( 1.0 - EPS ) .LT. 1.0,
+*
+*          where fl denotes the computed value.
+*
+*  EMIN    (output) INTEGER
+*          The minimum exponent before (gradual) underflow occurs.
+*
+*  RMIN    (output) DOUBLE PRECISION
+*          The smallest normalized number for the machine, given by
+*          BASE**( EMIN - 1 ), where  BASE  is the floating point value
+*          of BETA.
+*
+*  EMAX    (output) INTEGER
+*          The maximum exponent before overflow occurs.
+*
+*  RMAX    (output) DOUBLE PRECISION
+*          The largest positive number for the machine, given by
+*          BASE**EMAX * ( 1 - EPS ), where  BASE  is the floating point
+*          value of BETA.
+*
+*  Further Details
+*  ===============
+*
+*  The computation of  EPS  is based on a routine PARANOIA by
+*  W. Kahan of the University of California at Berkeley.
+*
+* =====================================================================
+*
+*     .. Local Scalars ..
+      LOGICAL            FIRST, IEEE, IWARN, LIEEE1, LRND
+      INTEGER            GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT,
+     $                   NGNMIN, NGPMIN
+      DOUBLE PRECISION   A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
+     $                   SIXTH, SMALL, THIRD, TWO, ZERO
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   PDLAMC3
+      EXTERNAL           PDLAMC3
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           PDLAMC1, PDLAMC4, PDLAMC5
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, MAX, MIN
+*     ..
+*     .. Save statement ..
+      SAVE               FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX,
+     $                   LRMIN, LT
+*     ..
+*     .. Data statements ..
+      DATA               FIRST / .TRUE. / , IWARN / .FALSE. /
+*     ..
+*     .. Executable Statements ..
+*
+      IF( FIRST ) THEN
+         FIRST = .FALSE.
+         ZERO = 0
+         ONE = 1
+         TWO = 2
+*
+*        LBETA, LT, LRND, LEPS, LEMIN and LRMIN  are the local values of
+*        BETA, T, RND, EPS, EMIN and RMIN.
+*
+*        Throughout this routine  we use the function  DLAMC3  to ensure
+*        that relevant values are stored  and not held in registers,  or
+*        are not affected by optimizers.
+*
+*        DLAMC1 returns the parameters  LBETA, LT, LRND and LIEEE1.
+*
+         CALL PDLAMC1( LBETA, LT, LRND, LIEEE1 )
+*
+*        Start to find EPS.
+*
+         B = LBETA
+         A = B**( -LT )
+         LEPS = A
+*
+*        Try some tricks to see whether or not this is the correct  EPS.
+*
+         B = TWO / 3
+         HALF = ONE / 2
+         SIXTH = PDLAMC3( B, -HALF )
+         THIRD = PDLAMC3( SIXTH, SIXTH )
+         B = PDLAMC3( THIRD, -HALF )
+         B = PDLAMC3( B, SIXTH )
+         B = ABS( B )
+         IF( B.LT.LEPS )
+     $      B = LEPS
+*
+         LEPS = 1
+*
+*+       WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
+   10    CONTINUE
+         IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN
+            LEPS = B
+            C = PDLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) )
+            C = PDLAMC3( HALF, -C )
+            B = PDLAMC3( HALF, C )
+            C = PDLAMC3( HALF, -B )
+            B = PDLAMC3( HALF, C )
+            GO TO 10
+         END IF
+*+       END WHILE
+*
+         IF( A.LT.LEPS )
+     $      LEPS = A
+*
+*        Computation of EPS complete.
+*
+*        Now find  EMIN.  Let A = + or - 1, and + or - (1 + BASE**(-3)).
+*        Keep dividing  A by BETA until (gradual) underflow occurs. This
+*        is detected when we cannot recover the previous A.
+*
+         RBASE = ONE / LBETA
+         SMALL = ONE
+         DO 20 I = 1, 3
+            SMALL = PDLAMC3( SMALL*RBASE, ZERO )
+   20    CONTINUE
+         A = PDLAMC3( ONE, SMALL )
+         CALL PDLAMC4( NGPMIN, ONE, LBETA )
+         CALL PDLAMC4( NGNMIN, -ONE, LBETA )
+         CALL PDLAMC4( GPMIN, A, LBETA )
+         CALL PDLAMC4( GNMIN, -A, LBETA )
+         IEEE = .FALSE.
+*
+         IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN
+            IF( NGPMIN.EQ.GPMIN ) THEN
+               LEMIN = NGPMIN
+*            ( Non twos-complement machines, no gradual underflow;
+*              e.g.,  VAX )
+            ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN
+               LEMIN = NGPMIN - 1 + LT
+               IEEE = .TRUE.
+*            ( Non twos-complement machines, with gradual underflow;
+*              e.g., IEEE standard followers )
+            ELSE
+               LEMIN = MIN( NGPMIN, GPMIN )
+*            ( A guess; no known machine )
+               IWARN = .TRUE.
+            END IF
+*
+         ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN
+            IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN
+               LEMIN = MAX( NGPMIN, NGNMIN )
+*            ( Twos-complement machines, no gradual underflow;
+*              e.g., CYBER 205 )
+            ELSE
+               LEMIN = MIN( NGPMIN, NGNMIN )
+*            ( A guess; no known machine )
+               IWARN = .TRUE.
+            END IF
+*
+         ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND.
+     $            ( GPMIN.EQ.GNMIN ) ) THEN
+            IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN
+               LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT
+*            ( Twos-complement machines with gradual underflow;
+*              no known machine )
+            ELSE
+               LEMIN = MIN( NGPMIN, NGNMIN )
+*            ( A guess; no known machine )
+               IWARN = .TRUE.
+            END IF
+*
+         ELSE
+            LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN )
+*         ( A guess; no known machine )
+            IWARN = .TRUE.
+         END IF
+***
+* Comment out this if block if EMIN is ok
+         IF( IWARN ) THEN
+            FIRST = .TRUE.
+            WRITE( 6, FMT = 9999 )LEMIN
+         END IF
+***
+*
+*        Assume IEEE arithmetic if we found denormalised  numbers above,
+*        or if arithmetic seems to round in the  IEEE style,  determined
+*        in routine DLAMC1. A true IEEE machine should have both  things
+*        true; however, faulty machines may have one or the other.
+*
+         IEEE = IEEE .OR. LIEEE1
+*
+*        Compute  RMIN by successive division by  BETA. We could compute
+*        RMIN as BASE**( EMIN - 1 ),  but some machines underflow during
+*        this computation.
+*
+         LRMIN = 1
+         DO 30 I = 1, 1 - LEMIN
+            LRMIN = PDLAMC3( LRMIN*RBASE, ZERO )
+   30    CONTINUE
+*
+*        Finally, call DLAMC5 to compute EMAX and RMAX.
+*
+         CALL PDLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX )
+      END IF
+*
+      BETA = LBETA
+      T = LT
+      RND = LRND
+      EPS = LEPS
+      EMIN = LEMIN
+      RMIN = LRMIN
+      EMAX = LEMAX
+      RMAX = LRMAX
+*
+      RETURN
+*
+ 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',
+     $      '  EMIN = ', I8, /
+     $      ' If, after inspection, the value EMIN looks',
+     $      ' acceptable please comment out ',
+     $      / ' the IF block as marked within the code of routine',
+     $      ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / )
+*
+*     End of DLAMC2
+*
+      END
+*
+************************************************************************
+*
+      DOUBLE PRECISION FUNCTION PDLAMC3( A, B )
+*
+*  -- LAPACK auxiliary routine (version 3.0) --
+*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+*     Courant Institute, Argonne National Lab, and Rice University
+*     October 31, 1992
+*
+*     .. Scalar Arguments ..
+      DOUBLE PRECISION   A, B
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DLAMC3  is intended to force  A  and  B  to be stored prior to doing
+*  the addition of  A  and  B ,  for use in situations where optimizers
+*  might hold one of these in a register.
+*
+*  Arguments
+*  =========
+*
+*  A, B    (input) DOUBLE PRECISION
+*          The values A and B.
+*
+* =====================================================================
+*
+*     .. Executable Statements ..
+*
+      PDLAMC3 = A + B
+*
+      RETURN
+*
+*     End of DLAMC3
+*
+      END
+*
+************************************************************************
+*
+      SUBROUTINE PDLAMC4( EMIN, START, BASE )
+*
+*  -- LAPACK auxiliary routine (version 3.0) --
+*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+*     Courant Institute, Argonne National Lab, and Rice University
+*     October 31, 1992
+*
+*     .. Scalar Arguments ..
+      INTEGER            BASE, EMIN
+      DOUBLE PRECISION   START
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DLAMC4 is a service routine for DLAMC2.
+*
+*  Arguments
+*  =========
+*
+*  EMIN    (output) EMIN
+*          The minimum exponent before (gradual) underflow, computed by
+*          setting A = START and dividing by BASE until the previous A
+*          can not be recovered.
+*
+*  START   (input) DOUBLE PRECISION
+*          The starting point for determining EMIN.
+*
+*  BASE    (input) INTEGER
+*          The base of the machine.
+*
+* =====================================================================
+*
+*     .. Local Scalars ..
+      INTEGER            I
+      DOUBLE PRECISION   A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   PDLAMC3
+      EXTERNAL           PDLAMC3
+*     ..
+*     .. Executable Statements ..
+*
+      A = START
+      ONE = 1
+      RBASE = ONE / BASE
+      ZERO = 0
+      EMIN = 1
+      B1 = PDLAMC3( A*RBASE, ZERO )
+      C1 = A
+      C2 = A
+      D1 = A
+      D2 = A
+*+    WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
+*    $       ( D1.EQ.A ).AND.( D2.EQ.A )      )LOOP
+   10 CONTINUE
+      IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND.
+     $    ( D2.EQ.A ) ) THEN
+         EMIN = EMIN - 1
+         A = B1
+         B1 = PDLAMC3( A / BASE, ZERO )
+         C1 = PDLAMC3( B1*BASE, ZERO )
+         D1 = ZERO
+         DO 20 I = 1, BASE
+            D1 = D1 + B1
+   20    CONTINUE
+         B2 = PDLAMC3( A*RBASE, ZERO )
+         C2 = PDLAMC3( B2 / RBASE, ZERO )
+         D2 = ZERO
+         DO 30 I = 1, BASE
+            D2 = D2 + B2
+   30    CONTINUE
+         GO TO 10
+      END IF
+*+    END WHILE
+*
+      RETURN
+*
+*     End of DLAMC4
+*
+      END
+*
+************************************************************************
+*
+      SUBROUTINE PDLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX )
+*
+*  -- LAPACK auxiliary routine (version 3.0) --
+*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+*     Courant Institute, Argonne National Lab, and Rice University
+*     October 31, 1992
+*
+*     .. Scalar Arguments ..
+      LOGICAL            IEEE
+      INTEGER            BETA, EMAX, EMIN, P
+      DOUBLE PRECISION   RMAX
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DLAMC5 attempts to compute RMAX, the largest machine floating-point
+*  number, without overflow.  It assumes that EMAX + abs(EMIN) sum
+*  approximately to a power of 2.  It will fail on machines where this
+*  assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
+*  EMAX = 28718).  It will also fail if the value supplied for EMIN is
+*  too large (i.e. too close to zero), probably with overflow.
+*
+*  Arguments
+*  =========
+*
+*  BETA    (input) INTEGER
+*          The base of floating-point arithmetic.
+*
+*  P       (input) INTEGER
+*          The number of base BETA digits in the mantissa of a
+*          floating-point value.
+*
+*  EMIN    (input) INTEGER
+*          The minimum exponent before (gradual) underflow.
+*
+*  IEEE    (input) LOGICAL
+*          A logical flag specifying whether or not the arithmetic
+*          system is thought to comply with the IEEE standard.
+*
+*  EMAX    (output) INTEGER
+*          The largest exponent before overflow
+*
+*  RMAX    (output) DOUBLE PRECISION
+*          The largest machine floating-point number.
+*
+* =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
+      DOUBLE PRECISION   OLDY, RECBAS, Y, Z
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   PDLAMC3
+      EXTERNAL           PDLAMC3
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MOD
+*     ..
+*     .. Executable Statements ..
+*
+*     First compute LEXP and UEXP, two powers of 2 that bound
+*     abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
+*     approximately to the bound that is closest to abs(EMIN).
+*     (EMAX is the exponent of the required number RMAX).
+*
+      LEXP = 1
+      EXBITS = 1
+   10 CONTINUE
+      TRY = LEXP*2
+      IF( TRY.LE.( -EMIN ) ) THEN
+         LEXP = TRY
+         EXBITS = EXBITS + 1
+         GO TO 10
+      END IF
+      IF( LEXP.EQ.-EMIN ) THEN
+         UEXP = LEXP
+      ELSE
+         UEXP = TRY
+         EXBITS = EXBITS + 1
+      END IF
+*
+*     Now -LEXP is less than or equal to EMIN, and -UEXP is greater
+*     than or equal to EMIN. EXBITS is the number of bits needed to
+*     store the exponent.
+*
+      IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN
+         EXPSUM = 2*LEXP
+      ELSE
+         EXPSUM = 2*UEXP
+      END IF
+*
+*     EXPSUM is the exponent range, approximately equal to
+*     EMAX - EMIN + 1 .
+*
+      EMAX = EXPSUM + EMIN - 1
+      NBITS = 1 + EXBITS + P
+*
+*     NBITS is the total number of bits needed to store a
+*     floating-point number.
+*
+      IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN
+*
+*        Either there are an odd number of bits used to store a
+*        floating-point number, which is unlikely, or some bits are
+*        not used in the representation of numbers, which is possible,
+*        (e.g. Cray machines) or the mantissa has an implicit bit,
+*        (e.g. IEEE machines, Dec Vax machines), which is perhaps the
+*        most likely. We have to assume the last alternative.
+*        If this is true, then we need to reduce EMAX by one because
+*        there must be some way of representing zero in an implicit-bit
+*        system. On machines like Cray, we are reducing EMAX by one
+*        unnecessarily.
+*
+         EMAX = EMAX - 1
+      END IF
+*
+      IF( IEEE ) THEN
+*
+*        Assume we are on an IEEE machine which reserves one exponent
+*        for infinity and NaN.
+*
+         EMAX = EMAX - 1
+      END IF
+*
+*     Now create RMAX, the largest machine number, which should
+*     be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
+*
+*     First compute 1.0 - BETA**(-P), being careful that the
+*     result is less than 1.0 .
+*
+      RECBAS = ONE / BETA
+      Z = BETA - ONE
+      Y = ZERO
+      DO 20 I = 1, P
+         Z = Z*RECBAS
+         IF( Y.LT.ONE )
+     $      OLDY = Y
+         Y = PDLAMC3( Y, Z )
+   20 CONTINUE
+      IF( Y.GE.ONE )
+     $   Y = OLDY
+*
+*     Now multiply by BETA**EMAX to get RMAX.
+*
+      DO 30 I = 1, EMAX
+         Y = PDLAMC3( Y*BETA, ZERO )
+   30 CONTINUE
+*
+      RMAX = Y
+      RETURN
+*
+*     End of DLAMC5
+*
+      END
+      DOUBLE PRECISION FUNCTION PDLANST( NORM, N, D, E )
+*
+*  -- LAPACK auxiliary routine (version 3.0) --
+*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+*     Courant Institute, Argonne National Lab, and Rice University
+*     February 29, 1992
+*
+*     .. Scalar Arguments ..
+      CHARACTER          NORM
+      INTEGER            N
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   D( * ), E( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DLANST  returns the value of the one norm,  or the Frobenius norm, or
+*  the  infinity norm,  or the  element of  largest absolute value  of a
+*  real symmetric tridiagonal matrix A.
+*
+*  Description
+*  ===========
+*
+*  DLANST returns the value
+*
+*     DLANST = ( max(abs(A(i,j))), NORM = 'M' or 'm'
+*              (
+*              ( norm1(A),         NORM = '1', 'O' or 'o'
+*              (
+*              ( normI(A),         NORM = 'I' or 'i'
+*              (
+*              ( normF(A),         NORM = 'F', 'f', 'E' or 'e'
+*
+*  where  norm1  denotes the  one norm of a matrix (maximum column sum),
+*  normI  denotes the  infinity norm  of a matrix  (maximum row sum) and
+*  normF  denotes the  Frobenius norm of a matrix (square root of sum of
+*  squares).  Note that  max(abs(A(i,j)))  is not a  matrix norm.
+*
+*  Arguments
+*  =========
+*
+*  NORM    (input) CHARACTER*1
+*          Specifies the value to be returned in DLANST as described
+*          above.
+*
+*  N       (input) INTEGER
+*          The order of the matrix A.  N >= 0.  When N = 0, DLANST is
+*          set to zero.
+*
+*  D       (input) DOUBLE PRECISION array, dimension (N)
+*          The diagonal elements of A.
+*
+*  E       (input) DOUBLE PRECISION array, dimension (N-1)
+*          The (n-1) sub-diagonal or super-diagonal elements of A.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ONE, ZERO
+      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            I
+      DOUBLE PRECISION   ANORM, SCALE, SUM
+*     ..
+*     .. External Functions ..
+      LOGICAL            PLSAME
+      EXTERNAL           PLSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           PDLASSQ
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, MAX, SQRT
+*     ..
+*     .. Executable Statements ..
+*
+      IF( N.LE.0 ) THEN
+         ANORM = ZERO
+      ELSE IF( PLSAME( NORM, 'M' ) ) THEN
+*
+*        Find max(abs(A(i,j))).
+*
+         ANORM = ABS( D( N ) )
+         DO 10 I = 1, N - 1
+            ANORM = MAX( ANORM, ABS( D( I ) ) )
+            ANORM = MAX( ANORM, ABS( E( I ) ) )
+   10    CONTINUE
+      ELSE IF( PLSAME( NORM, 'O' ) .OR. NORM.EQ.'1' .OR.
+     $         PLSAME( NORM, 'I' ) ) THEN
+*
+*        Find norm1(A).
+*
+         IF( N.EQ.1 ) THEN
+            ANORM = ABS( D( 1 ) )
+         ELSE
+            ANORM = MAX( ABS( D( 1 ) )+ABS( E( 1 ) ),
+     $              ABS( E( N-1 ) )+ABS( D( N ) ) )
+            DO 20 I = 2, N - 1
+               ANORM = MAX( ANORM, ABS( D( I ) )+ABS( E( I ) )+
+     $                 ABS( E( I-1 ) ) )
+   20       CONTINUE
+         END IF
+      ELSE IF( ( PLSAME( NORM, 'F' ) ) .OR. ( PLSAME( NORM, 'E'))) THEN
+*
+*        Find normF(A).
+*
+         SCALE = ZERO
+         SUM = ONE
+         IF( N.GT.1 ) THEN
+            CALL PDLASSQ( N-1, E, 1, SCALE, SUM )
+            SUM = 2*SUM
+         END IF
+         CALL PDLASSQ( N, D, 1, SCALE, SUM )
+         ANORM = SCALE*SQRT( SUM )
+      END IF
+*
+      PDLANST = ANORM
+      RETURN
+*
+*     End of DLANST
+*
+      END
+      DOUBLE PRECISION FUNCTION PDLAPY2( X, Y )
+*
+*  -- LAPACK auxiliary routine (version 3.0) --
+*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+*     Courant Institute, Argonne National Lab, and Rice University
+*     October 31, 1992
+*
+*     .. Scalar Arguments ..
+      DOUBLE PRECISION   X, Y
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DLAPY2 returns sqrt(x**2+y**2), taking care not to cause unnecessary
+*  overflow.
+*
+*  Arguments
+*  =========
+*
+*  X       (input) DOUBLE PRECISION
+*  Y       (input) DOUBLE PRECISION
+*          X and Y specify the values x and y.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO
+      PARAMETER          ( ZERO = 0.0D0 )
+      DOUBLE PRECISION   ONE
+      PARAMETER          ( ONE = 1.0D0 )
+*     ..
+*     .. Local Scalars ..
+      DOUBLE PRECISION   W, XABS, YABS, Z
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, MAX, MIN, SQRT
+*     ..
+*     .. Executable Statements ..
+*
+      XABS = ABS( X )
+      YABS = ABS( Y )
+      W = MAX( XABS, YABS )
+      Z = MIN( XABS, YABS )
+      IF( Z.EQ.ZERO ) THEN
+         PDLAPY2 = W
+      ELSE
+         PDLAPY2 = W*SQRT( ONE+( Z / W )**2 )
+      END IF
+      RETURN
+*
+*     End of DLAPY2
+*
+      END
+      SUBROUTINE PDLARTG( 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
+*     September 30, 1994
+*
+*     .. Scalar Arguments ..
+      DOUBLE PRECISION   CS, F, G, R, SN
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DLARTG 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 DROTG,
+*  with the following other differences:
+*     F and G are unchanged on return.
+*     If G=0, then CS=1 and SN=0.
+*     If F=0 and (G .ne. 0), then CS=0 and SN=1 without doing any
+*        floating point operations (saves work in DBDSQR when
+*        there are zeros on the diagonal).
+*
+*  If F exceeds G in magnitude, CS will be positive.
+*
+*  Arguments
+*  =========
+*
+*  F       (input) DOUBLE PRECISION
+*          The first component of vector to be rotated.
+*
+*  G       (input) DOUBLE PRECISION
+*          The second component of vector to be rotated.
+*
+*  CS      (output) DOUBLE PRECISION
+*          The cosine of the rotation.
+*
+*  SN      (output) DOUBLE PRECISION
+*          The sine of the rotation.
+*
+*  R       (output) DOUBLE PRECISION
+*          The nonzero component of the rotated vector.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO
+      PARAMETER          ( ZERO = 0.0D0 )
+      DOUBLE PRECISION   ONE
+      PARAMETER          ( ONE = 1.0D0 )
+      DOUBLE PRECISION   TWO
+      PARAMETER          ( TWO = 2.0D0 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            FIRST
+      INTEGER            COUNT, I
+      DOUBLE PRECISION   EPS, F1, G1, SAFMIN, SAFMN2, SAFMX2, SCALE
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   PDLAMCH
+      EXTERNAL           PDLAMCH
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, INT, LOG, MAX, SQRT
+*     ..
+*     .. Save statement ..
+      SAVE               FIRST, SAFMX2, SAFMIN, SAFMN2
+*     ..
+*     .. Data statements ..
+      DATA               FIRST / .TRUE. /
+*     ..
+*     .. Executable Statements ..
+*
+      IF( FIRST ) THEN
+         FIRST = .FALSE.
+         SAFMIN = PDLAMCH( 'S' )
+         EPS = PDLAMCH( 'E' )
+         SAFMN2 = PDLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) /
+     $            LOG( PDLAMCH( 'B' ) ) / TWO )
+         SAFMX2 = ONE / SAFMN2
+      END IF
+      IF( G.EQ.ZERO ) THEN
+         CS = ONE
+         SN = ZERO
+         R = F
+      ELSE IF( F.EQ.ZERO ) THEN
+         CS = ZERO
+         SN = ONE
+         R = G
+      ELSE
+         F1 = F
+         G1 = G
+         SCALE = MAX( ABS( F1 ), ABS( G1 ) )
+         IF( SCALE.GE.SAFMX2 ) THEN
+            COUNT = 0
+   10       CONTINUE
+            COUNT = COUNT + 1
+            F1 = F1*SAFMN2
+            G1 = G1*SAFMN2
+            SCALE = MAX( ABS( F1 ), ABS( G1 ) )
+            IF( SCALE.GE.SAFMX2 )
+     $         GO TO 10
+            R = SQRT( F1**2+G1**2 )
+            CS = F1 / R
+            SN = G1 / R
+            DO 20 I = 1, COUNT
+               R = R*SAFMX2
+   20       CONTINUE
+         ELSE IF( SCALE.LE.SAFMN2 ) THEN
+            COUNT = 0
+   30       CONTINUE
+            COUNT = COUNT + 1
+            F1 = F1*SAFMX2
+            G1 = G1*SAFMX2
+            SCALE = MAX( ABS( F1 ), ABS( G1 ) )
+            IF( SCALE.LE.SAFMN2 )
+     $         GO TO 30
+            R = SQRT( F1**2+G1**2 )
+            CS = F1 / R
+            SN = G1 / R
+            DO 40 I = 1, COUNT
+               R = R*SAFMN2
+   40       CONTINUE
+         ELSE
+            R = SQRT( F1**2+G1**2 )
+            CS = F1 / R
+            SN = G1 / R
+         END IF
+         IF( ABS( F ).GT.ABS( G ) .AND. CS.LT.ZERO ) THEN
+            CS = -CS
+            SN = -SN
+            R = -R
+         END IF
+      END IF
+      RETURN
+*
+*     End of DLARTG
+*
+      END
+      SUBROUTINE PDLASCL( TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO )
+*
+*  -- LAPACK auxiliary routine (version 3.0) --
+*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+*     Courant Institute, Argonne National Lab, and Rice University
+*     February 29, 1992
+*
+*     .. Scalar Arguments ..
+      CHARACTER          TYPE
+      INTEGER            INFO, KL, KU, LDA, M, N
+      DOUBLE PRECISION   CFROM, CTO
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   A( LDA, * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DLASCL multiplies the M by N real matrix A by the real scalar
+*  CTO/CFROM.  This is done without over/underflow as long as the final
+*  result CTO*A(I,J)/CFROM does not over/underflow. TYPE specifies that
+*  A may be full, upper triangular, lower triangular, upper Hessenberg,
+*  or banded.
+*
+*  Arguments
+*  =========
+*
+*  TYPE    (input) CHARACTER*1
+*          TYPE indices the storage type of the input matrix.
+*          = 'G':  A is a full matrix.
+*          = 'L':  A is a lower triangular matrix.
+*          = 'U':  A is an upper triangular matrix.
+*          = 'H':  A is an upper Hessenberg matrix.
+*          = 'B':  A is a symmetric band matrix with lower bandwidth KL
+*                  and upper bandwidth KU and with the only the lower
+*                  half stored.
+*          = 'Q':  A is a symmetric band matrix with lower bandwidth KL
+*                  and upper bandwidth KU and with the only the upper
+*                  half stored.
+*          = 'Z':  A is a band matrix with lower bandwidth KL and upper
+*                  bandwidth KU.
+*
+*  KL      (input) INTEGER
+*          The lower bandwidth of A.  Referenced only if TYPE = 'B',
+*          'Q' or 'Z'.
+*
+*  KU      (input) INTEGER
+*          The upper bandwidth of A.  Referenced only if TYPE = 'B',
+*          'Q' or 'Z'.
+*
+*  CFROM   (input) DOUBLE PRECISION
+*  CTO     (input) DOUBLE PRECISION
+*          The matrix A is multiplied by CTO/CFROM. A(I,J) is computed
+*          without over/underflow if the final result CTO*A(I,J)/CFROM
+*          can be represented without over/underflow.  CFROM must be
+*          nonzero.
+*
+*  M       (input) INTEGER
+*          The number of rows of the matrix A.  M >= 0.
+*
+*  N       (input) INTEGER
+*          The number of columns of the matrix A.  N >= 0.
+*
+*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,M)
+*          The matrix to be multiplied by CTO/CFROM.  See TYPE for the
+*          storage type.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,M).
+*
+*  INFO    (output) INTEGER
+*          0  - successful exit
+*          <0 - if INFO = -i, the i-th argument had an illegal value.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            DONE
+      INTEGER            I, ITYPE, J, K1, K2, K3, K4
+      DOUBLE PRECISION   BIGNUM, CFROM1, CFROMC, CTO1, CTOC, MUL, SMLNUM
+*     ..
+*     .. External Functions ..
+      LOGICAL            PLSAME
+      DOUBLE PRECISION   PDLAMCH
+      EXTERNAL           PLSAME, PDLAMCH
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, MAX, MIN
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           PXERBLA
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input arguments
+*
+      INFO = 0
+*
+      IF( PLSAME( TYPE, 'G' ) ) THEN
+         ITYPE = 0
+      ELSE IF( PLSAME( TYPE, 'L' ) ) THEN
+         ITYPE = 1
+      ELSE IF( PLSAME( TYPE, 'U' ) ) THEN
+         ITYPE = 2
+      ELSE IF( PLSAME( TYPE, 'H' ) ) THEN
+         ITYPE = 3
+      ELSE IF( PLSAME( TYPE, 'B' ) ) THEN
+         ITYPE = 4
+      ELSE IF( PLSAME( TYPE, 'Q' ) ) THEN
+         ITYPE = 5
+      ELSE IF( PLSAME( TYPE, 'Z' ) ) THEN
+         ITYPE = 6
+      ELSE
+         ITYPE = -1
+      END IF
+*
+      IF( ITYPE.EQ.-1 ) THEN
+         INFO = -1
+      ELSE IF( CFROM.EQ.ZERO ) THEN
+         INFO = -4
+      ELSE IF( M.LT.0 ) THEN
+         INFO = -6
+      ELSE IF( N.LT.0 .OR. ( ITYPE.EQ.4 .AND. N.NE.M ) .OR.
+     $         ( ITYPE.EQ.5 .AND. N.NE.M ) ) THEN
+         INFO = -7
+      ELSE IF( ITYPE.LE.3 .AND. LDA.LT.MAX( 1, M ) ) THEN
+         INFO = -9
+      ELSE IF( ITYPE.GE.4 ) THEN
+         IF( KL.LT.0 .OR. KL.GT.MAX( M-1, 0 ) ) THEN
+            INFO = -2
+         ELSE IF( KU.LT.0 .OR. KU.GT.MAX( N-1, 0 ) .OR.
+     $            ( ( ITYPE.EQ.4 .OR. ITYPE.EQ.5 ) .AND. KL.NE.KU ) )
+     $             THEN
+            INFO = -3
+         ELSE IF( ( ITYPE.EQ.4 .AND. LDA.LT.KL+1 ) .OR.
+     $            ( ITYPE.EQ.5 .AND. LDA.LT.KU+1 ) .OR.
+     $            ( ITYPE.EQ.6 .AND. LDA.LT.2*KL+KU+1 ) ) THEN
+            INFO = -9
+         END IF
+      END IF
+*
+      IF( INFO.NE.0 ) THEN
+         CALL PXERBLA( 'PDLASCL', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 .OR. M.EQ.0 )
+     $   RETURN
+*
+*     Get machine parameters
+*
+      SMLNUM = PDLAMCH( 'S' )
+      BIGNUM = ONE / SMLNUM
+*
+      CFROMC = CFROM
+      CTOC = CTO
+*
+   10 CONTINUE
+      CFROM1 = CFROMC*SMLNUM
+      CTO1 = CTOC / BIGNUM
+      IF( ABS( CFROM1 ).GT.ABS( CTOC ) .AND. CTOC.NE.ZERO ) THEN
+         MUL = SMLNUM
+         DONE = .FALSE.
+         CFROMC = CFROM1
+      ELSE IF( ABS( CTO1 ).GT.ABS( CFROMC ) ) THEN
+         MUL = BIGNUM
+         DONE = .FALSE.
+         CTOC = CTO1
+      ELSE
+         MUL = CTOC / CFROMC
+         DONE = .TRUE.
+      END IF
+*
+      IF( ITYPE.EQ.0 ) THEN
+*
+*        Full matrix
+*
+         DO 30 J = 1, N
+            DO 20 I = 1, M
+               A( I, J ) = A( I, J )*MUL
+   20       CONTINUE
+   30    CONTINUE
+*
+      ELSE IF( ITYPE.EQ.1 ) THEN
+*
+*        Lower triangular matrix
+*
+         DO 50 J = 1, N
+            DO 40 I = J, M
+               A( I, J ) = A( I, J )*MUL
+   40       CONTINUE
+   50    CONTINUE
+*
+      ELSE IF( ITYPE.EQ.2 ) THEN
+*
+*        Upper triangular matrix
+*
+         DO 70 J = 1, N
+            DO 60 I = 1, MIN( J, M )
+               A( I, J ) = A( I, J )*MUL
+   60       CONTINUE
+   70    CONTINUE
+*
+      ELSE IF( ITYPE.EQ.3 ) THEN
+*
+*        Upper Hessenberg matrix
+*
+         DO 90 J = 1, N
+            DO 80 I = 1, MIN( J+1, M )
+               A( I, J ) = A( I, J )*MUL
+   80       CONTINUE
+   90    CONTINUE
+*
+      ELSE IF( ITYPE.EQ.4 ) THEN
+*
+*        Lower half of a symmetric band matrix
+*
+         K3 = KL + 1
+         K4 = N + 1
+         DO 110 J = 1, N
+            DO 100 I = 1, MIN( K3, K4-J )
+               A( I, J ) = A( I, J )*MUL
+  100       CONTINUE
+  110    CONTINUE
+*
+      ELSE IF( ITYPE.EQ.5 ) THEN
+*
+*        Upper half of a symmetric band matrix
+*
+         K1 = KU + 2
+         K3 = KU + 1
+         DO 130 J = 1, N
+            DO 120 I = MAX( K1-J, 1 ), K3
+               A( I, J ) = A( I, J )*MUL
+  120       CONTINUE
+  130    CONTINUE
+*
+      ELSE IF( ITYPE.EQ.6 ) THEN
+*
+*        Band matrix
+*
+         K1 = KL + KU + 2
+         K2 = KL + 1
+         K3 = 2*KL + KU + 1
+         K4 = KL + KU + 1 + M
+         DO 150 J = 1, N
+            DO 140 I = MAX( K1-J, K2 ), MIN( K3, K4-J )
+               A( I, J ) = A( I, J )*MUL
+  140       CONTINUE
+  150    CONTINUE
+*
+      END IF
+*
+      IF( .NOT.DONE )
+     $   GO TO 10
+*
+      RETURN
+*
+*     End of DLASCL
+*
+      END
+      SUBROUTINE PDLASET( UPLO, M, N, ALPHA, BETA, A, LDA )
+*
+*  -- LAPACK auxiliary routine (version 3.0) --
+*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+*     Courant Institute, Argonne National Lab, and Rice University
+*     October 31, 1992
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            LDA, M, N
+      DOUBLE PRECISION   ALPHA, BETA
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   A( LDA, * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DLASET initializes an m-by-n matrix A to BETA on the diagonal and
+*  ALPHA on the offdiagonals.
+*
+*  Arguments
+*  =========
+*
+*  UPLO    (input) CHARACTER*1
+*          Specifies the part of the matrix A to be set.
+*          = 'U':      Upper triangular part is set; the strictly lower
+*                      triangular part of A is not changed.
+*          = 'L':      Lower triangular part is set; the strictly upper
+*                      triangular part of A is not changed.
+*          Otherwise:  All of the matrix A is set.
+*
+*  M       (input) INTEGER
+*          The number of rows of the matrix A.  M >= 0.
+*
+*  N       (input) INTEGER
+*          The number of columns of the matrix A.  N >= 0.
+*
+*  ALPHA   (input) DOUBLE PRECISION
+*          The constant to which the offdiagonal elements are to be set.
+*
+*  BETA    (input) DOUBLE PRECISION
+*          The constant to which the diagonal elements are to be set.
+*
+*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
+*          On exit, the leading m-by-n submatrix of A is set as follows:
+*
+*          if UPLO = 'U', A(i,j) = ALPHA, 1<=i<=j-1, 1<=j<=n,
+*          if UPLO = 'L', A(i,j) = ALPHA, j+1<=i<=m, 1<=j<=n,
+*          otherwise,     A(i,j) = ALPHA, 1<=i<=m, 1<=j<=n, i.ne.j,
+*
+*          and, for all UPLO, A(i,i) = BETA, 1<=i<=min(m,n).
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,M).
+*
+* =====================================================================
+*
+*     .. Local Scalars ..
+      INTEGER            I, J
+*     ..
+*     .. External Functions ..
+      LOGICAL            PLSAME
+      EXTERNAL           PLSAME
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MIN
+*     ..
+*     .. Executable Statements ..
+*
+      IF( PLSAME( UPLO, 'U' ) ) THEN
+*
+*        Set the strictly upper triangular or trapezoidal part of the
+*        array to ALPHA.
+*
+         DO 20 J = 2, N
+            DO 10 I = 1, MIN( J-1, M )
+               A( I, J ) = ALPHA
+   10       CONTINUE
+   20    CONTINUE
+*
+      ELSE IF( PLSAME( UPLO, 'L' ) ) THEN
+*
+*        Set the strictly lower triangular or trapezoidal part of the
+*        array to ALPHA.
+*
+         DO 40 J = 1, MIN( M, N )
+            DO 30 I = J + 1, M
+               A( I, J ) = ALPHA
+   30       CONTINUE
+   40    CONTINUE
+*
+      ELSE
+*
+*        Set the leading m-by-n submatrix to ALPHA.
+*
+         DO 60 J = 1, N
+            DO 50 I = 1, M
+               A( I, J ) = ALPHA
+   50       CONTINUE
+   60    CONTINUE
+      END IF
+*
+*     Set the first min(M,N) diagonal elements to BETA.
+*
+      DO 70 I = 1, MIN( M, N )
+         A( I, I ) = BETA
+   70 CONTINUE
+*
+      RETURN
+*
+*     End of DLASET
+*
+      END
+      SUBROUTINE PDLASR( SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA )
+*
+*  -- LAPACK auxiliary routine (version 3.0) --
+*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+*     Courant Institute, Argonne National Lab, and Rice University
+*     October 31, 1992
+*
+*     .. Scalar Arguments ..
+      CHARACTER          DIRECT, PIVOT, SIDE
+      INTEGER            LDA, M, N
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   A( LDA, * ), C( * ), S( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DLASR   performs the transformation
+*
+*     A := P*A,   when SIDE = 'L' or 'l'  (  Left-hand side )
+*
+*     A := A*P',  when SIDE = 'R' or 'r'  ( Right-hand side )
+*
+*  where A is an m by n real matrix and P is an orthogonal matrix,
+*  consisting of a sequence of plane rotations determined by the
+*  parameters PIVOT and DIRECT as follows ( z = m when SIDE = 'L' or 'l'
+*  and z = n when SIDE = 'R' or 'r' ):
+*
+*  When  DIRECT = 'F' or 'f'  ( Forward sequence ) then
+*
+*     P = P( z - 1 )*...*P( 2 )*P( 1 ),
+*
+*  and when DIRECT = 'B' or 'b'  ( Backward sequence ) then
+*
+*     P = P( 1 )*P( 2 )*...*P( z - 1 ),
+*
+*  where  P( k ) is a plane rotation matrix for the following planes:
+*
+*     when  PIVOT = 'V' or 'v'  ( Variable pivot ),
+*        the plane ( k, k + 1 )
+*
+*     when  PIVOT = 'T' or 't'  ( Top pivot ),
+*        the plane ( 1, k + 1 )
+*
+*     when  PIVOT = 'B' or 'b'  ( Bottom pivot ),
+*        the plane ( k, z )
+*
+*  c( k ) and s( k )  must contain the  cosine and sine that define the
+*  matrix  P( k ).  The two by two plane rotation part of the matrix
+*  P( k ), R( k ), is assumed to be of the form
+*
+*     R( k ) = (  c( k )  s( k ) ).
+*              ( -s( k )  c( k ) )
+*
+*  This version vectorises across rows of the array A when SIDE = 'L'.
+*
+*  Arguments
+*  =========
+*
+*  SIDE    (input) CHARACTER*1
+*          Specifies whether the plane rotation matrix P is applied to
+*          A on the left or the right.
+*          = 'L':  Left, compute A := P*A
+*          = 'R':  Right, compute A:= A*P'
+*
+*  DIRECT  (input) CHARACTER*1
+*          Specifies whether P is a forward or backward sequence of
+*          plane rotations.
+*          = 'F':  Forward, P = P( z - 1 )*...*P( 2 )*P( 1 )
+*          = 'B':  Backward, P = P( 1 )*P( 2 )*...*P( z - 1 )
+*
+*  PIVOT   (input) CHARACTER*1
+*          Specifies the plane for which P(k) is a plane rotation
+*          matrix.
+*          = 'V':  Variable pivot, the plane (k,k+1)
+*          = 'T':  Top pivot, the plane (1,k+1)
+*          = 'B':  Bottom pivot, the plane (k,z)
+*
+*  M       (input) INTEGER
+*          The number of rows of the matrix A.  If m <= 1, an immediate
+*          return is effected.
+*
+*  N       (input) INTEGER
+*          The number of columns of the matrix A.  If n <= 1, an
+*          immediate return is effected.
+*
+*  C, S    (input) DOUBLE PRECISION arrays, dimension
+*                  (M-1) if SIDE = 'L'
+*                  (N-1) if SIDE = 'R'
+*          c(k) and s(k) contain the cosine and sine that define the
+*          matrix P(k).  The two by two plane rotation part of the
+*          matrix P(k), R(k), is assumed to be of the form
+*          R( k ) = (  c( k )  s( k ) ).
+*                   ( -s( k )  c( k ) )
+*
+*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
+*          The m by n matrix A.  On exit, A is overwritten by P*A if
+*          SIDE = 'R' or by A*P' if SIDE = 'L'.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,M).
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ONE, ZERO
+      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            I, INFO, J
+      DOUBLE PRECISION   CTEMP, STEMP, TEMP
+*     ..
+*     .. External Functions ..
+      LOGICAL            PLSAME
+      EXTERNAL           PLSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           PXERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input parameters
+*
+      INFO = 0
+      IF( .NOT.( PLSAME( SIDE, 'L' ) .OR. PLSAME( SIDE, 'R' ) ) ) THEN
+         INFO = 1
+      ELSE IF( .NOT.( PLSAME( PIVOT, 'V' ) .OR. PLSAME( PIVOT,
+     $         'T' ) .OR. PLSAME( PIVOT, 'B' ) ) ) THEN
+         INFO = 2
+      ELSE IF( .NOT.( PLSAME( DIRECT, 'F' ) .OR. PLSAME( DIRECT, 'B' )))
+     $          THEN
+         INFO = 3
+      ELSE IF( M.LT.0 ) THEN
+         INFO = 4
+      ELSE IF( N.LT.0 ) THEN
+         INFO = 5
+      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
+         INFO = 9
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL PXERBLA( 'PDLASR ', INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) )
+     $   RETURN
+      IF( PLSAME( SIDE, 'L' ) ) THEN
+*
+*        Form  P * A
+*
+         IF( PLSAME( PIVOT, 'V' ) ) THEN
+            IF( PLSAME( DIRECT, 'F' ) ) THEN
+               DO 20 J = 1, M - 1
+                  CTEMP = C( J )
+                  STEMP = S( J )
+                  IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
+                     DO 10 I = 1, N
+                        TEMP = A( J+1, I )
+                        A( J+1, I ) = CTEMP*TEMP - STEMP*A( J, I )
+                        A( J, I ) = STEMP*TEMP + CTEMP*A( J, I )
+   10                CONTINUE
+                  END IF
+   20          CONTINUE
+            ELSE IF( PLSAME( DIRECT, 'B' ) ) THEN
+               DO 40 J = M - 1, 1, -1
+                  CTEMP = C( J )
+                  STEMP = S( J )
+                  IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
+                     DO 30 I = 1, N
+                        TEMP = A( J+1, I )
+                        A( J+1, I ) = CTEMP*TEMP - STEMP*A( J, I )
+                        A( J, I ) = STEMP*TEMP + CTEMP*A( J, I )
+   30                CONTINUE
+                  END IF
+   40          CONTINUE
+            END IF
+         ELSE IF( PLSAME( PIVOT, 'T' ) ) THEN
+            IF( PLSAME( DIRECT, 'F' ) ) THEN
+               DO 60 J = 2, M
+                  CTEMP = C( J-1 )
+                  STEMP = S( J-1 )
+                  IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
+                     DO 50 I = 1, N
+                        TEMP = A( J, I )
+                        A( J, I ) = CTEMP*TEMP - STEMP*A( 1, I )
+                        A( 1, I ) = STEMP*TEMP + CTEMP*A( 1, I )
+   50                CONTINUE
+                  END IF
+   60          CONTINUE
+            ELSE IF( PLSAME( DIRECT, 'B' ) ) THEN
+               DO 80 J = M, 2, -1
+                  CTEMP = C( J-1 )
+                  STEMP = S( J-1 )
+                  IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
+                     DO 70 I = 1, N
+                        TEMP = A( J, I )
+                        A( J, I ) = CTEMP*TEMP - STEMP*A( 1, I )
+                        A( 1, I ) = STEMP*TEMP + CTEMP*A( 1, I )
+   70                CONTINUE
+                  END IF
+   80          CONTINUE
+            END IF
+         ELSE IF( PLSAME( PIVOT, 'B' ) ) THEN
+            IF( PLSAME( DIRECT, 'F' ) ) THEN
+               DO 100 J = 1, M - 1
+                  CTEMP = C( J )
+                  STEMP = S( J )
+                  IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
+                     DO 90 I = 1, N
+                        TEMP = A( J, I )
+                        A( J, I ) = STEMP*A( M, I ) + CTEMP*TEMP
+                        A( M, I ) = CTEMP*A( M, I ) - STEMP*TEMP
+   90                CONTINUE
+                  END IF
+  100          CONTINUE
+            ELSE IF( PLSAME( DIRECT, 'B' ) ) THEN
+               DO 120 J = M - 1, 1, -1
+                  CTEMP = C( J )
+                  STEMP = S( J )
+                  IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
+                     DO 110 I = 1, N
+                        TEMP = A( J, I )
+                        A( J, I ) = STEMP*A( M, I ) + CTEMP*TEMP
+                        A( M, I ) = CTEMP*A( M, I ) - STEMP*TEMP
+  110                CONTINUE
+                  END IF
+  120          CONTINUE
+            END IF
+         END IF
+      ELSE IF( PLSAME( SIDE, 'R' ) ) THEN
+*
+*        Form A * P'
+*
+         IF( PLSAME( PIVOT, 'V' ) ) THEN
+            IF( PLSAME( DIRECT, 'F' ) ) THEN
+               DO 140 J = 1, N - 1
+                  CTEMP = C( J )
+                  STEMP = S( J )
+                  IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
+                     DO 130 I = 1, M
+                        TEMP = A( I, J+1 )
+                        A( I, J+1 ) = CTEMP*TEMP - STEMP*A( I, J )
+                        A( I, J ) = STEMP*TEMP + CTEMP*A( I, J )
+  130                CONTINUE
+                  END IF
+  140          CONTINUE
+            ELSE IF( PLSAME( DIRECT, 'B' ) ) THEN
+               DO 160 J = N - 1, 1, -1
+                  CTEMP = C( J )
+                  STEMP = S( J )
+                  IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
+                     DO 150 I = 1, M
+                        TEMP = A( I, J+1 )
+                        A( I, J+1 ) = CTEMP*TEMP - STEMP*A( I, J )
+                        A( I, J ) = STEMP*TEMP + CTEMP*A( I, J )
+  150                CONTINUE
+                  END IF
+  160          CONTINUE
+            END IF
+         ELSE IF( PLSAME( PIVOT, 'T' ) ) THEN
+            IF( PLSAME( DIRECT, 'F' ) ) THEN
+               DO 180 J = 2, N
+                  CTEMP = C( J-1 )
+                  STEMP = S( J-1 )
+                  IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
+                     DO 170 I = 1, M
+                        TEMP = A( I, J )
+                        A( I, J ) = CTEMP*TEMP - STEMP*A( I, 1 )
+                        A( I, 1 ) = STEMP*TEMP + CTEMP*A( I, 1 )
+  170                CONTINUE
+                  END IF
+  180          CONTINUE
+            ELSE IF( PLSAME( DIRECT, 'B' ) ) THEN
+               DO 200 J = N, 2, -1
+                  CTEMP = C( J-1 )
+                  STEMP = S( J-1 )
+                  IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
+                     DO 190 I = 1, M
+                        TEMP = A( I, J )
+                        A( I, J ) = CTEMP*TEMP - STEMP*A( I, 1 )
+                        A( I, 1 ) = STEMP*TEMP + CTEMP*A( I, 1 )
+  190                CONTINUE
+                  END IF
+  200          CONTINUE
+            END IF
+         ELSE IF( PLSAME( PIVOT, 'B' ) ) THEN
+            IF( PLSAME( DIRECT, 'F' ) ) THEN
+               DO 220 J = 1, N - 1
+                  CTEMP = C( J )
+                  STEMP = S( J )
+                  IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
+                     DO 210 I = 1, M
+                        TEMP = A( I, J )
+                        A( I, J ) = STEMP*A( I, N ) + CTEMP*TEMP
+                        A( I, N ) = CTEMP*A( I, N ) - STEMP*TEMP
+  210                CONTINUE
+                  END IF
+  220          CONTINUE
+            ELSE IF( PLSAME( DIRECT, 'B' ) ) THEN
+               DO 240 J = N - 1, 1, -1
+                  CTEMP = C( J )
+                  STEMP = S( J )
+                  IF( ( CTEMP.NE.ONE ) .OR. ( STEMP.NE.ZERO ) ) THEN
+                     DO 230 I = 1, M
+                        TEMP = A( I, J )
+                        A( I, J ) = STEMP*A( I, N ) + CTEMP*TEMP
+                        A( I, N ) = CTEMP*A( I, N ) - STEMP*TEMP
+  230                CONTINUE
+                  END IF
+  240          CONTINUE
+            END IF
+         END IF
+      END IF
+*
+      RETURN
+*
+*     End of DLASR
+*
+      END
+      SUBROUTINE PDLASRT( ID, N, D, INFO )
+*
+*  -- LAPACK routine (version 3.0) --
+*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+*     Courant Institute, Argonne National Lab, and Rice University
+*     September 30, 1994
+*
+*     .. Scalar Arguments ..
+      CHARACTER          ID
+      INTEGER            INFO, N
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   D( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  Sort the numbers in D in increasing order (if ID = 'I') or
+*  in decreasing order (if ID = 'D' ).
+*
+*  Use Quick Sort, reverting to Insertion sort on arrays of
+*  size <= 20. Dimension of STACK limits N to about 2**32.
+*
+*  Arguments
+*  =========
+*
+*  ID      (input) CHARACTER*1
+*          = 'I': sort D in increasing order;
+*          = 'D': sort D in decreasing order.
+*
+*  N       (input) INTEGER
+*          The length of the array D.
+*
+*  D       (input/output) DOUBLE PRECISION array, dimension (N)
+*          On entry, the array to be sorted.
+*          On exit, D has been sorted into increasing order
+*          (D(1) <= ... <= D(N) ) or into decreasing order
+*          (D(1) >= ... >= D(N) ), depending on ID.
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      INTEGER            SELECT
+      PARAMETER          ( SELECT = 20 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            DIR, ENDD, I, J, START, STKPNT
+      DOUBLE PRECISION   D1, D2, D3, DMNMX, TMP
+*     ..
+*     .. Local Arrays ..
+      INTEGER            STACK( 2, 32 )
+*     ..
+*     .. External Functions ..
+      LOGICAL            PLSAME
+      EXTERNAL           PLSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           PXERBLA
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input paramters.
+*
+      INFO = 0
+      DIR = -1
+      IF( PLSAME( ID, 'D' ) ) THEN
+         DIR = 0
+      ELSE IF( PLSAME( ID, 'I' ) ) THEN
+         DIR = 1
+      END IF
+      IF( DIR.EQ.-1 ) THEN
+         INFO = -1
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -2
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL PXERBLA( 'PDLASRT', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.LE.1 )
+     $   RETURN
+*
+      STKPNT = 1
+      STACK( 1, 1 ) = 1
+      STACK( 2, 1 ) = N
+   10 CONTINUE
+      START = STACK( 1, STKPNT )
+      ENDD = STACK( 2, STKPNT )
+      STKPNT = STKPNT - 1
+      IF( ENDD-START.LE.SELECT .AND. ENDD-START.GT.0 ) THEN
+*
+*        Do Insertion sort on D( START:ENDD )
+*
+         IF( DIR.EQ.0 ) THEN
+*
+*           Sort into decreasing order
+*
+            DO 30 I = START + 1, ENDD
+               DO 20 J = I, START + 1, -1
+                  IF( D( J ).GT.D( J-1 ) ) THEN
+                     DMNMX = D( J )
+                     D( J ) = D( J-1 )
+                     D( J-1 ) = DMNMX
+                  ELSE
+                     GO TO 30
+                  END IF
+   20          CONTINUE
+   30       CONTINUE
+*
+         ELSE
+*
+*           Sort into increasing order
+*
+            DO 50 I = START + 1, ENDD
+               DO 40 J = I, START + 1, -1
+                  IF( D( J ).LT.D( J-1 ) ) THEN
+                     DMNMX = D( J )
+                     D( J ) = D( J-1 )
+                     D( J-1 ) = DMNMX
+                  ELSE
+                     GO TO 50
+                  END IF
+   40          CONTINUE
+   50       CONTINUE
+*
+         END IF
+*
+      ELSE IF( ENDD-START.GT.SELECT ) THEN
+*
+*        Partition D( START:ENDD ) and stack parts, largest one first
+*
+*        Choose partition entry as median of 3
+*
+         D1 = D( START )
+         D2 = D( ENDD )
+         I = ( START+ENDD ) / 2
+         D3 = D( I )
+         IF( D1.LT.D2 ) THEN
+            IF( D3.LT.D1 ) THEN
+               DMNMX = D1
+            ELSE IF( D3.LT.D2 ) THEN
+               DMNMX = D3
+            ELSE
+               DMNMX = D2
+            END IF
+         ELSE
+            IF( D3.LT.D2 ) THEN
+               DMNMX = D2
+            ELSE IF( D3.LT.D1 ) THEN
+               DMNMX = D3
+            ELSE
+               DMNMX = D1
+            END IF
+         END IF
+*
+         IF( DIR.EQ.0 ) THEN
+*
+*           Sort into decreasing order
+*
+            I = START - 1
+            J = ENDD + 1
+   60       CONTINUE
+   70       CONTINUE
+            J = J - 1
+            IF( D( J ).LT.DMNMX )
+     $         GO TO 70
+   80       CONTINUE
+            I = I + 1
+            IF( D( I ).GT.DMNMX )
+     $         GO TO 80
+            IF( I.LT.J ) THEN
+               TMP = D( I )
+               D( I ) = D( J )
+               D( J ) = TMP
+               GO TO 60
+            END IF
+            IF( J-START.GT.ENDD-J-1 ) THEN
+               STKPNT = STKPNT + 1
+               STACK( 1, STKPNT ) = START
+               STACK( 2, STKPNT ) = J
+               STKPNT = STKPNT + 1
+               STACK( 1, STKPNT ) = J + 1
+               STACK( 2, STKPNT ) = ENDD
+            ELSE
+               STKPNT = STKPNT + 1
+               STACK( 1, STKPNT ) = J + 1
+               STACK( 2, STKPNT ) = ENDD
+               STKPNT = STKPNT + 1
+               STACK( 1, STKPNT ) = START
+               STACK( 2, STKPNT ) = J
+            END IF
+         ELSE
+*
+*           Sort into increasing order
+*
+            I = START - 1
+            J = ENDD + 1
+   90       CONTINUE
+  100       CONTINUE
+            J = J - 1
+            IF( D( J ).GT.DMNMX )
+     $         GO TO 100
+  110       CONTINUE
+            I = I + 1
+            IF( D( I ).LT.DMNMX )
+     $         GO TO 110
+            IF( I.LT.J ) THEN
+               TMP = D( I )
+               D( I ) = D( J )
+               D( J ) = TMP
+               GO TO 90
+            END IF
+            IF( J-START.GT.ENDD-J-1 ) THEN
+               STKPNT = STKPNT + 1
+               STACK( 1, STKPNT ) = START
+               STACK( 2, STKPNT ) = J
+               STKPNT = STKPNT + 1
+               STACK( 1, STKPNT ) = J + 1
+               STACK( 2, STKPNT ) = ENDD
+            ELSE
+               STKPNT = STKPNT + 1
+               STACK( 1, STKPNT ) = J + 1
+               STACK( 2, STKPNT ) = ENDD
+               STKPNT = STKPNT + 1
+               STACK( 1, STKPNT ) = START
+               STACK( 2, STKPNT ) = J
+            END IF
+         END IF
+      END IF
+      IF( STKPNT.GT.0 )
+     $   GO TO 10
+      RETURN
+*
+*     End of PDLASRT
+*
+      END
+      SUBROUTINE PDLASSQ( N, X, INCX, SCALE, SUMSQ )
+*
+*  -- LAPACK auxiliary routine (version 3.0) --
+*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+*     Courant Institute, Argonne National Lab, and Rice University
+*     June 30, 1999
+*
+*     .. Scalar Arguments ..
+      INTEGER            INCX, N
+      DOUBLE PRECISION   SCALE, SUMSQ
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   X( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DLASSQ  returns the values  scl  and  smsq  such that
+*
+*     ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq,
+*
+*  where  x( i ) = X( 1 + ( i - 1 )*INCX ). The value of  sumsq  is
+*  assumed to be non-negative and  scl  returns the value
+*
+*     scl = max( scale, abs( x( i ) ) ).
+*
+*  scale and sumsq must be supplied in SCALE and SUMSQ and
+*  scl and smsq are overwritten on SCALE and SUMSQ respectively.
+*
+*  The routine makes only one pass through the vector x.
+*
+*  Arguments
+*  =========
+*
+*  N       (input) INTEGER
+*          The number of elements to be used from the vector X.
+*
+*  X       (input) DOUBLE PRECISION array, dimension (N)
+*          The vector for which a scaled sum of squares is computed.
+*             x( i )  = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n.
+*
+*  INCX    (input) INTEGER
+*          The increment between successive values of the vector X.
+*          INCX > 0.
+*
+*  SCALE   (input/output) DOUBLE PRECISION
+*          On entry, the value  scale  in the equation above.
+*          On exit, SCALE is overwritten with  scl , the scaling factor
+*          for the sum of squares.
+*
+*  SUMSQ   (input/output) DOUBLE PRECISION
+*          On entry, the value  sumsq  in the equation above.
+*          On exit, SUMSQ is overwritten with  smsq , the basic sum of
+*          squares from which  scl  has been factored out.
+*
+* =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO
+      PARAMETER          ( ZERO = 0.0D+0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            IX
+      DOUBLE PRECISION   ABSXI
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS
+*     ..
+*     .. Executable Statements ..
+*
+      IF( N.GT.0 ) THEN
+         DO 10 IX = 1, 1 + ( N-1 )*INCX, INCX
+            IF( X( IX ).NE.ZERO ) THEN
+               ABSXI = ABS( X( IX ) )
+               IF( SCALE.LT.ABSXI ) THEN
+                  SUMSQ = 1 + SUMSQ*( SCALE / ABSXI )**2
+                  SCALE = ABSXI
+               ELSE
+                  SUMSQ = SUMSQ + ( ABSXI / SCALE )**2
+               END IF
+            END IF
+   10    CONTINUE
+      END IF
+      RETURN
+*
+*     End of DLASSQ
+*
+      END
+      LOGICAL          FUNCTION PLSAME( CA, CB )
+*
+*  -- LAPACK auxiliary routine (version 3.0) --
+*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+*     Courant Institute, Argonne National Lab, and Rice University
+*     September 30, 1994
+*
+*     .. Scalar Arguments ..
+      CHARACTER          CA, CB
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  LSAME returns .TRUE. if CA is the same letter as CB regardless of
+*  case.
+*
+*  Arguments
+*  =========
+*
+*  CA      (input) CHARACTER*1
+*  CB      (input) CHARACTER*1
+*          CA and CB specify the single characters to be compared.
+*
+* =====================================================================
+*
+*     .. Intrinsic Functions ..
+      INTRINSIC          ICHAR
+*     ..
+*     .. Local Scalars ..
+      INTEGER            INTA, INTB, ZCODE
+*     ..
+*     .. Executable Statements ..
+*
+*     Test if the characters are equal
+*
+      PLSAME = CA.EQ.CB
+      IF( PLSAME )
+     $   RETURN
+*
+*     Now test for equivalence if both characters are alphabetic.
+*
+      ZCODE = ICHAR( 'Z' )
+*
+*     Use 'Z' rather than 'A' so that ASCII can be detected on Prime
+*     machines, on which ICHAR returns a value with bit 8 set.
+*     ICHAR('A') on Prime machines returns 193 which is the same as
+*     ICHAR('A') on an EBCDIC machine.
+*
+      INTA = ICHAR( CA )
+      INTB = ICHAR( CB )
+*
+      IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN
+*
+*        ASCII is assumed - ZCODE is the ASCII code of either lower or
+*        upper case 'Z'.
+*
+         IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32
+         IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32
+*
+      ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN
+*
+*        EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or
+*        upper case 'Z'.
+*
+         IF( INTA.GE.129 .AND. INTA.LE.137 .OR.
+     $       INTA.GE.145 .AND. INTA.LE.153 .OR.
+     $       INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64
+         IF( INTB.GE.129 .AND. INTB.LE.137 .OR.
+     $       INTB.GE.145 .AND. INTB.LE.153 .OR.
+     $       INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64
+*
+      ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN
+*
+*        ASCII is assumed, on Prime machines - ZCODE is the ASCII code
+*        plus 128 of either lower or upper case 'Z'.
+*
+         IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32
+         IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32
+      END IF
+      PLSAME = INTA.EQ.INTB
+*
+*     RETURN
+*
+*     End of LSAME
+*
+      END
+      SUBROUTINE PXERBLA( SRNAME, INFO )
+*
+*  -- LAPACK auxiliary routine (version 3.0) --
+*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+*     Courant Institute, Argonne National Lab, and Rice University
+*     September 30, 1994
+*
+*     .. Scalar Arguments ..
+      CHARACTER*6        SRNAME
+      INTEGER            INFO
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  XERBLA  is an error handler for the LAPACK routines.
+*  It is called by an LAPACK routine if an input parameter has an
+*  invalid value.  A message is printed and execution stops.
+*
+*  Installers may consider modifying the STOP statement in order to
+*  call system-specific exception-handling facilities.
+*
+*  Arguments
+*  =========
+*
+*  SRNAME  (input) CHARACTER*6
+*          The name of the routine which called XERBLA.
+*
+*  INFO    (input) INTEGER
+*          The position of the invalid parameter in the parameter list
+*          of the calling routine.
+*
+* =====================================================================
+*
+*     .. Executable Statements ..
+*
+      WRITE( *, FMT = 9999 )SRNAME, INFO
+*
+      STOP
+*
+ 9999 FORMAT( ' ** On entry to ', A6, ' parameter number ', I2, ' had ',
+     $      'an illegal value' )
+*
+*     End of XERBLA
+*
+      END
