SRC/clahqr.f

   1 *> \brief \b CLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CLAHQR + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clahqr.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clahqr.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clahqr.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
  22 *                          IHIZ, Z, LDZ, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
  26 *       LOGICAL            WANTT, WANTZ
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       COMPLEX            H( LDH, * ), W( * ), Z( LDZ, * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *>    CLAHQR is an auxiliary routine called by CHSEQR to update the
  39 *>    eigenvalues and Schur decomposition already computed by CHSEQR, by
  40 *>    dealing with the Hessenberg submatrix in rows and columns ILO to
  41 *>    IHI.
  42 *> \endverbatim
  43 *
  44 *  Arguments:
  45 *  ==========
  46 *
  47 *> \param[in] WANTT
  48 *> \verbatim
  49 *>          WANTT is LOGICAL
  50 *>          = .TRUE. : the full Schur form T is required;
  51 *>          = .FALSE.: only eigenvalues are required.
  52 *> \endverbatim
  53 *>
  54 *> \param[in] WANTZ
  55 *> \verbatim
  56 *>          WANTZ is LOGICAL
  57 *>          = .TRUE. : the matrix of Schur vectors Z is required;
  58 *>          = .FALSE.: Schur vectors are not required.
  59 *> \endverbatim
  60 *>
  61 *> \param[in] N
  62 *> \verbatim
  63 *>          N is INTEGER
  64 *>          The order of the matrix H.  N >= 0.
  65 *> \endverbatim
  66 *>
  67 *> \param[in] ILO
  68 *> \verbatim
  69 *>          ILO is INTEGER
  70 *> \endverbatim
  71 *>
  72 *> \param[in] IHI
  73 *> \verbatim
  74 *>          IHI is INTEGER
  75 *>          It is assumed that H is already upper triangular in rows and
  76 *>          columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1).
  77 *>          CLAHQR works primarily with the Hessenberg submatrix in rows
  78 *>          and columns ILO to IHI, but applies transformations to all of
  79 *>          H if WANTT is .TRUE..
  80 *>          1 <= ILO <= max(1,IHI); IHI <= N.
  81 *> \endverbatim
  82 *>
  83 *> \param[in,out] H
  84 *> \verbatim
  85 *>          H is COMPLEX array, dimension (LDH,N)
  86 *>          On entry, the upper Hessenberg matrix H.
  87 *>          On exit, if INFO is zero and if WANTT is .TRUE., then H
  88 *>          is upper triangular in rows and columns ILO:IHI.  If INFO
  89 *>          is zero and if WANTT is .FALSE., then the contents of H
  90 *>          are unspecified on exit.  The output state of H in case
  91 *>          INF is positive is below under the description of INFO.
  92 *> \endverbatim
  93 *>
  94 *> \param[in] LDH
  95 *> \verbatim
  96 *>          LDH is INTEGER
  97 *>          The leading dimension of the array H. LDH >= max(1,N).
  98 *> \endverbatim
  99 *>
 100 *> \param[out] W
 101 *> \verbatim
 102 *>          W is COMPLEX array, dimension (N)
 103 *>          The computed eigenvalues ILO to IHI are stored in the
 104 *>          corresponding elements of W. If WANTT is .TRUE., the
 105 *>          eigenvalues are stored in the same order as on the diagonal
 106 *>          of the Schur form returned in H, with W(i) = H(i,i).
 107 *> \endverbatim
 108 *>
 109 *> \param[in] ILOZ
 110 *> \verbatim
 111 *>          ILOZ is INTEGER
 112 *> \endverbatim
 113 *>
 114 *> \param[in] IHIZ
 115 *> \verbatim
 116 *>          IHIZ is INTEGER
 117 *>          Specify the rows of Z to which transformations must be
 118 *>          applied if WANTZ is .TRUE..
 119 *>          1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
 120 *> \endverbatim
 121 *>
 122 *> \param[in,out] Z
 123 *> \verbatim
 124 *>          Z is COMPLEX array, dimension (LDZ,N)
 125 *>          If WANTZ is .TRUE., on entry Z must contain the current
 126 *>          matrix Z of transformations accumulated by CHSEQR, and on
 127 *>          exit Z has been updated; transformations are applied only to
 128 *>          the submatrix Z(ILOZ:IHIZ,ILO:IHI).
 129 *>          If WANTZ is .FALSE., Z is not referenced.
 130 *> \endverbatim
 131 *>
 132 *> \param[in] LDZ
 133 *> \verbatim
 134 *>          LDZ is INTEGER
 135 *>          The leading dimension of the array Z. LDZ >= max(1,N).
 136 *> \endverbatim
 137 *>
 138 *> \param[out] INFO
 139 *> \verbatim
 140 *>          INFO is INTEGER
 141 *>           =   0: successful exit
 142 *>          .GT. 0: if INFO = i, CLAHQR failed to compute all the
 143 *>                  eigenvalues ILO to IHI in a total of 30 iterations
 144 *>                  per eigenvalue; elements i+1:ihi of W contain
 145 *>                  those eigenvalues which have been successfully
 146 *>                  computed.
 147 *>
 148 *>                  If INFO .GT. 0 and WANTT is .FALSE., then on exit,
 149 *>                  the remaining unconverged eigenvalues are the
 150 *>                  eigenvalues of the upper Hessenberg matrix
 151 *>                  rows and columns ILO thorugh INFO of the final,
 152 *>                  output value of H.
 153 *>
 154 *>                  If INFO .GT. 0 and WANTT is .TRUE., then on exit
 155 *>          (*)       (initial value of H)*U  = U*(final value of H)
 156 *>                  where U is an orthognal matrix.    The final
 157 *>                  value of H is upper Hessenberg and triangular in
 158 *>                  rows and columns INFO+1 through IHI.
 159 *>
 160 *>                  If INFO .GT. 0 and WANTZ is .TRUE., then on exit
 161 *>                      (final value of Z)  = (initial value of Z)*U
 162 *>                  where U is the orthogonal matrix in (*)
 163 *>                  (regardless of the value of WANTT.)
 164 *> \endverbatim
 165 *
 166 *  Authors:
 167 *  ========
 168 *
 169 *> \author Univ. of Tennessee
 170 *> \author Univ. of California Berkeley
 171 *> \author Univ. of Colorado Denver
 172 *> \author NAG Ltd.
 173 *
 174 *> \date November 2015
 175 *
 176 *> \ingroup complexOTHERauxiliary
 177 *
 178 *> \par Contributors:
 179 *  ==================
 180 *>
 181 *> \verbatim
 182 *>
 183 *>     02-96 Based on modifications by
 184 *>     David Day, Sandia National Laboratory, USA
 185 *>
 186 *>     12-04 Further modifications by
 187 *>     Ralph Byers, University of Kansas, USA
 188 *>     This is a modified version of CLAHQR from LAPACK version 3.0.
 189 *>     It is (1) more robust against overflow and underflow and
 190 *>     (2) adopts the more conservative Ahues & Tisseur stopping
 191 *>     criterion (LAWN 122, 1997).
 192 *> \endverbatim
 193 *>
 194 *  =====================================================================
 195       SUBROUTINE CLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
 196      $                   IHIZ, Z, LDZ, INFO )
 197 *
 198 *  -- LAPACK auxiliary routine (version 3.6.0) --
 199 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 200 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 201 *     November 2015
 202 *
 203 *     .. Scalar Arguments ..
 204       INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
 205       LOGICAL            WANTT, WANTZ
 206 *     ..
 207 *     .. Array Arguments ..
 208       COMPLEX            H( LDH, * ), W( * ), Z( LDZ, * )
 209 *     ..
 210 *
 211 *  =========================================================
 212 *
 213 *     .. Parameters ..
 214       COMPLEX            ZERO, ONE
 215       PARAMETER          ( ZERO = ( 0.0e0, 0.0e0 ),
 216      $                   ONE = ( 1.0e0, 0.0e0 ) )
 217       REAL               RZERO, RONE, HALF
 218       PARAMETER          ( RZERO = 0.0e0, RONE = 1.0e0, HALF = 0.5e0 )
 219       REAL               DAT1
 220       PARAMETER          ( DAT1 = 3.0e0 / 4.0e0 )
 221 *     ..
 222 *     .. Local Scalars ..
 223       COMPLEX            CDUM, H11, H11S, H22, SC, SUM, T, T1, TEMP, U,
 224      $                   V2, X, Y
 225       REAL               AA, AB, BA, BB, H10, H21, RTEMP, S, SAFMAX,
 226      $                   SAFMIN, SMLNUM, SX, T2, TST, ULP
 227       INTEGER            I, I1, I2, ITS, ITMAX, J, JHI, JLO, K, L, M,
 228      $                   NH, NZ
 229 *     ..
 230 *     .. Local Arrays ..
 231       COMPLEX            V( 2 )
 232 *     ..
 233 *     .. External Functions ..
 234       COMPLEX            CLADIV
 235       REAL               SLAMCH
 236       EXTERNAL           CLADIV, SLAMCH
 237 *     ..
 238 *     .. External Subroutines ..
 239       EXTERNAL           CCOPY, CLARFG, CSCAL, SLABAD
 240 *     ..
 241 *     .. Statement Functions ..
 242       REAL               CABS1
 243 *     ..
 244 *     .. Intrinsic Functions ..
 245       INTRINSIC          ABS, AIMAG, CONJG, MAX, MIN, REAL, SQRT
 246 *     ..
 247 *     .. Statement Function definitions ..
 248       CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
 249 *     ..
 250 *     .. Executable Statements ..
 251 *
 252       INFO = 0
 253 *
 254 *     Quick return if possible
 255 *
 256       IF( N.EQ.0 )
 257      $   RETURN
 258       IF( ILO.EQ.IHI ) THEN
 259          W( ILO ) = H( ILO, ILO )
 260          RETURN
 261       END IF
 262 *
 263 *     ==== clear out the trash ====
 264       DO 10 J = ILO, IHI - 3
 265          H( J+2, J ) = ZERO
 266          H( J+3, J ) = ZERO
 267    10 CONTINUE
 268       IF( ILO.LE.IHI-2 )
 269      $   H( IHI, IHI-2 ) = ZERO
 270 *     ==== ensure that subdiagonal entries are real ====
 271       IF( WANTT ) THEN
 272          JLO = 1
 273          JHI = N
 274       ELSE
 275          JLO = ILO
 276          JHI = IHI
 277       END IF
 278       DO 20 I = ILO + 1, IHI
 279          IF( AIMAG( H( I, I-1 ) ).NE.RZERO ) THEN
 280 *           ==== The following redundant normalization
 281 *           .    avoids problems with both gradual and
 282 *           .    sudden underflow in ABS(H(I,I-1)) ====
 283             SC = H( I, I-1 ) / CABS1( H( I, I-1 ) )
 284             SC = CONJG( SC ) / ABS( SC )
 285             H( I, I-1 ) = ABS( H( I, I-1 ) )
 286             CALL CSCAL( JHI-I+1, SC, H( I, I ), LDH )
 287             CALL CSCAL( MIN( JHI, I+1 )-JLO+1, CONJG( SC ), H( JLO, I ),
 288      $                  1 )
 289             IF( WANTZ )
 290      $         CALL CSCAL( IHIZ-ILOZ+1, CONJG( SC ), Z( ILOZ, I ), 1 )
 291          END IF
 292    20 CONTINUE
 293 *
 294       NH = IHI - ILO + 1
 295       NZ = IHIZ - ILOZ + 1
 296 *
 297 *     Set machine-dependent constants for the stopping criterion.
 298 *
 299       SAFMIN = SLAMCH( 'SAFE MINIMUM' )
 300       SAFMAX = RONE / SAFMIN
 301       CALL SLABAD( SAFMIN, SAFMAX )
 302       ULP = SLAMCH( 'PRECISION' )
 303       SMLNUM = SAFMIN*( REAL( NH ) / ULP )
 304 *
 305 *     I1 and I2 are the indices of the first row and last column of H
 306 *     to which transformations must be applied. If eigenvalues only are
 307 *     being computed, I1 and I2 are set inside the main loop.
 308 *
 309       IF( WANTT ) THEN
 310          I1 = 1
 311          I2 = N
 312       END IF
 313 *
 314 *     ITMAX is the total number of QR iterations allowed.
 315 *
 316       ITMAX = 30 * MAX( 10, NH )
 317 *
 318 *     The main loop begins here. I is the loop index and decreases from
 319 *     IHI to ILO in steps of 1. Each iteration of the loop works
 320 *     with the active submatrix in rows and columns L to I.
 321 *     Eigenvalues I+1 to IHI have already converged. Either L = ILO, or
 322 *     H(L,L-1) is negligible so that the matrix splits.
 323 *
 324       I = IHI
 325    30 CONTINUE
 326       IF( I.LT.ILO )
 327      $   GO TO 150
 328 *
 329 *     Perform QR iterations on rows and columns ILO to I until a
 330 *     submatrix of order 1 splits off at the bottom because a
 331 *     subdiagonal element has become negligible.
 332 *
 333       L = ILO
 334       DO 130 ITS = 0, ITMAX
 335 *
 336 *        Look for a single small subdiagonal element.
 337 *
 338          DO 40 K = I, L + 1, -1
 339             IF( CABS1( H( K, K-1 ) ).LE.SMLNUM )
 340      $         GO TO 50
 341             TST = CABS1( H( K-1, K-1 ) ) + CABS1( H( K, K ) )
 342             IF( TST.EQ.ZERO ) THEN
 343                IF( K-2.GE.ILO )
 344      $            TST = TST + ABS( REAL( H( K-1, K-2 ) ) )
 345                IF( K+1.LE.IHI )
 346      $            TST = TST + ABS( REAL( H( K+1, K ) ) )
 347             END IF
 348 *           ==== The following is a conservative small subdiagonal
 349 *           .    deflation criterion due to Ahues & Tisseur (LAWN 122,
 350 *           .    1997). It has better mathematical foundation and
 351 *           .    improves accuracy in some examples.  ====
 352             IF( ABS( REAL( H( K, K-1 ) ) ).LE.ULP*TST ) THEN
 353                AB = MAX( CABS1( H( K, K-1 ) ), CABS1( H( K-1, K ) ) )
 354                BA = MIN( CABS1( H( K, K-1 ) ), CABS1( H( K-1, K ) ) )
 355                AA = MAX( CABS1( H( K, K ) ),
 356      $              CABS1( H( K-1, K-1 )-H( K, K ) ) )
 357                BB = MIN( CABS1( H( K, K ) ),
 358      $              CABS1( H( K-1, K-1 )-H( K, K ) ) )
 359                S = AA + AB
 360                IF( BA*( AB / S ).LE.MAX( SMLNUM,
 361      $             ULP*( BB*( AA / S ) ) ) )GO TO 50
 362             END IF
 363    40    CONTINUE
 364    50    CONTINUE
 365          L = K
 366          IF( L.GT.ILO ) THEN
 367 *
 368 *           H(L,L-1) is negligible
 369 *
 370             H( L, L-1 ) = ZERO
 371          END IF
 372 *
 373 *        Exit from loop if a submatrix of order 1 has split off.
 374 *
 375          IF( L.GE.I )
 376      $      GO TO 140
 377 *
 378 *        Now the active submatrix is in rows and columns L to I. If
 379 *        eigenvalues only are being computed, only the active submatrix
 380 *        need be transformed.
 381 *
 382          IF( .NOT.WANTT ) THEN
 383             I1 = L
 384             I2 = I
 385          END IF
 386 *
 387          IF( ITS.EQ.10 ) THEN
 388 *
 389 *           Exceptional shift.
 390 *
 391             S = DAT1*ABS( REAL( H( L+1, L ) ) )
 392             T = S + H( L, L )
 393          ELSE IF( ITS.EQ.20 ) THEN
 394 *
 395 *           Exceptional shift.
 396 *
 397             S = DAT1*ABS( REAL( H( I, I-1 ) ) )
 398             T = S + H( I, I )
 399          ELSE
 400 *
 401 *           Wilkinson's shift.
 402 *
 403             T = H( I, I )
 404             U = SQRT( H( I-1, I ) )*SQRT( H( I, I-1 ) )
 405             S = CABS1( U )
 406             IF( S.NE.RZERO ) THEN
 407                X = HALF*( H( I-1, I-1 )-T )
 408                SX = CABS1( X )
 409                S = MAX( S, CABS1( X ) )
 410                Y = S*SQRT( ( X / S )**2+( U / S )**2 )
 411                IF( SX.GT.RZERO ) THEN
 412                   IF( REAL( X / SX )*REAL( Y )+AIMAG( X / SX )*
 413      $                AIMAG( Y ).LT.RZERO )Y = -Y
 414                END IF
 415                T = T - U*CLADIV( U, ( X+Y ) )
 416             END IF
 417          END IF
 418 *
 419 *        Look for two consecutive small subdiagonal elements.
 420 *
 421          DO 60 M = I - 1, L + 1, -1
 422 *
 423 *           Determine the effect of starting the single-shift QR
 424 *           iteration at row M, and see if this would make H(M,M-1)
 425 *           negligible.
 426 *
 427             H11 = H( M, M )
 428             H22 = H( M+1, M+1 )
 429             H11S = H11 - T
 430             H21 = REAL( H( M+1, M ) )
 431             S = CABS1( H11S ) + ABS( H21 )
 432             H11S = H11S / S
 433             H21 = H21 / S
 434             V( 1 ) = H11S
 435             V( 2 ) = H21
 436             H10 = REAL( H( M, M-1 ) )
 437             IF( ABS( H10 )*ABS( H21 ).LE.ULP*
 438      $          ( CABS1( H11S )*( CABS1( H11 )+CABS1( H22 ) ) ) )
 439      $          GO TO 70
 440    60    CONTINUE
 441          H11 = H( L, L )
 442          H22 = H( L+1, L+1 )
 443          H11S = H11 - T
 444          H21 = REAL( H( L+1, L ) )
 445          S = CABS1( H11S ) + ABS( H21 )
 446          H11S = H11S / S
 447          H21 = H21 / S
 448          V( 1 ) = H11S
 449          V( 2 ) = H21
 450    70    CONTINUE
 451 *
 452 *        Single-shift QR step
 453 *
 454          DO 120 K = M, I - 1
 455 *
 456 *           The first iteration of this loop determines a reflection G
 457 *           from the vector V and applies it from left and right to H,
 458 *           thus creating a nonzero bulge below the subdiagonal.
 459 *
 460 *           Each subsequent iteration determines a reflection G to
 461 *           restore the Hessenberg form in the (K-1)th column, and thus
 462 *           chases the bulge one step toward the bottom of the active
 463 *           submatrix.
 464 *
 465 *           V(2) is always real before the call to CLARFG, and hence
 466 *           after the call T2 ( = T1*V(2) ) is also real.
 467 *
 468             IF( K.GT.M )
 469      $         CALL CCOPY( 2, H( K, K-1 ), 1, V, 1 )
 470             CALL CLARFG( 2, V( 1 ), V( 2 ), 1, T1 )
 471             IF( K.GT.M ) THEN
 472                H( K, K-1 ) = V( 1 )
 473                H( K+1, K-1 ) = ZERO
 474             END IF
 475             V2 = V( 2 )
 476             T2 = REAL( T1*V2 )
 477 *
 478 *           Apply G from the left to transform the rows of the matrix
 479 *           in columns K to I2.
 480 *
 481             DO 80 J = K, I2
 482                SUM = CONJG( T1 )*H( K, J ) + T2*H( K+1, J )
 483                H( K, J ) = H( K, J ) - SUM
 484                H( K+1, J ) = H( K+1, J ) - SUM*V2
 485    80       CONTINUE
 486 *
 487 *           Apply G from the right to transform the columns of the
 488 *           matrix in rows I1 to min(K+2,I).
 489 *
 490             DO 90 J = I1, MIN( K+2, I )
 491                SUM = T1*H( J, K ) + T2*H( J, K+1 )
 492                H( J, K ) = H( J, K ) - SUM
 493                H( J, K+1 ) = H( J, K+1 ) - SUM*CONJG( V2 )
 494    90       CONTINUE
 495 *
 496             IF( WANTZ ) THEN
 497 *
 498 *              Accumulate transformations in the matrix Z
 499 *
 500                DO 100 J = ILOZ, IHIZ
 501                   SUM = T1*Z( J, K ) + T2*Z( J, K+1 )
 502                   Z( J, K ) = Z( J, K ) - SUM
 503                   Z( J, K+1 ) = Z( J, K+1 ) - SUM*CONJG( V2 )
 504   100          CONTINUE
 505             END IF
 506 *
 507             IF( K.EQ.M .AND. M.GT.L ) THEN
 508 *
 509 *              If the QR step was started at row M > L because two
 510 *              consecutive small subdiagonals were found, then extra
 511 *              scaling must be performed to ensure that H(M,M-1) remains
 512 *              real.
 513 *
 514                TEMP = ONE - T1
 515                TEMP = TEMP / ABS( TEMP )
 516                H( M+1, M ) = H( M+1, M )*CONJG( TEMP )
 517                IF( M+2.LE.I )
 518      $            H( M+2, M+1 ) = H( M+2, M+1 )*TEMP
 519                DO 110 J = M, I
 520                   IF( J.NE.M+1 ) THEN
 521                      IF( I2.GT.J )
 522      $                  CALL CSCAL( I2-J, TEMP, H( J, J+1 ), LDH )
 523                      CALL CSCAL( J-I1, CONJG( TEMP ), H( I1, J ), 1 )
 524                      IF( WANTZ ) THEN
 525                         CALL CSCAL( NZ, CONJG( TEMP ), Z( ILOZ, J ), 1 )
 526                      END IF
 527                   END IF
 528   110          CONTINUE
 529             END IF
 530   120    CONTINUE
 531 *
 532 *        Ensure that H(I,I-1) is real.
 533 *
 534          TEMP = H( I, I-1 )
 535          IF( AIMAG( TEMP ).NE.RZERO ) THEN
 536             RTEMP = ABS( TEMP )
 537             H( I, I-1 ) = RTEMP
 538             TEMP = TEMP / RTEMP
 539             IF( I2.GT.I )
 540      $         CALL CSCAL( I2-I, CONJG( TEMP ), H( I, I+1 ), LDH )
 541             CALL CSCAL( I-I1, TEMP, H( I1, I ), 1 )
 542             IF( WANTZ ) THEN
 543                CALL CSCAL( NZ, TEMP, Z( ILOZ, I ), 1 )
 544             END IF
 545          END IF
 546 *
 547   130 CONTINUE
 548 *
 549 *     Failure to converge in remaining number of iterations
 550 *
 551       INFO = I
 552       RETURN
 553 *
 554   140 CONTINUE
 555 *
 556 *     H(I,I-1) is negligible: one eigenvalue has converged.
 557 *
 558       W( I ) = H( I, I )
 559 *
 560 *     return to start of the main loop with new value of I.
 561 *
 562       I = L - 1
 563       GO TO 30
 564 *
 565   150 CONTINUE
 566       RETURN
 567 *
 568 *     End of CLAHQR
 569 *
 570       END