3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download ZTGSEN + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztgsen.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztgsen.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztgsen.f">
21 * SUBROUTINE ZTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
22 * ALPHA, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF,
23 * WORK, LWORK, IWORK, LIWORK, INFO )
25 * .. Scalar Arguments ..
26 * LOGICAL WANTQ, WANTZ
27 * INTEGER IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK,
29 * DOUBLE PRECISION PL, PR
31 * .. Array Arguments ..
34 * DOUBLE PRECISION DIF( * )
35 * COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
36 * $ BETA( * ), Q( LDQ, * ), WORK( * ), Z( LDZ, * )
45 *> ZTGSEN reorders the generalized Schur decomposition of a complex
46 *> matrix pair (A, B) (in terms of an unitary equivalence trans-
47 *> formation Q**H * (A, B) * Z), so that a selected cluster of eigenvalues
48 *> appears in the leading diagonal blocks of the pair (A,B). The leading
49 *> columns of Q and Z form unitary bases of the corresponding left and
50 *> right eigenspaces (deflating subspaces). (A, B) must be in
51 *> generalized Schur canonical form, that is, A and B are both upper
54 *> ZTGSEN also computes the generalized eigenvalues
56 *> w(j)= ALPHA(j) / BETA(j)
58 *> of the reordered matrix pair (A, B).
60 *> Optionally, the routine computes estimates of reciprocal condition
61 *> numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11),
62 *> (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s)
63 *> between the matrix pairs (A11, B11) and (A22,B22) that correspond to
64 *> the selected cluster and the eigenvalues outside the cluster, resp.,
65 *> and norms of "projections" onto left and right eigenspaces w.r.t.
66 *> the selected cluster in the (1,1)-block.
76 *> Specifies whether condition numbers are required for the
77 *> cluster of eigenvalues (PL and PR) or the deflating subspaces
79 *> =0: Only reorder w.r.t. SELECT. No extras.
80 *> =1: Reciprocal of norms of "projections" onto left and right
81 *> eigenspaces w.r.t. the selected cluster (PL and PR).
82 *> =2: Upper bounds on Difu and Difl. F-norm-based estimate
84 *> =3: Estimate of Difu and Difl. 1-norm-based estimate
86 *> About 5 times as expensive as IJOB = 2.
87 *> =4: Compute PL, PR and DIF (i.e. 0, 1 and 2 above): Economic
88 *> version to get it all.
89 *> =5: Compute PL, PR and DIF (i.e. 0, 1 and 3 above)
95 *> .TRUE. : update the left transformation matrix Q;
96 *> .FALSE.: do not update Q.
102 *> .TRUE. : update the right transformation matrix Z;
103 *> .FALSE.: do not update Z.
108 *> SELECT is LOGICAL array, dimension (N)
109 *> SELECT specifies the eigenvalues in the selected cluster. To
110 *> select an eigenvalue w(j), SELECT(j) must be set to
117 *> The order of the matrices A and B. N >= 0.
122 *> A is COMPLEX*16 array, dimension(LDA,N)
123 *> On entry, the upper triangular matrix A, in generalized
124 *> Schur canonical form.
125 *> On exit, A is overwritten by the reordered matrix A.
131 *> The leading dimension of the array A. LDA >= max(1,N).
136 *> B is COMPLEX*16 array, dimension(LDB,N)
137 *> On entry, the upper triangular matrix B, in generalized
138 *> Schur canonical form.
139 *> On exit, B is overwritten by the reordered matrix B.
145 *> The leading dimension of the array B. LDB >= max(1,N).
150 *> ALPHA is COMPLEX*16 array, dimension (N)
155 *> BETA is COMPLEX*16 array, dimension (N)
157 *> The diagonal elements of A and B, respectively,
158 *> when the pair (A,B) has been reduced to generalized Schur
159 *> form. ALPHA(i)/BETA(i) i=1,...,N are the generalized
165 *> Q is COMPLEX*16 array, dimension (LDQ,N)
166 *> On entry, if WANTQ = .TRUE., Q is an N-by-N matrix.
167 *> On exit, Q has been postmultiplied by the left unitary
168 *> transformation matrix which reorder (A, B); The leading M
169 *> columns of Q form orthonormal bases for the specified pair of
170 *> left eigenspaces (deflating subspaces).
171 *> If WANTQ = .FALSE., Q is not referenced.
177 *> The leading dimension of the array Q. LDQ >= 1.
178 *> If WANTQ = .TRUE., LDQ >= N.
183 *> Z is COMPLEX*16 array, dimension (LDZ,N)
184 *> On entry, if WANTZ = .TRUE., Z is an N-by-N matrix.
185 *> On exit, Z has been postmultiplied by the left unitary
186 *> transformation matrix which reorder (A, B); The leading M
187 *> columns of Z form orthonormal bases for the specified pair of
188 *> left eigenspaces (deflating subspaces).
189 *> If WANTZ = .FALSE., Z is not referenced.
195 *> The leading dimension of the array Z. LDZ >= 1.
196 *> If WANTZ = .TRUE., LDZ >= N.
202 *> The dimension of the specified pair of left and right
203 *> eigenspaces, (deflating subspaces) 0 <= M <= N.
208 *> PL is DOUBLE PRECISION
213 *> PR is DOUBLE PRECISION
215 *> If IJOB = 1, 4 or 5, PL, PR are lower bounds on the
216 *> reciprocal of the norm of "projections" onto left and right
217 *> eigenspace with respect to the selected cluster.
219 *> If M = 0 or M = N, PL = PR = 1.
220 *> If IJOB = 0, 2 or 3 PL, PR are not referenced.
225 *> DIF is DOUBLE PRECISION array, dimension (2).
226 *> If IJOB >= 2, DIF(1:2) store the estimates of Difu and Difl.
227 *> If IJOB = 2 or 4, DIF(1:2) are F-norm-based upper bounds on
228 *> Difu and Difl. If IJOB = 3 or 5, DIF(1:2) are 1-norm-based
229 *> estimates of Difu and Difl, computed using reversed
230 *> communication with ZLACN2.
231 *> If M = 0 or N, DIF(1:2) = F-norm([A, B]).
232 *> If IJOB = 0 or 1, DIF is not referenced.
237 *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
238 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
244 *> The dimension of the array WORK. LWORK >= 1
245 *> If IJOB = 1, 2 or 4, LWORK >= 2*M*(N-M)
246 *> If IJOB = 3 or 5, LWORK >= 4*M*(N-M)
248 *> If LWORK = -1, then a workspace query is assumed; the routine
249 *> only calculates the optimal size of the WORK array, returns
250 *> this value as the first entry of the WORK array, and no error
251 *> message related to LWORK is issued by XERBLA.
256 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
257 *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
263 *> The dimension of the array IWORK. LIWORK >= 1.
264 *> If IJOB = 1, 2 or 4, LIWORK >= N+2;
265 *> If IJOB = 3 or 5, LIWORK >= MAX(N+2, 2*M*(N-M));
267 *> If LIWORK = -1, then a workspace query is assumed; the
268 *> routine only calculates the optimal size of the IWORK array,
269 *> returns this value as the first entry of the IWORK array, and
270 *> no error message related to LIWORK is issued by XERBLA.
276 *> =0: Successful exit.
277 *> <0: If INFO = -i, the i-th argument had an illegal value.
278 *> =1: Reordering of (A, B) failed because the transformed
279 *> matrix pair (A, B) would be too far from generalized
280 *> Schur form; the problem is very ill-conditioned.
281 *> (A, B) may have been partially reordered.
282 *> If requested, 0 is returned in DIF(*), PL and PR.
288 *> \author Univ. of Tennessee
289 *> \author Univ. of California Berkeley
290 *> \author Univ. of Colorado Denver
295 *> \ingroup complex16OTHERcomputational
297 *> \par Further Details:
298 * =====================
302 *> ZTGSEN first collects the selected eigenvalues by computing unitary
303 *> U and W that move them to the top left corner of (A, B). In other
304 *> words, the selected eigenvalues are the eigenvalues of (A11, B11) in
306 *> U**H*(A, B)*W = (A11 A12) (B11 B12) n1
307 *> ( 0 A22),( 0 B22) n2
310 *> where N = n1+n2 and U**H means the conjugate transpose of U. The first
311 *> n1 columns of U and W span the specified pair of left and right
312 *> eigenspaces (deflating subspaces) of (A, B).
314 *> If (A, B) has been obtained from the generalized real Schur
315 *> decomposition of a matrix pair (C, D) = Q*(A, B)*Z**H, then the
316 *> reordered generalized Schur form of (C, D) is given by
318 *> (C, D) = (Q*U)*(U**H *(A, B)*W)*(Z*W)**H,
320 *> and the first n1 columns of Q*U and Z*W span the corresponding
321 *> deflating subspaces of (C, D) (Q and Z store Q*U and Z*W, resp.).
323 *> Note that if the selected eigenvalue is sufficiently ill-conditioned,
324 *> then its value may differ significantly from its value before
327 *> The reciprocal condition numbers of the left and right eigenspaces
328 *> spanned by the first n1 columns of U and W (or Q*U and Z*W) may
329 *> be returned in DIF(1:2), corresponding to Difu and Difl, resp.
331 *> The Difu and Difl are defined as:
333 *> Difu[(A11, B11), (A22, B22)] = sigma-min( Zu )
335 *> Difl[(A11, B11), (A22, B22)] = Difu[(A22, B22), (A11, B11)],
337 *> where sigma-min(Zu) is the smallest singular value of the
338 *> (2*n1*n2)-by-(2*n1*n2) matrix
340 *> Zu = [ kron(In2, A11) -kron(A22**H, In1) ]
341 *> [ kron(In2, B11) -kron(B22**H, In1) ].
343 *> Here, Inx is the identity matrix of size nx and A22**H is the
344 *> conjugate transpose of A22. kron(X, Y) is the Kronecker product between
345 *> the matrices X and Y.
347 *> When DIF(2) is small, small changes in (A, B) can cause large changes
348 *> in the deflating subspace. An approximate (asymptotic) bound on the
349 *> maximum angular error in the computed deflating subspaces is
351 *> EPS * norm((A, B)) / DIF(2),
353 *> where EPS is the machine precision.
355 *> The reciprocal norm of the projectors on the left and right
356 *> eigenspaces associated with (A11, B11) may be returned in PL and PR.
357 *> They are computed as follows. First we compute L and R so that
358 *> P*(A, B)*Q is block diagonal, where
360 *> P = ( I -L ) n1 Q = ( I R ) n1
361 *> ( 0 I ) n2 and ( 0 I ) n2
364 *> and (L, R) is the solution to the generalized Sylvester equation
366 *> A11*R - L*A22 = -A12
367 *> B11*R - L*B22 = -B12
369 *> Then PL = (F-norm(L)**2+1)**(-1/2) and PR = (F-norm(R)**2+1)**(-1/2).
370 *> An approximate (asymptotic) bound on the average absolute error of
371 *> the selected eigenvalues is
373 *> EPS * norm((A, B)) / PL.
375 *> There are also global error bounds which valid for perturbations up
376 *> to a certain restriction: A lower bound (x) on the smallest
377 *> F-norm(E,F) for which an eigenvalue of (A11, B11) may move and
378 *> coalesce with an eigenvalue of (A22, B22) under perturbation (E,F),
379 *> (i.e. (A + E, B + F), is
381 *> x = min(Difu,Difl)/((1/(PL*PL)+1/(PR*PR))**(1/2)+2*max(1/PL,1/PR)).
383 *> An approximate bound on x can be computed from DIF(1:2), PL and PR.
385 *> If y = ( F-norm(E,F) / x) <= 1, the angles between the perturbed
386 *> (L', R') and unperturbed (L, R) left and right deflating subspaces
387 *> associated with the selected cluster in the (1,1)-blocks can be
390 *> max-angle(L, L') <= arctan( y * PL / (1 - y * (1 - PL * PL)**(1/2))
391 *> max-angle(R, R') <= arctan( y * PR / (1 - y * (1 - PR * PR)**(1/2))
393 *> See LAPACK User's Guide section 4.11 or the following references
394 *> for more information.
396 *> Note that if the default method for computing the Frobenius-norm-
397 *> based estimate DIF is not wanted (see ZLATDF), then the parameter
398 *> IDIFJB (see below) should be changed from 3 to 4 (routine ZLATDF
399 *> (IJOB = 2 will be used)). See ZTGSYL for more details.
402 *> \par Contributors:
405 *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
406 *> Umea University, S-901 87 Umea, Sweden.
411 *> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
412 *> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
413 *> M.S. Moonen et al (eds), Linear Algebra for Large Scale and
414 *> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
416 *> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
417 *> Eigenvalues of a Regular Matrix Pair (A, B) and Condition
418 *> Estimation: Theory, Algorithms and Software, Report
419 *> UMINF - 94.04, Department of Computing Science, Umea University,
420 *> S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87.
421 *> To appear in Numerical Algorithms, 1996.
423 *> [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
424 *> for Solving the Generalized Sylvester Equation and Estimating the
425 *> Separation between Regular Matrix Pairs, Report UMINF - 93.23,
426 *> Department of Computing Science, Umea University, S-901 87 Umea,
427 *> Sweden, December 1993, Revised April 1994, Also as LAPACK working
428 *> Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1,
431 * =====================================================================
432 SUBROUTINE ZTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
433 $ ALPHA, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF,
434 $ WORK, LWORK, IWORK, LIWORK, INFO )
436 * -- LAPACK computational routine (version 3.6.1) --
437 * -- LAPACK is a software package provided by Univ. of Tennessee, --
438 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
441 * .. Scalar Arguments ..
443 INTEGER IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK,
445 DOUBLE PRECISION PL, PR
447 * .. Array Arguments ..
450 DOUBLE PRECISION DIF( * )
451 COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
452 $ BETA( * ), Q( LDQ, * ), WORK( * ), Z( LDZ, * )
455 * =====================================================================
459 PARAMETER ( IDIFJB = 3 )
460 DOUBLE PRECISION ZERO, ONE
461 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
463 * .. Local Scalars ..
464 LOGICAL LQUERY, SWAP, WANTD, WANTD1, WANTD2, WANTP
465 INTEGER I, IERR, IJB, K, KASE, KS, LIWMIN, LWMIN, MN2,
467 DOUBLE PRECISION DSCALE, DSUM, RDSCAL, SAFMIN
468 COMPLEX*16 TEMP1, TEMP2
473 * .. External Subroutines ..
474 EXTERNAL XERBLA, ZLACN2, ZLACPY, ZLASSQ, ZSCAL, ZTGEXC,
477 * .. Intrinsic Functions ..
478 INTRINSIC ABS, DCMPLX, DCONJG, MAX, SQRT
480 * .. External Functions ..
481 DOUBLE PRECISION DLAMCH
484 * .. Executable Statements ..
486 * Decode and test the input parameters
489 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
491 IF( IJOB.LT.0 .OR. IJOB.GT.5 ) THEN
493 ELSE IF( N.LT.0 ) THEN
495 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
497 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
499 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
501 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
506 CALL XERBLA( 'ZTGSEN', -INFO )
512 WANTP = IJOB.EQ.1 .OR. IJOB.GE.4
513 WANTD1 = IJOB.EQ.2 .OR. IJOB.EQ.4
514 WANTD2 = IJOB.EQ.3 .OR. IJOB.EQ.5
515 WANTD = WANTD1 .OR. WANTD2
517 * Set M to the dimension of the specified pair of deflating
521 IF( .NOT.LQUERY .OR. IJOB.NE.0 ) THEN
523 ALPHA( K ) = A( K, K )
524 BETA( K ) = B( K, K )
535 IF( IJOB.EQ.1 .OR. IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
536 LWMIN = MAX( 1, 2*M*( N-M ) )
537 LIWMIN = MAX( 1, N+2 )
538 ELSE IF( IJOB.EQ.3 .OR. IJOB.EQ.5 ) THEN
539 LWMIN = MAX( 1, 4*M*( N-M ) )
540 LIWMIN = MAX( 1, 2*M*( N-M ), N+2 )
549 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
551 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
556 CALL XERBLA( 'ZTGSEN', -INFO )
558 ELSE IF( LQUERY ) THEN
562 * Quick return if possible.
564 IF( M.EQ.N .OR. M.EQ.0 ) THEN
573 CALL ZLASSQ( N, A( 1, I ), 1, DSCALE, DSUM )
574 CALL ZLASSQ( N, B( 1, I ), 1, DSCALE, DSUM )
576 DIF( 1 ) = DSCALE*SQRT( DSUM )
582 * Get machine constant
584 SAFMIN = DLAMCH( 'S' )
586 * Collect the selected blocks at the top-left corner of (A, B).
594 * Swap the K-th block to position KS. Compute unitary Q
595 * and Z that will swap adjacent diagonal blocks in (A, B).
598 $ CALL ZTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
603 * Swap is rejected: exit.
620 * Solve generalized Sylvester equation for R and L:
621 * A11 * R - L * A22 = A12
622 * B11 * R - L * B22 = B12
627 CALL ZLACPY( 'Full', N1, N2, A( 1, I ), LDA, WORK, N1 )
628 CALL ZLACPY( 'Full', N1, N2, B( 1, I ), LDB, WORK( N1*N2+1 ),
631 CALL ZTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK,
632 $ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ), N1,
633 $ DSCALE, DIF( 1 ), WORK( N1*N2*2+1 ),
634 $ LWORK-2*N1*N2, IWORK, IERR )
636 * Estimate the reciprocal of norms of "projections" onto
637 * left and right eigenspaces
641 CALL ZLASSQ( N1*N2, WORK, 1, RDSCAL, DSUM )
642 PL = RDSCAL*SQRT( DSUM )
643 IF( PL.EQ.ZERO ) THEN
646 PL = DSCALE / ( SQRT( DSCALE*DSCALE / PL+PL )*SQRT( PL ) )
650 CALL ZLASSQ( N1*N2, WORK( N1*N2+1 ), 1, RDSCAL, DSUM )
651 PR = RDSCAL*SQRT( DSUM )
652 IF( PR.EQ.ZERO ) THEN
655 PR = DSCALE / ( SQRT( DSCALE*DSCALE / PR+PR )*SQRT( PR ) )
660 * Compute estimates Difu and Difl.
668 * Frobenius norm-based Difu estimate.
670 CALL ZTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK,
671 $ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ),
672 $ N1, DSCALE, DIF( 1 ), WORK( N1*N2*2+1 ),
673 $ LWORK-2*N1*N2, IWORK, IERR )
675 * Frobenius norm-based Difl estimate.
677 CALL ZTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA, WORK,
678 $ N2, B( I, I ), LDB, B, LDB, WORK( N1*N2+1 ),
679 $ N2, DSCALE, DIF( 2 ), WORK( N1*N2*2+1 ),
680 $ LWORK-2*N1*N2, IWORK, IERR )
683 * Compute 1-norm-based estimates of Difu and Difl using
684 * reversed communication with ZLACN2. In each step a
685 * generalized Sylvester equation or a transposed variant
695 * 1-norm-based estimate of Difu.
698 CALL ZLACN2( MN2, WORK( MN2+1 ), WORK, DIF( 1 ), KASE,
703 * Solve generalized Sylvester equation
705 CALL ZTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA,
706 $ WORK, N1, B, LDB, B( I, I ), LDB,
707 $ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ),
708 $ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,
712 * Solve the transposed variant.
714 CALL ZTGSYL( 'C', IJB, N1, N2, A, LDA, A( I, I ), LDA,
715 $ WORK, N1, B, LDB, B( I, I ), LDB,
716 $ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ),
717 $ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,
722 DIF( 1 ) = DSCALE / DIF( 1 )
724 * 1-norm-based estimate of Difl.
727 CALL ZLACN2( MN2, WORK( MN2+1 ), WORK, DIF( 2 ), KASE,
732 * Solve generalized Sylvester equation
734 CALL ZTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA,
735 $ WORK, N2, B( I, I ), LDB, B, LDB,
736 $ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ),
737 $ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,
741 * Solve the transposed variant.
743 CALL ZTGSYL( 'C', IJB, N2, N1, A( I, I ), LDA, A, LDA,
744 $ WORK, N2, B, LDB, B( I, I ), LDB,
745 $ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ),
746 $ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,
751 DIF( 2 ) = DSCALE / DIF( 2 )
755 * If B(K,K) is complex, make it real and positive (normalization
756 * of the generalized Schur form) and Store the generalized
757 * eigenvalues of reordered pair (A, B)
760 DSCALE = ABS( B( K, K ) )
761 IF( DSCALE.GT.SAFMIN ) THEN
762 TEMP1 = DCONJG( B( K, K ) / DSCALE )
763 TEMP2 = B( K, K ) / DSCALE
765 CALL ZSCAL( N-K, TEMP1, B( K, K+1 ), LDB )
766 CALL ZSCAL( N-K+1, TEMP1, A( K, K ), LDA )
768 $ CALL ZSCAL( N, TEMP2, Q( 1, K ), 1 )
770 B( K, K ) = DCMPLX( ZERO, ZERO )
773 ALPHA( K ) = A( K, K )
774 BETA( K ) = B( K, K )