arpack.xlahqr2.patch

diff -rupN arpack-ng-3.1.5/SRC/Makefile.am arpack-ng-cbb0bf599d53ea0ef5a5056f833dd04b2fae6b63/SRC/Makefile.am
--- arpack-ng-3.1.5/SRC/Makefile.am	2014-02-14 16:25:43.000000000 +0530
+++ arpack-ng-cbb0bf599d53ea0ef5a5056f833dd04b2fae6b63/SRC/Makefile.am	2014-09-15 04:10:16.000000000 +0530
@@ -7,10 +7,10 @@ libarpacksrc_la_SOURCES = \
 	dgetv0.f dlaqrb.f dstqrb.f dsortc.f dsortr.f dstatn.f dstats.f \
 	dnaitr.f dnapps.f dnaup2.f dnaupd.f dnconv.f dneigh.f dngets.f \
 	dsaitr.f dsapps.f dsaup2.f dsaupd.f dsconv.f dseigt.f dsgets.f \
-	dneupd.f dseupd.f dsesrt.f \
+	dneupd.f dseupd.f dsesrt.f dlahqr2.f slahqr2.f \
 	cnaitr.f cnapps.f cnaup2.f cnaupd.f cneigh.f cneupd.f cngets.f \
         cgetv0.f csortc.f cstatn.f \
 	znaitr.f znapps.f znaup2.f znaupd.f zneigh.f zneupd.f zngets.f \
         zgetv0.f zsortc.f zstatn.f
 
-EXTRA_DIST = debug.h stat.h version.h
\ No newline at end of file
+EXTRA_DIST = debug.h stat.h version.h
diff -rupN arpack-ng-3.1.5/SRC/Makefile.in arpack-ng-cbb0bf599d53ea0ef5a5056f833dd04b2fae6b63/SRC/Makefile.in
--- arpack-ng-3.1.5/SRC/Makefile.in	2014-02-15 19:16:02.000000000 +0530
+++ arpack-ng-cbb0bf599d53ea0ef5a5056f833dd04b2fae6b63/SRC/Makefile.in	2014-09-15 04:10:16.000000000 +0530
@@ -101,10 +101,10 @@ am_libarpacksrc_la_OBJECTS = sgetv0.lo s
 	dsortr.lo dstatn.lo dstats.lo dnaitr.lo dnapps.lo dnaup2.lo \
 	dnaupd.lo dnconv.lo dneigh.lo dngets.lo dsaitr.lo dsapps.lo \
 	dsaup2.lo dsaupd.lo dsconv.lo dseigt.lo dsgets.lo dneupd.lo \
-	dseupd.lo dsesrt.lo cnaitr.lo cnapps.lo cnaup2.lo cnaupd.lo \
-	cneigh.lo cneupd.lo cngets.lo cgetv0.lo csortc.lo cstatn.lo \
-	znaitr.lo znapps.lo znaup2.lo znaupd.lo zneigh.lo zneupd.lo \
-	zngets.lo zgetv0.lo zsortc.lo zstatn.lo
+	dseupd.lo dsesrt.lo dlahqr2.lo slahqr2.lo cnaitr.lo cnapps.lo \
+	cnaup2.lo cnaupd.lo cneigh.lo cneupd.lo cngets.lo cgetv0.lo \
+	csortc.lo cstatn.lo znaitr.lo znapps.lo znaup2.lo znaupd.lo \
+	zneigh.lo zneupd.lo zngets.lo zgetv0.lo zsortc.lo zstatn.lo
 libarpacksrc_la_OBJECTS = $(am_libarpacksrc_la_OBJECTS)
 AM_V_lt = $(am__v_lt_@AM_V@)
 am__v_lt_ = $(am__v_lt_@AM_DEFAULT_V@)
@@ -298,7 +298,7 @@ libarpacksrc_la_SOURCES = \
 	dgetv0.f dlaqrb.f dstqrb.f dsortc.f dsortr.f dstatn.f dstats.f \
 	dnaitr.f dnapps.f dnaup2.f dnaupd.f dnconv.f dneigh.f dngets.f \
 	dsaitr.f dsapps.f dsaup2.f dsaupd.f dsconv.f dseigt.f dsgets.f \
-	dneupd.f dseupd.f dsesrt.f \
+	dneupd.f dseupd.f dsesrt.f dlahqr2.f slahqr2.f \
 	cnaitr.f cnapps.f cnaup2.f cnaupd.f cneigh.f cneupd.f cngets.f \
         cgetv0.f csortc.f cstatn.f \
 	znaitr.f znapps.f znaup2.f znaupd.f zneigh.f zneupd.f zngets.f \
