diff --git a/SRC/dlaqz0.f b/SRC/dlaqz0.f
index 5bfc78bcc..1714480f4 100644
--- a/SRC/dlaqz0.f
+++ b/SRC/dlaqz0.f
@@ -409,8 +409,8 @@
 
 *     rcost = ilaenv( 17,'DLAQZ0',jbcmpz,n,ilo,ihi,lwork )
       rcost = 10
-      itemp1 = int(nshifts/sqrt( 1+2*nshifts/(dble(rcost)/100*n) ))
-      itemp1 = ((k-1)/4)*4+4
+      itemp1 = int(nsr/sqrt( 1+2*nsr/(dble(rcost)/100*n) ))
+      itemp1 = ((itemp1-1)/4)*4+4
       nbr = nsr+itemp1
 
       if( n .lt. nmin ) then
@@ -442,7 +442,7 @@
          info =-19
       end if
       if( info.NE.0 ) then
-         call xerbla( 'DLAQZ0',-info )
+         call xerbla( 'DLAQZ0',info )
          return
       end if
 *
@@ -462,7 +462,7 @@
       istop = ihi
       maxit = 3*(ihi-ilo+1)
       ld = 0
-
+ 
       do iiter = 1,maxit
          if(iiter .ge. maxit) then
             info = istop
@@ -473,28 +473,30 @@
          end if
 
 *        Check deflations at the end
-         if (abs(A(istop-1,istop-2)) .le. ulp*(abs(A(istop-1,
-     $      istop-1))+abs(A(istop-2,istop-2)))) then
+         if (abs(A(istop-1,istop-2)) .le. max(smlnum,ulp*(abs(A(istop-1,
+     $      istop-1))+abs(A(istop-2,istop-2))))) then
             A(istop-1,istop-2) = zero
             istop = istop-2
             ld = 0
             eshift = zero
-         else if (abs(A(istop,istop-1)) .le. ulp*(abs(A(istop,
-     $      istop))+abs(A(istop-1,istop-1)))) then
+         else if (abs(A(istop,istop-1)) .le. max(smlnum,
+     $      ulp*(abs(A(istop,istop))+abs(A(istop-1,istop-1))))) then
             A(istop,istop-1) = zero
             istop = istop-1
             ld = 0
             eshift = zero
          end if
 *        Check deflations at the start
-         if (abs(A(istart+2,istart+1)) .le. ulp*(abs(A(istart+1,
-     $      istart+1))+abs(A(istart+2,istart+2)))) then
+         if (abs(A(istart+2,istart+1)) .le. max(smlnum,
+     $      ulp*(abs(A(istart+1,istart+1))+abs(A(istart+2,
+     $      istart+2))))) then
             A(istart+2,istart+1) = zero
             istart = istart+2
             ld = 0
             eshift = zero
-         else if (abs(A(istart+1,istart)) .le. ulp*(abs(A(istart,
-     $      istart))+abs(A(istart+1,istart+1)))) then
+         else if (abs(A(istart+1,istart)) .le. max(smlnum,
+     $      ulp*(abs(A(istart,istart))+abs(A(istart+1,
+     $      istart+1))))) then
             A(istart+1,istart) = zero
             istart = istart+1
             ld = 0
@@ -508,8 +510,8 @@
 *        Check interior deflations
          istart2 = istart
          do k = istop,istart+1,-1
-            if (abs(A(k,k-1)) .le. ulp*(abs(A(k,k))+abs(A(k-1,
-     $         k-1)))) then
+            if (abs(A(k,k-1)) .le. max(smlnum,ulp*(abs(A(k,k))+abs(A(k-
+     $         1,k-1))))) then
                A(k,k-1) = zero
                istart2 = k
                exit
@@ -671,7 +673,7 @@
 
       end do
 