diff -rupN arpack-ng-3.1.5/SRC/dlahqr2.f arpack-ng-cbb0bf599d53ea0ef5a5056f833dd04b2fae6b63/SRC/dlahqr2.f
--- arpack-ng-3.1.5/SRC/dlahqr2.f	1970-01-01 05:30:00.000000000 +0530
+++ arpack-ng-cbb0bf599d53ea0ef5a5056f833dd04b2fae6b63/SRC/dlahqr2.f	2014-09-15 04:10:16.000000000 +0530
@@ -0,0 +1,410 @@
+      SUBROUTINE DLAHQR2( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
+     $                   ILOZ, IHIZ, Z, LDZ, INFO )
+*
+*  -- LAPACK auxiliary routine (version 2.0) --
+*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+*     Courant Institute, Argonne National Lab, and Rice University
+*     October 31, 1992
+*
+*     .. Scalar Arguments ..
+      LOGICAL            WANTT, WANTZ
+      INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DLAHQR2 is an auxiliary routine called by DHSEQR to update the
+*  eigenvalues and Schur decomposition already computed by DHSEQR, by
+*  dealing with the Hessenberg submatrix in rows and columns ILO to IHI.
+*
+*  Arguments
+*  =========
+*
+*  WANTT   (input) LOGICAL
+*          = .TRUE. : the full Schur form T is required;
+*          = .FALSE.: only eigenvalues are required.
+*
+*  WANTZ   (input) LOGICAL
+*          = .TRUE. : the matrix of Schur vectors Z is required;
+*          = .FALSE.: Schur vectors are not required.
+*
+*  N       (input) INTEGER
+*          The order of the matrix H.  N >= 0.
+*
+*  ILO     (input) INTEGER
+*  IHI     (input) INTEGER
+*          It is assumed that H is already upper quasi-triangular in
+*          rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless
+*          ILO = 1). DLAHQR works primarily with the Hessenberg
+*          submatrix in rows and columns ILO to IHI, but applies
+*          transformations to all of H if WANTT is .TRUE..
+*          1 <= ILO <= max(1,IHI); IHI <= N.
+*
+*  H       (input/output) DOUBLE PRECISION array, dimension (LDH,N)
+*          On entry, the upper Hessenberg matrix H.
+*          On exit, if WANTT is .TRUE., H is upper quasi-triangular in
+*          rows and columns ILO:IHI, with any 2-by-2 diagonal blocks in
+*          standard form. If WANTT is .FALSE., the contents of H are
+*          unspecified on exit.
+*
+*  LDH     (input) INTEGER
+*          The leading dimension of the array H. LDH >= max(1,N).
+*
+*  WR      (output) DOUBLE PRECISION array, dimension (N)
+*  WI      (output) DOUBLE PRECISION array, dimension (N)
+*          The real and imaginary parts, respectively, of the computed
+*          eigenvalues ILO to IHI are stored in the corresponding
+*          elements of WR and WI. If two eigenvalues are computed as a
+*          complex conjugate pair, they are stored in consecutive
+*          elements of WR and WI, say the i-th and (i+1)th, with
+*          WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the
+*          eigenvalues are stored in the same order as on the diagonal
+*          of the Schur form returned in H, with WR(i) = H(i,i), and, if
+*          H(i:i+1,i:i+1) is a 2-by-2 diagonal block,
+*          WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i).
+*
+*  ILOZ    (input) INTEGER
+*  IHIZ    (input) INTEGER
+*          Specify the rows of Z to which transformations must be
+*          applied if WANTZ is .TRUE..
+*          1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
+*
+*  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
+*          If WANTZ is .TRUE., on entry Z must contain the current
+*          matrix Z of transformations accumulated by DHSEQR, and on
+*          exit Z has been updated; transformations are applied only to
+*          the submatrix Z(ILOZ:IHIZ,ILO:IHI).
+*          If WANTZ is .FALSE., Z is not referenced.
+*
+*  LDZ     (input) INTEGER
+*          The leading dimension of the array Z. LDZ >= max(1,N).
+*
+*  INFO    (output) INTEGER
+*          = 0: successful exit
+*          > 0: DLAHQR failed to compute all the eigenvalues ILO to IHI
+*               in a total of 30*(IHI-ILO+1) iterations; if INFO = i,
+*               elements i+1:ihi of WR and WI contain those eigenvalues
+*               which have been successfully computed.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
+      DOUBLE PRECISION   DAT1, DAT2
+      PARAMETER          ( DAT1 = 0.75D+0, DAT2 = -0.4375D+0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            I, I1, I2, ITN, ITS, J, K, L, M, NH, NR, NZ
+      DOUBLE PRECISION   CS, H00, H10, H11, H12, H21, H22, H33, H33S,
+     $                   H43H34, H44, H44S, OVFL, S, SMLNUM, SN, SUM,
+     $                   T1, T2, T3, TST1, ULP, UNFL, V1, V2, V3
+*     ..
+*     .. Local Arrays ..
+      DOUBLE PRECISION   V( 3 ), WORK( 1 )
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   DLAMCH, DLANHS
+      EXTERNAL           DLAMCH, DLANHS
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           DCOPY, DLABAD, DLANV2, DLARFG, DROT
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, MAX, MIN
+*     ..
+*     .. Executable Statements ..
+*
+      INFO = 0
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+      IF( ILO.EQ.IHI ) THEN
+         WR( ILO ) = H( ILO, ILO )
+         WI( ILO ) = ZERO
+         RETURN
+      END IF
+*
+      NH = IHI - ILO + 1
+      NZ = IHIZ - ILOZ + 1
+*
+*     Set machine-dependent constants for the stopping criterion.
+*     If norm(H) <= sqrt(OVFL), overflow should not occur.
+*
+      UNFL = DLAMCH( 'Safe minimum' )
+      OVFL = ONE / UNFL
+      CALL DLABAD( UNFL, OVFL )
+      ULP = DLAMCH( 'Precision' )
+      SMLNUM = UNFL*( NH / ULP )
+*
+*     I1 and I2 are the indices of the first row and last column of H
+*     to which transformations must be applied. If eigenvalues only are
+*     being computed, I1 and I2 are set inside the main loop.
+*
+      IF( WANTT ) THEN
+         I1 = 1
+         I2 = N
+      END IF
+*
+*     ITN is the total number of QR iterations allowed.
+*
+      ITN = 30*NH
+*
+*     The main loop begins here. I is the loop index and decreases from
+*     IHI to ILO in steps of 1 or 2. Each iteration of the loop works
+*     with the active submatrix in rows and columns L to I.
+*     Eigenvalues I+1 to IHI have already converged. Either L = ILO or
+*     H(L,L-1) is negligible so that the matrix splits.
+*
+      I = IHI
+   10 CONTINUE
+      L = ILO
+      IF( I.LT.ILO )
+     $   GO TO 150
+*
+*     Perform QR iterations on rows and columns ILO to I until a
+*     submatrix of order 1 or 2 splits off at the bottom because a
+*     subdiagonal element has become negligible.
+*
+      DO 130 ITS = 0, ITN
+*
+*        Look for a single small subdiagonal element.
+*
+         DO 20 K = I, L + 1, -1
+            TST1 = ABS( H( K-1, K-1 ) ) + ABS( H( K, K ) )
+            IF( TST1.EQ.ZERO )
+     $         TST1 = DLANHS( '1', I-L+1, H( L, L ), LDH, WORK )
+            IF( ABS( H( K, K-1 ) ).LE.MAX( ULP*TST1, SMLNUM ) )
+     $         GO TO 30
+   20    CONTINUE
+   30    CONTINUE
+         L = K
+         IF( L.GT.ILO ) THEN
+*
+*           H(L,L-1) is negligible
+*
+            H( L, L-1 ) = ZERO
+         END IF
+*
+*        Exit from loop if a submatrix of order 1 or 2 has split off.
+*
+         IF( L.GE.I-1 )
+     $      GO TO 140
+*
+*        Now the active submatrix is in rows and columns L to I. If
+*        eigenvalues only are being computed, only the active submatrix
+*        need be transformed.
+*
+         IF( .NOT.WANTT ) THEN
+            I1 = L
+            I2 = I
+         END IF
+*
+         IF( ITS.EQ.10 .OR. ITS.EQ.20 ) THEN
+*
+*           Exceptional shift.
+*
+            S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) )
+            H44 = DAT1*S
+            H33 = H44
+            H43H34 = DAT2*S*S
+         ELSE
+*
+*           Prepare to use Wilkinson's double shift
+*
+            H44 = H( I, I )
+            H33 = H( I-1, I-1 )
+            H43H34 = H( I, I-1 )*H( I-1, I )
+         END IF
+*
+*        Look for two consecutive small subdiagonal elements.
+*
+         DO 40 M = I - 2, L, -1
+*
+*           Determine the effect of starting the double-shift QR
+*           iteration at row M, and see if this would make H(M,M-1)
+*           negligible.
+*
+            H11 = H( M, M )
+            H22 = H( M+1, M+1 )
+            H21 = H( M+1, M )
+            H12 = H( M, M+1 )
+            H44S = H44 - H11
+            H33S = H33 - H11
+            V1 = ( H33S*H44S-H43H34 ) / H21 + H12
+            V2 = H22 - H11 - H33S - H44S
+            V3 = H( M+2, M+1 )
+            S = ABS( V1 ) + ABS( V2 ) + ABS( V3 )
+            V1 = V1 / S
+            V2 = V2 / S
+            V3 = V3 / S
+            V( 1 ) = V1
+            V( 2 ) = V2
+            V( 3 ) = V3
+            IF( M.EQ.L )
+     $         GO TO 50
+            H00 = H( M-1, M-1 )
+            H10 = H( M, M-1 )
+            TST1 = ABS( V1 )*( ABS( H00 )+ABS( H11 )+ABS( H22 ) )
+            IF( ABS( H10 )*( ABS( V2 )+ABS( V3 ) ).LE.ULP*TST1 )
+     $         GO TO 50
+   40    CONTINUE
+   50    CONTINUE
+*
+*        Double-shift QR step
+*
+         DO 120 K = M, I - 1
+*
+*           The first iteration of this loop determines a reflection G
+*           from the vector V and applies it from left and right to H,
+*           thus creating a nonzero bulge below the subdiagonal.
+*
+*           Each subsequent iteration determines a reflection G to
+*           restore the Hessenberg form in the (K-1)th column, and thus
+*           chases the bulge one step toward the bottom of the active
+*           submatrix. NR is the order of G.
+*
+            NR = MIN( 3, I-K+1 )
+            IF( K.GT.M )
+     $         CALL DCOPY( NR, H( K, K-1 ), 1, V, 1 )
+            CALL DLARFG( NR, V( 1 ), V( 2 ), 1, T1 )
+            IF( K.GT.M ) THEN
+               H( K, K-1 ) = V( 1 )
+               H( K+1, K-1 ) = ZERO
+               IF( K.LT.I-1 )
+     $            H( K+2, K-1 ) = ZERO
+            ELSE IF( M.GT.L ) THEN
+               H( K, K-1 ) = -H( K, K-1 )
+            END IF
+            V2 = V( 2 )
+            T2 = T1*V2
+            IF( NR.EQ.3 ) THEN
+               V3 = V( 3 )
+               T3 = T1*V3
+*
+*              Apply G from the left to transform the rows of the matrix
+*              in columns K to I2.
+*
+               DO 60 J = K, I2
+                  SUM = H( K, J ) + V2*H( K+1, J ) + V3*H( K+2, J )
+                  H( K, J ) = H( K, J ) - SUM*T1
+                  H( K+1, J ) = H( K+1, J ) - SUM*T2
+                  H( K+2, J ) = H( K+2, J ) - SUM*T3
+   60          CONTINUE
+*
+*              Apply G from the right to transform the columns of the
+*              matrix in rows I1 to min(K+3,I).
+*
+               DO 70 J = I1, MIN( K+3, I )
+                  SUM = H( J, K ) + V2*H( J, K+1 ) + V3*H( J, K+2 )
+                  H( J, K ) = H( J, K ) - SUM*T1
+                  H( J, K+1 ) = H( J, K+1 ) - SUM*T2
+                  H( J, K+2 ) = H( J, K+2 ) - SUM*T3
+   70          CONTINUE
+*
+               IF( WANTZ ) THEN
+*
+*                 Accumulate transformations in the matrix Z
+*
+                  DO 80 J = ILOZ, IHIZ
+                     SUM = Z( J, K ) + V2*Z( J, K+1 ) + V3*Z( J, K+2 )
+                     Z( J, K ) = Z( J, K ) - SUM*T1
+                     Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2
+                     Z( J, K+2 ) = Z( J, K+2 ) - SUM*T3
+   80             CONTINUE
+               END IF
+            ELSE IF( NR.EQ.2 ) THEN
+*
+*              Apply G from the left to transform the rows of the matrix
+*              in columns K to I2.
+*
+               DO 90 J = K, I2
+                  SUM = H( K, J ) + V2*H( K+1, J )
+                  H( K, J ) = H( K, J ) - SUM*T1
+                  H( K+1, J ) = H( K+1, J ) - SUM*T2
+   90          CONTINUE
+*
+*              Apply G from the right to transform the columns of the
+*              matrix in rows I1 to min(K+3,I).
+*
+               DO 100 J = I1, I
+                  SUM = H( J, K ) + V2*H( J, K+1 )
+                  H( J, K ) = H( J, K ) - SUM*T1
+                  H( J, K+1 ) = H( J, K+1 ) - SUM*T2
+  100          CONTINUE
+*
+               IF( WANTZ ) THEN
+*
+*                 Accumulate transformations in the matrix Z
+*
+                  DO 110 J = ILOZ, IHIZ
+                     SUM = Z( J, K ) + V2*Z( J, K+1 )
+                     Z( J, K ) = Z( J, K ) - SUM*T1
+                     Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2
+  110             CONTINUE
+               END IF
+            END IF
+  120    CONTINUE
+*
+  130 CONTINUE
+*
+*     Failure to converge in remaining number of iterations
+*
+      INFO = I
+      RETURN
+*
+  140 CONTINUE
+*
+      IF( L.EQ.I ) THEN
+*
+*        H(I,I-1) is negligible: one eigenvalue has converged.
+*
+         WR( I ) = H( I, I )
+         WI( I ) = ZERO
+      ELSE IF( L.EQ.I-1 ) THEN
+*
+*        H(I-1,I-2) is negligible: a pair of eigenvalues have converged.
+*
+*        Transform the 2-by-2 submatrix to standard Schur form,
+*        and compute and store the eigenvalues.
+*
+         CALL DLANV2( H( I-1, I-1 ), H( I-1, I ), H( I, I-1 ),
+     $                H( I, I ), WR( I-1 ), WI( I-1 ), WR( I ), WI( I ),
+     $                CS, SN )
+*
+         IF( WANTT ) THEN
+*
+*           Apply the transformation to the rest of H.
+*
+            IF( I2.GT.I )
+     $         CALL DROT( I2-I, H( I-1, I+1 ), LDH, H( I, I+1 ), LDH,
+     $                    CS, SN )
+            CALL DROT( I-I1-1, H( I1, I-1 ), 1, H( I1, I ), 1, CS, SN )
+         END IF
+         IF( WANTZ ) THEN
+*
+*           Apply the transformation to Z.
+*
+            CALL DROT( NZ, Z( ILOZ, I-1 ), 1, Z( ILOZ, I ), 1, CS, SN )
+         END IF
+      END IF
+*
+*     Decrement number of remaining iterations, and return to start of
+*     the main loop with new value of I.
+*
+      ITN = ITN - ITS
+      I = L - 1
+      GO TO 10
+*
+  150 CONTINUE
+      RETURN
+*
+*     End of DLAHQR
+*
+      END
diff -rupN arpack-ng-3.1.5/SRC/dneupd.f arpack-ng-cbb0bf599d53ea0ef5a5056f833dd04b2fae6b63/SRC/dneupd.f
--- arpack-ng-3.1.5/SRC/dneupd.f	2014-02-14 16:25:43.000000000 +0530
+++ arpack-ng-cbb0bf599d53ea0ef5a5056f833dd04b2fae6b63/SRC/dneupd.f	2014-09-15 04:10:16.000000000 +0530
@@ -181,7 +181,7 @@ c          Error flag on output.
 c
 c          =  0: Normal exit.
 c
-c          =  1: The Schur form computed by LAPACK routine dlahqr 
+c          =  1: The Schur form computed by LAPACK routine dlahqr2 
 c                could not be reordered by LAPACK routine dtrsen .
 c                Re-enter subroutine dneupd  with IPARAM(5)=NCV and 
 c                increase the size of the arrays DR and DI to have 
@@ -197,7 +197,7 @@ c          = -5: WHICH must be one of 'L
 c          = -6: BMAT must be one of 'I' or 'G'.
 c          = -7: Length of private work WORKL array is not sufficient.
 c          = -8: Error return from calculation of a real Schur form.
-c                Informational error from LAPACK routine dlahqr .
+c                Informational error from LAPACK routine dlahqr2 .
 c          = -9: Error return from calculation of eigenvectors.
 c                Informational error from LAPACK routine dtrevc .
 c          = -10: IPARAM(7) must be 1,2,3,4.
@@ -232,7 +232,7 @@ c     dvout    ARPACK utility routine th
 c     dgeqr2   LAPACK routine that computes the QR factorization of 
 c             a matrix.
 c     dlacpy   LAPACK matrix copy routine.
-c     dlahqr   LAPACK routine to compute the real Schur form of an
+c     dlahqr2   LAPACK routine to compute the real Schur form of an
 c             upper Hessenberg matrix.
 c     dlamch   LAPACK routine that determines machine constants.
 c     dlapy2   LAPACK routine to compute sqrt(x**2+y**2) carefully.