-   80 call dlaqz3(ilschur,ilq,ilz,n,ilo,ihi,A,ldA,B,ldB,Q,ldQ,Z,ldZ,
+   80 call dlaqz3(ilschur,ilq,ilz,n,1,n,A,ldA,B,ldB,Q,ldQ,Z,ldZ,
      $alphar,alphai,beta,norm_info)
 
       end subroutine
diff --git a/SRC/dlaqz1.f b/SRC/dlaqz1.f
index 6406d5cdd..b7886f7b1 100644
--- a/SRC/dlaqz1.f
+++ b/SRC/dlaqz1.f
@@ -119,7 +119,13 @@
          parameter(zero=0.0d0,one=1.0d0,half=0.5d0)
 
 *        Local scalars
-         double precision :: w(2)
+         double precision :: w(2),safmin,safmax
+
+*        External Functions
+         double precision,external :: dlamch
+
+         safmin = dlamch('SAFE MINIMUM')
+         safmax = one/safmin
 
 *        Calculate first shifted vector
          w(1) = beta1*A(1,1)-sr1*B(1,1)
@@ -140,4 +146,12 @@
 *        Account for imaginary part
          v(1) = v(1)+si*si*B(1,1)
 
+         if( abs(v(1)).gt.safmax .or. abs(v(2)) .gt. safmax .or. abs(v(3
+     $      )).gt.safmax .or. isnan(v(1)) .or. isnan(v(2)) .or. isnan(v(
+     $      3)) ) then
+            v(1) = zero
+            v(2) = zero
+            v(3) = zero
+         end if
+
       end subroutine
diff --git a/SRC/dlaqz6.f b/SRC/dlaqz6.f
index 6d6731378..cf5525744 100644
--- a/SRC/dlaqz6.f
+++ b/SRC/dlaqz6.f
@@ -324,7 +324,7 @@
          info =-17
       end if
       if( info.NE.0 ) then
-         call xerbla( 'DLAQZ0',-info )
+         call xerbla( 'DLAQZ6',-info )
          return
       end if
    
@@ -343,22 +343,22 @@
       jbcmpz(2:2)=wantQ
       jbcmpz(3:3)=wantZ
 
-      nmin = ilaenv( 12,'DLAQZ0',jbcmpz,n,ilo,ihi,lwork )
+      nmin = ilaenv( 12,'DLAQZ6',jbcmpz,n,ilo,ihi,lwork )
 
-      nwr = ilaenv( 13,'DLAQZ0',jbcmpz,n,ilo,ihi,lwork )
+      nwr = ilaenv( 13,'DLAQZ6',jbcmpz,n,ilo,ihi,lwork )
       nwr = max( 2,nwr )
       nwr = min( ihi-ilo+1,( n-1 ) / 3,nwr )
 
-      nibble = ilaenv( 14,'DLAQZ0',jbcmpz,n,ilo,ihi,lwork )
+      nibble = ilaenv( 14,'DLAQZ6',jbcmpz,n,ilo,ihi,lwork )
       
-      nsr = ilaenv( 15,'DLAQZ0',jbcmpz,n,ilo,ihi,lwork )
+      nsr = ilaenv( 15,'DLAQZ6',jbcmpz,n,ilo,ihi,lwork )
       nsr = min( nsr,( n+6 ) / 9,ihi-ilo )
       nsr = max( 2,nsr-mod( nsr,2 ) )
 
 *     rcost = ilaenv( 17,'DLAQZ0',jbcmpz,n,ilo,ihi,lwork )
       rcost = 10
-      itemp1 = int(nshifts/sqrt( 1+2*nshifts/(dble(rcost)/100*n) ))
-      itemp1 = ((k-1)/4)*4+4
+      itemp1 = int(nsr/sqrt( 1+2*nsr/(dble(rcost)/100*n) ))
+      itemp1 = ((itemp1-1)/4)*4+4
       nbr = nsr+itemp1
 
       if( n .lt. nmin ) then
@@ -390,7 +390,7 @@
          info =-19
       end if
       if( info.NE.0 ) then
-         call xerbla( 'DLAQZ0',-info )
+         call xerbla( 'DLAQZ6',info )
          return
       end if
 *