@@ -364,7 +364,7 @@ c     | External Subroutines |
 c     %----------------------%
 c
       external   dcopy  , dger   , dgeqr2 , dlacpy , 
-     &           dlahqr , dlaset , dmout  , dorm2r , 
+     &           dlahqr2 , dlaset , dmout  , dorm2r , 
      &           dtrevc , dtrmm  , dtrsen , dscal  , 
      &           dvout  , ivout
 c
@@ -612,18 +612,18 @@ c
             go to 9000
          end if
 c
-c        %-----------------------------------------------------------%
-c        | Call LAPACK routine dlahqr  to compute the real Schur form |
-c        | of the upper Hessenberg matrix returned by DNAUPD .        |
-c        | Make a copy of the upper Hessenberg matrix.               |
-c        | Initialize the Schur vector matrix Q to the identity.     |
-c        %-----------------------------------------------------------%
+c        %-------------------------------------------------------------%
+c        | Call LAPACK routine dlahqr2  to compute the real Schur form |
+c        | of the upper Hessenberg matrix returned by DNAUPD .         |
+c        | Make a copy of the upper Hessenberg matrix.                 |
+c        | Initialize the Schur vector matrix Q to the identity.       |
+c        %-------------------------------------------------------------%
 c     
          call dcopy (ldh*ncv, workl(ih), 1, workl(iuptri), 1)
          call dlaset ('All', ncv, ncv, 
      &                zero , one, workl(invsub),
      &                ldq)
-         call dlahqr (.true., .true.       , ncv, 
+         call dlahqr2 (.true., .true.       , ncv,
      &                1     , ncv          , workl(iuptri), 
      &                ldh   , workl(iheigr), workl(iheigi),
      &                1     , ncv          , workl(invsub), 
diff -rupN arpack-ng-3.1.5/SRC/slahqr2.f arpack-ng-cbb0bf599d53ea0ef5a5056f833dd04b2fae6b63/SRC/slahqr2.f
--- arpack-ng-3.1.5/SRC/slahqr2.f	1970-01-01 05:30:00.000000000 +0530
+++ arpack-ng-cbb0bf599d53ea0ef5a5056f833dd04b2fae6b63/SRC/slahqr2.f	2014-09-15 04:10:16.000000000 +0530
@@ -0,0 +1,410 @@
+      SUBROUTINE SLAHQR2( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
+     $                   ILOZ, IHIZ, Z, LDZ, INFO )
+*
+*  -- LAPACK auxiliary routine (version 2.0) --
+*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+*     Courant Institute, Argonne National Lab, and Rice University
+*     October 31, 1992
+*
+*     .. Scalar Arguments ..
+      LOGICAL            WANTT, WANTZ
+      INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
+*     ..
+*     .. Array Arguments ..
+      REAL               H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  SLAHQR2 is an auxiliary routine called by SHSEQR to update the
+*  eigenvalues and Schur decomposition already computed by SHSEQR, by
+*  dealing with the Hessenberg submatrix in rows and columns ILO to IHI.
+*
+*  Arguments
+*  =========
+*
+*  WANTT   (input) LOGICAL
+*          = .TRUE. : the full Schur form T is required;
+*          = .FALSE.: only eigenvalues are required.
+*
+*  WANTZ   (input) LOGICAL
+*          = .TRUE. : the matrix of Schur vectors Z is required;
+*          = .FALSE.: Schur vectors are not required.
+*
+*  N       (input) INTEGER
+*          The order of the matrix H.  N >= 0.
+*
+*  ILO     (input) INTEGER
+*  IHI     (input) INTEGER
+*          It is assumed that H is already upper quasi-triangular in
+*          rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless
+*          ILO = 1). SLAHQR works primarily with the Hessenberg
+*          submatrix in rows and columns ILO to IHI, but applies
+*          transformations to all of H if WANTT is .TRUE..
+*          1 <= ILO <= max(1,IHI); IHI <= N.
+*
+*  H       (input/output) REAL array, dimension (LDH,N)
+*          On entry, the upper Hessenberg matrix H.
+*          On exit, if WANTT is .TRUE., H is upper quasi-triangular in
+*          rows and columns ILO:IHI, with any 2-by-2 diagonal blocks in
+*          standard form. If WANTT is .FALSE., the contents of H are
+*          unspecified on exit.
+*
+*  LDH     (input) INTEGER
+*          The leading dimension of the array H. LDH >= max(1,N).
+*
+*  WR      (output) REAL array, dimension (N)
+*  WI      (output) REAL array, dimension (N)
+*          The real and imaginary parts, respectively, of the computed
+*          eigenvalues ILO to IHI are stored in the corresponding
+*          elements of WR and WI. If two eigenvalues are computed as a
+*          complex conjugate pair, they are stored in consecutive
+*          elements of WR and WI, say the i-th and (i+1)th, with
+*          WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the
+*          eigenvalues are stored in the same order as on the diagonal
+*          of the Schur form returned in H, with WR(i) = H(i,i), and, if
+*          H(i:i+1,i:i+1) is a 2-by-2 diagonal block,
+*          WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i).
+*
+*  ILOZ    (input) INTEGER
+*  IHIZ    (input) INTEGER
+*          Specify the rows of Z to which transformations must be
+*          applied if WANTZ is .TRUE..
+*          1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
+*
+*  Z       (input/output) REAL array, dimension (LDZ,N)
+*          If WANTZ is .TRUE., on entry Z must contain the current
+*          matrix Z of transformations accumulated by SHSEQR, and on
+*          exit Z has been updated; transformations are applied only to
+*          the submatrix Z(ILOZ:IHIZ,ILO:IHI).
+*          If WANTZ is .FALSE., Z is not referenced.
+*
+*  LDZ     (input) INTEGER
+*          The leading dimension of the array Z. LDZ >= max(1,N).
+*
+*  INFO    (output) INTEGER
+*          = 0: successful exit
+*          > 0: SLAHQR failed to compute all the eigenvalues ILO to IHI
+*               in a total of 30*(IHI-ILO+1) iterations; if INFO = i,
+*               elements i+1:ihi of WR and WI contain those eigenvalues
+*               which have been successfully computed.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ZERO, ONE
+      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
+      REAL               DAT1, DAT2
+      PARAMETER          ( DAT1 = 0.75E+0, DAT2 = -0.4375E+0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            I, I1, I2, ITN, ITS, J, K, L, M, NH, NR, NZ
+      REAL               CS, H00, H10, H11, H12, H21, H22, H33, H33S,
+     $                   H43H34, H44, H44S, OVFL, S, SMLNUM, SN, SUM,
+     $                   T1, T2, T3, TST1, ULP, UNFL, V1, V2, V3
+*     ..
+*     .. Local Arrays ..
+      REAL               V( 3 ), WORK( 1 )
+*     ..
+*     .. External Functions ..
+      REAL               SLAMCH, SLANHS
+      EXTERNAL           SLAMCH, SLANHS
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           SCOPY, SLABAD, SLANV2, SLARFG, SROT
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, MAX, MIN
+*     ..
+*     .. Executable Statements ..
+*
+      INFO = 0
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+      IF( ILO.EQ.IHI ) THEN
+         WR( ILO ) = H( ILO, ILO )
+         WI( ILO ) = ZERO
+         RETURN
+      END IF
+*
+      NH = IHI - ILO + 1
+      NZ = IHIZ - ILOZ + 1
+*
+*     Set machine-dependent constants for the stopping criterion.
+*     If norm(H) <= sqrt(OVFL), overflow should not occur.
+*
+      UNFL = SLAMCH( 'Safe minimum' )
+      OVFL = ONE / UNFL
+      CALL SLABAD( UNFL, OVFL )
+      ULP = SLAMCH( 'Precision' )
+      SMLNUM = UNFL*( NH / ULP )
+*
+*     I1 and I2 are the indices of the first row and last column of H
+*     to which transformations must be applied. If eigenvalues only are
+*     being computed, I1 and I2 are set inside the main loop.
+*
+      IF( WANTT ) THEN
+         I1 = 1
+         I2 = N
+      END IF
+*
+*     ITN is the total number of QR iterations allowed.
+*
+      ITN = 30*NH
+*
+*     The main loop begins here. I is the loop index and decreases from
+*     IHI to ILO in steps of 1 or 2. Each iteration of the loop works
+*     with the active submatrix in rows and columns L to I.
+*     Eigenvalues I+1 to IHI have already converged. Either L = ILO or
+*     H(L,L-1) is negligible so that the matrix splits.
+*
+      I = IHI
+   10 CONTINUE
+      L = ILO
+      IF( I.LT.ILO )
+     $   GO TO 150
+*
+*     Perform QR iterations on rows and columns ILO to I until a
+*     submatrix of order 1 or 2 splits off at the bottom because a
+*     subdiagonal element has become negligible.
+*
+      DO 130 ITS = 0, ITN
+*
+*        Look for a single small subdiagonal element.
+*
+         DO 20 K = I, L + 1, -1
+            TST1 = ABS( H( K-1, K-1 ) ) + ABS( H( K, K ) )
+            IF( TST1.EQ.ZERO )
+     $         TST1 = SLANHS( '1', I-L+1, H( L, L ), LDH, WORK )
+            IF( ABS( H( K, K-1 ) ).LE.MAX( ULP*TST1, SMLNUM ) )
+     $         GO TO 30
+   20    CONTINUE
+   30    CONTINUE
+         L = K
+         IF( L.GT.ILO ) THEN
+*
+*           H(L,L-1) is negligible
+*
+            H( L, L-1 ) = ZERO
+         END IF
+*
+*        Exit from loop if a submatrix of order 1 or 2 has split off.
+*
+         IF( L.GE.I-1 )
+     $      GO TO 140
+*
+*        Now the active submatrix is in rows and columns L to I. If
+*        eigenvalues only are being computed, only the active submatrix
+*        need be transformed.
+*
+         IF( .NOT.WANTT ) THEN
+            I1 = L
+            I2 = I
+         END IF
+*
+         IF( ITS.EQ.10 .OR. ITS.EQ.20 ) THEN
+*
+*           Exceptional shift.
+*
+            S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) )
+            H44 = DAT1*S
+            H33 = H44
+            H43H34 = DAT2*S*S
+         ELSE
+*
+*           Prepare to use Wilkinson's double shift
+*
+            H44 = H( I, I )
+            H33 = H( I-1, I-1 )
+            H43H34 = H( I, I-1 )*H( I-1, I )
+         END IF
+*
+*        Look for two consecutive small subdiagonal elements.
+*
+         DO 40 M = I - 2, L, -1
+*
+*           Determine the effect of starting the double-shift QR
+*           iteration at row M, and see if this would make H(M,M-1)
+*           negligible.
+*
+            H11 = H( M, M )
+            H22 = H( M+1, M+1 )
+            H21 = H( M+1, M )
+            H12 = H( M, M+1 )
+            H44S = H44 - H11
+            H33S = H33 - H11
+            V1 = ( H33S*H44S-H43H34 ) / H21 + H12
+            V2 = H22 - H11 - H33S - H44S
+            V3 = H( M+2, M+1 )
+            S = ABS( V1 ) + ABS( V2 ) + ABS( V3 )
+            V1 = V1 / S
+            V2 = V2 / S
+            V3 = V3 / S
+            V( 1 ) = V1
+            V( 2 ) = V2
+            V( 3 ) = V3
+            IF( M.EQ.L )
+     $         GO TO 50
+            H00 = H( M-1, M-1 )
+            H10 = H( M, M-1 )
+            TST1 = ABS( V1 )*( ABS( H00 )+ABS( H11 )+ABS( H22 ) )
+            IF( ABS( H10 )*( ABS( V2 )+ABS( V3 ) ).LE.ULP*TST1 )
+     $         GO TO 50
+   40    CONTINUE
+   50    CONTINUE
+*
+*        Double-shift QR step
+*
+         DO 120 K = M, I - 1
+*
+*           The first iteration of this loop determines a reflection G
+*           from the vector V and applies it from left and right to H,
+*           thus creating a nonzero bulge below the subdiagonal.
+*
+*           Each subsequent iteration determines a reflection G to
+*           restore the Hessenberg form in the (K-1)th column, and thus
+*           chases the bulge one step toward the bottom of the active
+*           submatrix. NR is the order of G.
+*
+            NR = MIN( 3, I-K+1 )
+            IF( K.GT.M )
+     $         CALL SCOPY( NR, H( K, K-1 ), 1, V, 1 )
+            CALL SLARFG( NR, V( 1 ), V( 2 ), 1, T1 )
+            IF( K.GT.M ) THEN
+               H( K, K-1 ) = V( 1 )
+               H( K+1, K-1 ) = ZERO
+               IF( K.LT.I-1 )
+     $            H( K+2, K-1 ) = ZERO
+            ELSE IF( M.GT.L ) THEN
+               H( K, K-1 ) = -H( K, K-1 )
+            END IF
+            V2 = V( 2 )
+            T2 = T1*V2
+            IF( NR.EQ.3 ) THEN
+               V3 = V( 3 )
+               T3 = T1*V3
+*
+*              Apply G from the left to transform the rows of the matrix
+*              in columns K to I2.
+*
+               DO 60 J = K, I2
+                  SUM = H( K, J ) + V2*H( K+1, J ) + V3*H( K+2, J )
+                  H( K, J ) = H( K, J ) - SUM*T1
+                  H( K+1, J ) = H( K+1, J ) - SUM*T2
+                  H( K+2, J ) = H( K+2, J ) - SUM*T3
+   60          CONTINUE
+*
+*              Apply G from the right to transform the columns of the
+*              matrix in rows I1 to min(K+3,I).
+*
+               DO 70 J = I1, MIN( K+3, I )
+                  SUM = H( J, K ) + V2*H( J, K+1 ) + V3*H( J, K+2 )
+                  H( J, K ) = H( J, K ) - SUM*T1
+                  H( J, K+1 ) = H( J, K+1 ) - SUM*T2
+                  H( J, K+2 ) = H( J, K+2 ) - SUM*T3
+   70          CONTINUE
+*
+               IF( WANTZ ) THEN
+*
+*                 Accumulate transformations in the matrix Z
+*
+                  DO 80 J = ILOZ, IHIZ
+                     SUM = Z( J, K ) + V2*Z( J, K+1 ) + V3*Z( J, K+2 )
+                     Z( J, K ) = Z( J, K ) - SUM*T1
+                     Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2
+                     Z( J, K+2 ) = Z( J, K+2 ) - SUM*T3
+   80             CONTINUE
+               END IF
+            ELSE IF( NR.EQ.2 ) THEN
+*
+*              Apply G from the left to transform the rows of the matrix
+*              in columns K to I2.
+*
+               DO 90 J = K, I2
+                  SUM = H( K, J ) + V2*H( K+1, J )
+                  H( K, J ) = H( K, J ) - SUM*T1
+                  H( K+1, J ) = H( K+1, J ) - SUM*T2
+   90          CONTINUE
+*
+*              Apply G from the right to transform the columns of the
+*              matrix in rows I1 to min(K+3,I).
+*
+               DO 100 J = I1, I
+                  SUM = H( J, K ) + V2*H( J, K+1 )
+                  H( J, K ) = H( J, K ) - SUM*T1
+                  H( J, K+1 ) = H( J, K+1 ) - SUM*T2
+  100          CONTINUE
+*
+               IF( WANTZ ) THEN
+*
+*                 Accumulate transformations in the matrix Z
+*
+                  DO 110 J = ILOZ, IHIZ
+                     SUM = Z( J, K ) + V2*Z( J, K+1 )
+                     Z( J, K ) = Z( J, K ) - SUM*T1
+                     Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2
+  110             CONTINUE
+               END IF
+            END IF
+  120    CONTINUE
+*
+  130 CONTINUE
+*
+*     Failure to converge in remaining number of iterations
+*
+      INFO = I
+      RETURN
+*
+  140 CONTINUE
+*
+      IF( L.EQ.I ) THEN
+*
+*        H(I,I-1) is negligible: one eigenvalue has converged.
+*
+         WR( I ) = H( I, I )
+         WI( I ) = ZERO
+      ELSE IF( L.EQ.I-1 ) THEN
+*
+*        H(I-1,I-2) is negligible: a pair of eigenvalues have converged.
+*
+*        Transform the 2-by-2 submatrix to standard Schur form,
+*        and compute and store the eigenvalues.
+*
+         CALL SLANV2( H( I-1, I-1 ), H( I-1, I ), H( I, I-1 ),
+     $                H( I, I ), WR( I-1 ), WI( I-1 ), WR( I ), WI( I ),
+     $                CS, SN )
+*
+         IF( WANTT ) THEN
+*
+*           Apply the transformation to the rest of H.
+*
+            IF( I2.GT.I )
+     $         CALL SROT( I2-I, H( I-1, I+1 ), LDH, H( I, I+1 ), LDH,
+     $                    CS, SN )
+            CALL SROT( I-I1-1, H( I1, I-1 ), 1, H( I1, I ), 1, CS, SN )
+         END IF
+         IF( WANTZ ) THEN
+*
+*           Apply the transformation to Z.
+*
+            CALL SROT( NZ, Z( ILOZ, I-1 ), 1, Z( ILOZ, I ), 1, CS, SN )
+         END IF
+      END IF
+*
+*     Decrement number of remaining iterations, and return to start of
+*     the main loop with new value of I.
+*
+      ITN = ITN - ITS
+      I = L - 1
+      GO TO 10
+*
+  150 CONTINUE
+      RETURN
+*
+*     End of SLAHQR
+*
+      END
diff -rupN arpack-ng-3.1.5/SRC/sneupd.f arpack-ng-cbb0bf599d53ea0ef5a5056f833dd04b2fae6b63/SRC/sneupd.f
--- arpack-ng-3.1.5/SRC/sneupd.f	2014-02-14 16:25:43.000000000 +0530
+++ arpack-ng-cbb0bf599d53ea0ef5a5056f833dd04b2fae6b63/SRC/sneupd.f	2014-09-15 04:10:16.000000000 +0530
@@ -181,7 +181,7 @@ c          Error flag on output.
 c
 c          =  0: Normal exit.
 c
-c          =  1: The Schur form computed by LAPACK routine slahqr
+c          =  1: The Schur form computed by LAPACK routine slahqr2
 c                could not be reordered by LAPACK routine strsen.
 c                Re-enter subroutine sneupd with IPARAM(5)=NCV and 
 c                increase the size of the arrays DR and DI to have 
@@ -197,7 +197,7 @@ c          = -5: WHICH must be one of 'L
 c          = -6: BMAT must be one of 'I' or 'G'.
 c          = -7: Length of private work WORKL array is not sufficient.
 c          = -8: Error return from calculation of a real Schur form.
-c                Informational error from LAPACK routine slahqr.
+c                Informational error from LAPACK routine slahqr2.
 c          = -9: Error return from calculation of eigenvectors.
 c                Informational error from LAPACK routine strevc.
 c          = -10: IPARAM(7) must be 1,2,3,4.
@@ -232,7 +232,7 @@ c     svout   ARPACK utility routine tha
 c     sgeqr2  LAPACK routine that computes the QR factorization of 
 c             a matrix.
 c     slacpy  LAPACK matrix copy routine.
-c     slahqr  LAPACK routine to compute the real Schur form of an
+c     slahqr2  LAPACK routine to compute the real Schur form of an
 c             upper Hessenberg matrix.
 c     slamch  LAPACK routine that determines machine constants.
 c     slapy2  LAPACK routine to compute sqrt(x**2+y**2) carefully.
@@ -364,7 +364,7 @@ c     | External Subroutines |
 c     %----------------------%
 c
       external   scopy , sger  , sgeqr2, slacpy, 
-     &           slahqr, slaset, smout , sorm2r, 
+     &           slahqr2, slaset, smout , sorm2r, 
      &           strevc, strmm , strsen, sscal , 
      &           svout , ivout
 c
@@ -613,7 +613,7 @@ c
          end if
 c
 c        %-----------------------------------------------------------%
-c        | Call LAPACK routine slahqr to compute the real Schur form |
+c        | Call LAPACK routine slahqr2 to compute the real Schur form |
 c        | of the upper Hessenberg matrix returned by SNAUPD.        |
 c        | Make a copy of the upper Hessenberg matrix.               |
 c        | Initialize the Schur vector matrix Q to the identity.     |
@@ -623,7 +623,7 @@ c     
          call slaset('All', ncv, ncv, 
      &                zero , one, workl(invsub),
      &                ldq)
-         call slahqr(.true., .true.       , ncv, 
+         call slahqr2(.true., .true.       , ncv, 
      &                1     , ncv          , workl(iuptri), 
      &                ldh   , workl(iheigr), workl(iheigi),
      &                1     , ncv          , workl(invsub),