LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlaqz3.f
Go to the documentation of this file.
1*> \brief \b DLAQZ3
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download DLAQZ3 + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqz3.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqz3.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqz3.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE DLAQZ3( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NW, A, LDA, B,
20* $ LDB, Q, LDQ, Z, LDZ, NS, ND, ALPHAR, ALPHAI, BETA, QC, LDQC,
21* $ ZC, LDZC, WORK, LWORK, REC, INFO )
22* IMPLICIT NONE
23*
24* Arguments
25* LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
26* INTEGER, INTENT( IN ) :: N, ILO, IHI, NW, LDA, LDB, LDQ, LDZ,
27* $ LDQC, LDZC, LWORK, REC
28*
29* DOUBLE PRECISION, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ),
30* $ Q( LDQ, * ), Z( LDZ, * ), ALPHAR( * ), ALPHAI( * ), BETA( * )
31* INTEGER, INTENT( OUT ) :: NS, ND, INFO
32* DOUBLE PRECISION :: QC( LDQC, * ), ZC( LDZC, * ), WORK( * )
33* ..
34*
35*
36*> \par Purpose:
37* =============
38*>
39*> \verbatim
40*>
41*> DLAQZ3 performs AED
42*> \endverbatim
43*
44* Arguments:
45* ==========
46*
47*> \param[in] ILSCHUR
48*> \verbatim
49*> ILSCHUR is LOGICAL
50*> Determines whether or not to update the full Schur form
51*> \endverbatim
52*>
53*> \param[in] ILQ
54*> \verbatim
55*> ILQ is LOGICAL
56*> Determines whether or not to update the matrix Q
57*> \endverbatim
58*>
59*> \param[in] ILZ
60*> \verbatim
61*> ILZ is LOGICAL
62*> Determines whether or not to update the matrix Z
63*> \endverbatim
64*>
65*> \param[in] N
66*> \verbatim
67*> N is INTEGER
68*> The order of the matrices A, B, Q, and Z. N >= 0.
69*> \endverbatim
70*>
71*> \param[in] ILO
72*> \verbatim
73*> ILO is INTEGER
74*> \endverbatim
75*>
76*> \param[in] IHI
77*> \verbatim
78*> IHI is INTEGER
79*> ILO and IHI mark the rows and columns of (A,B) which
80*> are to be normalized
81*> \endverbatim
82*>
83*> \param[in] NW
84*> \verbatim
85*> NW is INTEGER
86*> The desired size of the deflation window.
87*> \endverbatim
88*>
89*> \param[in,out] A
90*> \verbatim
91*> A is DOUBLE PRECISION array, dimension (LDA, N)
92*> \endverbatim
93*>
94*> \param[in] LDA
95*> \verbatim
96*> LDA is INTEGER
97*> The leading dimension of the array A. LDA >= max( 1, N ).
98*> \endverbatim
99*>
100*> \param[in,out] B
101*> \verbatim
102*> B is DOUBLE PRECISION array, dimension (LDB, N)
103*> \endverbatim
104*>
105*> \param[in] LDB
106*> \verbatim
107*> LDB is INTEGER
108*> The leading dimension of the array B. LDB >= max( 1, N ).
109*> \endverbatim
110*>
111*> \param[in,out] Q
112*> \verbatim
113*> Q is DOUBLE PRECISION array, dimension (LDQ, N)
114*> \endverbatim
115*>
116*> \param[in] LDQ
117*> \verbatim
118*> LDQ is INTEGER
119*> \endverbatim
120*>
121*> \param[in,out] Z
122*> \verbatim
123*> Z is DOUBLE PRECISION array, dimension (LDZ, N)
124*> \endverbatim
125*>
126*> \param[in] LDZ
127*> \verbatim
128*> LDZ is INTEGER
129*> \endverbatim
130*>
131*> \param[out] NS
132*> \verbatim
133*> NS is INTEGER
134*> The number of unconverged eigenvalues available to
135*> use as shifts.
136*> \endverbatim
137*>
138*> \param[out] ND
139*> \verbatim
140*> ND is INTEGER
141*> The number of converged eigenvalues found.
142*> \endverbatim
143*>
144*> \param[out] ALPHAR
145*> \verbatim
146*> ALPHAR is DOUBLE PRECISION array, dimension (N)
147*> The real parts of each scalar alpha defining an eigenvalue
148*> of GNEP.
149*> \endverbatim
150*>
151*> \param[out] ALPHAI
152*> \verbatim
153*> ALPHAI is DOUBLE PRECISION array, dimension (N)
154*> The imaginary parts of each scalar alpha defining an
155*> eigenvalue of GNEP.
156*> If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
157*> positive, then the j-th and (j+1)-st eigenvalues are a
158*> complex conjugate pair, with ALPHAI(j+1) = -ALPHAI(j).
159*> \endverbatim
160*>
161*> \param[out] BETA
162*> \verbatim
163*> BETA is DOUBLE PRECISION array, dimension (N)
164*> The scalars beta that define the eigenvalues of GNEP.
165*> Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
166*> beta = BETA(j) represent the j-th eigenvalue of the matrix
167*> pair (A,B), in one of the forms lambda = alpha/beta or
168*> mu = beta/alpha. Since either lambda or mu may overflow,
169*> they should not, in general, be computed.
170*> \endverbatim
171*>
172*> \param[in,out] QC
173*> \verbatim
174*> QC is DOUBLE PRECISION array, dimension (LDQC, NW)
175*> \endverbatim
176*>
177*> \param[in] LDQC
178*> \verbatim
179*> LDQC is INTEGER
180*> \endverbatim
181*>
182*> \param[in,out] ZC
183*> \verbatim
184*> ZC is DOUBLE PRECISION array, dimension (LDZC, NW)
185*> \endverbatim
186*>
187*> \param[in] LDZC
188*> \verbatim
189*> LDZ is INTEGER
190*> \endverbatim
191*>
192*> \param[out] WORK
193*> \verbatim
194*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
195*> On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
196*> \endverbatim
197*>
198*> \param[in] LWORK
199*> \verbatim
200*> LWORK is INTEGER
201*> The dimension of the array WORK. LWORK >= max(1,N).
202*>
203*> If LWORK = -1, then a workspace query is assumed; the routine
204*> only calculates the optimal size of the WORK array, returns
205*> this value as the first entry of the WORK array, and no error
206*> message related to LWORK is issued by XERBLA.
207*> \endverbatim
208*>
209*> \param[in] REC
210*> \verbatim
211*> REC is INTEGER
212*> REC indicates the current recursion level. Should be set
213*> to 0 on first call.
214*> \endverbatim
215*>
216*> \param[out] INFO
217*> \verbatim
218*> INFO is INTEGER
219*> = 0: successful exit
220*> < 0: if INFO = -i, the i-th argument had an illegal value
221*> \endverbatim
222*
223* Authors:
224* ========
225*
226*> \author Thijs Steel, KU Leuven
227*
228*> \date May 2020
229*
230*> \ingroup laqz3
231*>
232* =====================================================================
233 RECURSIVE SUBROUTINE dlaqz3( ILSCHUR, ILQ, ILZ, N, ILO, IHI,
234 $ NW,
235 $ A, LDA, B, LDB, Q, LDQ, Z, LDZ, NS,
236 $ ND, ALPHAR, ALPHAI, BETA, QC, LDQC,
237 $ ZC, LDZC, WORK, LWORK, REC, INFO )
238 IMPLICIT NONE
239
240* Arguments
241 LOGICAL, INTENT( IN ) :: ilschur, ilq, ilz
242 INTEGER, INTENT( IN ) :: n, ilo, ihi, nw, lda, ldb, ldq, ldz,
243 $ ldqc, ldzc, lwork, rec
244
245 DOUBLE PRECISION, INTENT( INOUT ) :: a( lda, * ), b( ldb, * ),
246 $ q( ldq, * ), z( ldz, * ), alphar( * ),
247 $ alphai( * ), beta( * )
248 INTEGER, INTENT( OUT ) :: ns, nd, info
249 DOUBLE PRECISION :: qc( ldqc, * ), zc( ldzc, * ), work( * )
250
251* Parameters
252 DOUBLE PRECISION :: zero, one, half
253 parameter( zero = 0.0d0, one = 1.0d0, half = 0.5d0 )
254
255* Local Scalars
256 LOGICAL :: bulge
257 INTEGER :: jw, kwtop, kwbot, istopm, istartm, k, k2, dtgexc_info,
258 $ ifst, ilst, lworkreq, qz_small_info
259 DOUBLE PRECISION :: s, smlnum, ulp, safmin, safmax, c1, s1, temp
260
261* External Functions
262 EXTERNAL :: xerbla, dtgexc, dlaqz0, dlacpy, dlaset,
264 DOUBLE PRECISION, EXTERNAL :: dlamch
265
266 info = 0
267
268* Set up deflation window
269 jw = min( nw, ihi-ilo+1 )
270 kwtop = ihi-jw+1
271 IF ( kwtop .EQ. ilo ) THEN
272 s = zero
273 ELSE
274 s = a( kwtop, kwtop-1 )
275 END IF
276
277* Determine required workspace
278 ifst = 1
279 ilst = jw
280 CALL dtgexc( .true., .true., jw, a, lda, b, ldb, qc, ldqc, zc,
281 $ ldzc, ifst, ilst, work, -1, dtgexc_info )
282 lworkreq = int( work( 1 ) )
283 CALL dlaqz0( 'S', 'V', 'V', jw, 1, jw, a( kwtop, kwtop ), lda,
284 $ b( kwtop, kwtop ), ldb, alphar, alphai, beta, qc,
285 $ ldqc, zc, ldzc, work, -1, rec+1, qz_small_info )
286 lworkreq = max( lworkreq, int( work( 1 ) )+2*jw**2 )
287 lworkreq = max( lworkreq, n*nw, 2*nw**2+n )
288 IF ( lwork .EQ.-1 ) THEN
289* workspace query, quick return
290 work( 1 ) = lworkreq
291 RETURN
292 ELSE IF ( lwork .LT. lworkreq ) THEN
293 info = -26
294 END IF
295
296 IF( info.NE.0 ) THEN
297 CALL xerbla( 'DLAQZ3', -info )
298 RETURN
299 END IF
300
301* Get machine constants
302 safmin = dlamch( 'SAFE MINIMUM' )
303 safmax = one/safmin
304 ulp = dlamch( 'PRECISION' )
305 smlnum = safmin*( dble( n )/ulp )
306
307 IF ( ihi .EQ. kwtop ) THEN
308* 1 by 1 deflation window, just try a regular deflation
309 alphar( kwtop ) = a( kwtop, kwtop )
310 alphai( kwtop ) = zero
311 beta( kwtop ) = b( kwtop, kwtop )
312 ns = 1
313 nd = 0
314 IF ( abs( s ) .LE. max( smlnum, ulp*abs( a( kwtop,
315 $ kwtop ) ) ) ) THEN
316 ns = 0
317 nd = 1
318 IF ( kwtop .GT. ilo ) THEN
319 a( kwtop, kwtop-1 ) = zero
320 END IF
321 END IF
322 END IF
323
324
325* Store window in case of convergence failure
326 CALL dlacpy( 'ALL', jw, jw, a( kwtop, kwtop ), lda, work, jw )
327 CALL dlacpy( 'ALL', jw, jw, b( kwtop, kwtop ), ldb,
328 $ work( jw**2+
329 $ 1 ), jw )
330
331* Transform window to real schur form
332 CALL dlaset( 'FULL', jw, jw, zero, one, qc, ldqc )
333 CALL dlaset( 'FULL', jw, jw, zero, one, zc, ldzc )
334 CALL dlaqz0( 'S', 'V', 'V', jw, 1, jw, a( kwtop, kwtop ), lda,
335 $ b( kwtop, kwtop ), ldb, alphar, alphai, beta, qc,
336 $ ldqc, zc, ldzc, work( 2*jw**2+1 ), lwork-2*jw**2,
337 $ rec+1, qz_small_info )
338
339 IF( qz_small_info .NE. 0 ) THEN
340* Convergence failure, restore the window and exit
341 nd = 0
342 ns = jw-qz_small_info
343 CALL dlacpy( 'ALL', jw, jw, work, jw, a( kwtop, kwtop ),
344 $ lda )
345 CALL dlacpy( 'ALL', jw, jw, work( jw**2+1 ), jw, b( kwtop,
346 $ kwtop ), ldb )
347 RETURN
348 END IF
349
350* Deflation detection loop
351 IF ( kwtop .EQ. ilo .OR. s .EQ. zero ) THEN
352 kwbot = kwtop-1
353 ELSE
354 kwbot = ihi
355 k = 1
356 k2 = 1
357 DO WHILE ( k .LE. jw )
358 bulge = .false.
359 IF ( kwbot-kwtop+1 .GE. 2 ) THEN
360 bulge = a( kwbot, kwbot-1 ) .NE. zero
361 END IF
362 IF ( bulge ) THEN
363
364* Try to deflate complex conjugate eigenvalue pair
365 temp = abs( a( kwbot, kwbot ) )+sqrt( abs( a( kwbot,
366 $ kwbot-1 ) ) )*sqrt( abs( a( kwbot-1, kwbot ) ) )
367 IF( temp .EQ. zero )THEN
368 temp = abs( s )
369 END IF
370 IF ( max( abs( s*qc( 1, kwbot-kwtop ) ), abs( s*qc( 1,
371 $ kwbot-kwtop+1 ) ) ) .LE. max( smlnum,
372 $ ulp*temp ) ) THEN
373* Deflatable
374 kwbot = kwbot-2
375 ELSE
376* Not deflatable, move out of the way
377 ifst = kwbot-kwtop+1
378 ilst = k2
379 CALL dtgexc( .true., .true., jw, a( kwtop, kwtop ),
380 $ lda, b( kwtop, kwtop ), ldb, qc, ldqc,
381 $ zc, ldzc, ifst, ilst, work, lwork,
382 $ dtgexc_info )
383 k2 = k2+2
384 END IF
385 k = k+2
386 ELSE
387
388* Try to deflate real eigenvalue
389 temp = abs( a( kwbot, kwbot ) )
390 IF( temp .EQ. zero ) THEN
391 temp = abs( s )
392 END IF
393 IF ( ( abs( s*qc( 1, kwbot-kwtop+1 ) ) ) .LE. max( ulp*
394 $ temp, smlnum ) ) THEN
395* Deflatable
396 kwbot = kwbot-1
397 ELSE
398* Not deflatable, move out of the way
399 ifst = kwbot-kwtop+1
400 ilst = k2
401 CALL dtgexc( .true., .true., jw, a( kwtop, kwtop ),
402 $ lda, b( kwtop, kwtop ), ldb, qc, ldqc,
403 $ zc, ldzc, ifst, ilst, work, lwork,
404 $ dtgexc_info )
405 k2 = k2+1
406 END IF
407
408 k = k+1
409
410 END IF
411 END DO
412 END IF
413
414* Store eigenvalues
415 nd = ihi-kwbot
416 ns = jw-nd
417 k = kwtop
418 DO WHILE ( k .LE. ihi )
419 bulge = .false.
420 IF ( k .LT. ihi ) THEN
421 IF ( a( k+1, k ) .NE. zero ) THEN
422 bulge = .true.
423 END IF
424 END IF
425 IF ( bulge ) THEN
426* 2x2 eigenvalue block
427 CALL dlag2( a( k, k ), lda, b( k, k ), ldb, safmin,
428 $ beta( k ), beta( k+1 ), alphar( k ),
429 $ alphar( k+1 ), alphai( k ) )
430 alphai( k+1 ) = -alphai( k )
431 k = k+2
432 ELSE
433* 1x1 eigenvalue block
434 alphar( k ) = a( k, k )
435 alphai( k ) = zero
436 beta( k ) = b( k, k )
437 k = k+1
438 END IF
439 END DO
440
441 IF ( kwtop .NE. ilo .AND. s .NE. zero ) THEN
442* Reflect spike back, this will create optimally packed bulges
443 a( kwtop:kwbot, kwtop-1 ) = a( kwtop, kwtop-1 )*qc( 1,
444 $ 1:jw-nd )
445 DO k = kwbot-1, kwtop, -1
446 CALL dlartg( a( k, kwtop-1 ), a( k+1, kwtop-1 ), c1, s1,
447 $ temp )
448 a( k, kwtop-1 ) = temp
449 a( k+1, kwtop-1 ) = zero
450 k2 = max( kwtop, k-1 )
451 CALL drot( ihi-k2+1, a( k, k2 ), lda, a( k+1, k2 ), lda,
452 $ c1,
453 $ s1 )
454 CALL drot( ihi-( k-1 )+1, b( k, k-1 ), ldb, b( k+1,
455 $ k-1 ),
456 $ ldb, c1, s1 )
457 CALL drot( jw, qc( 1, k-kwtop+1 ), 1, qc( 1,
458 $ k+1-kwtop+1 ),
459 $ 1, c1, s1 )
460 END DO
461
462* Chase bulges down
463 istartm = kwtop
464 istopm = ihi
465 k = kwbot-1
466 DO WHILE ( k .GE. kwtop )
467 IF ( ( k .GE. kwtop+1 ) .AND. a( k+1, k-1 ) .NE. zero ) THEN
468
469* Move double pole block down and remove it
470 DO k2 = k-1, kwbot-2
471 CALL dlaqz2( .true., .true., k2, kwtop, kwtop+jw-1,
472 $ kwbot, a, lda, b, ldb, jw, kwtop, qc,
473 $ ldqc, jw, kwtop, zc, ldzc )
474 END DO
475
476 k = k-2
477 ELSE
478
479* k points to single shift
480 DO k2 = k, kwbot-2
481
482* Move shift down
483 CALL dlartg( b( k2+1, k2+1 ), b( k2+1, k2 ), c1,
484 $ s1,
485 $ temp )
486 b( k2+1, k2+1 ) = temp
487 b( k2+1, k2 ) = zero
488 CALL drot( k2+2-istartm+1, a( istartm, k2+1 ), 1,
489 $ a( istartm, k2 ), 1, c1, s1 )
490 CALL drot( k2-istartm+1, b( istartm, k2+1 ), 1,
491 $ b( istartm, k2 ), 1, c1, s1 )
492 CALL drot( jw, zc( 1, k2+1-kwtop+1 ), 1, zc( 1,
493 $ k2-kwtop+1 ), 1, c1, s1 )
494
495 CALL dlartg( a( k2+1, k2 ), a( k2+2, k2 ), c1, s1,
496 $ temp )
497 a( k2+1, k2 ) = temp
498 a( k2+2, k2 ) = zero
499 CALL drot( istopm-k2, a( k2+1, k2+1 ), lda,
500 $ a( k2+2,
501 $ k2+1 ), lda, c1, s1 )
502 CALL drot( istopm-k2, b( k2+1, k2+1 ), ldb,
503 $ b( k2+2,
504 $ k2+1 ), ldb, c1, s1 )
505 CALL drot( jw, qc( 1, k2+1-kwtop+1 ), 1, qc( 1,
506 $ k2+2-kwtop+1 ), 1, c1, s1 )
507
508 END DO
509
510* Remove the shift
511 CALL dlartg( b( kwbot, kwbot ), b( kwbot, kwbot-1 ),
512 $ c1,
513 $ s1, temp )
514 b( kwbot, kwbot ) = temp
515 b( kwbot, kwbot-1 ) = zero
516 CALL drot( kwbot-istartm, b( istartm, kwbot ), 1,
517 $ b( istartm, kwbot-1 ), 1, c1, s1 )
518 CALL drot( kwbot-istartm+1, a( istartm, kwbot ), 1,
519 $ a( istartm, kwbot-1 ), 1, c1, s1 )
520 CALL drot( jw, zc( 1, kwbot-kwtop+1 ), 1, zc( 1,
521 $ kwbot-1-kwtop+1 ), 1, c1, s1 )
522
523 k = k-1
524 END IF
525 END DO
526
527 END IF
528
529* Apply Qc and Zc to rest of the matrix
530 IF ( ilschur ) THEN
531 istartm = 1
532 istopm = n
533 ELSE
534 istartm = ilo
535 istopm = ihi
536 END IF
537
538 IF ( istopm-ihi > 0 ) THEN
539 CALL dgemm( 'T', 'N', jw, istopm-ihi, jw, one, qc, ldqc,
540 $ a( kwtop, ihi+1 ), lda, zero, work, jw )
541 CALL dlacpy( 'ALL', jw, istopm-ihi, work, jw, a( kwtop,
542 $ ihi+1 ), lda )
543 CALL dgemm( 'T', 'N', jw, istopm-ihi, jw, one, qc, ldqc,
544 $ b( kwtop, ihi+1 ), ldb, zero, work, jw )
545 CALL dlacpy( 'ALL', jw, istopm-ihi, work, jw, b( kwtop,
546 $ ihi+1 ), ldb )
547 END IF
548 IF ( ilq ) THEN
549 CALL dgemm( 'N', 'N', n, jw, jw, one, q( 1, kwtop ), ldq,
550 $ qc,
551 $ ldqc, zero, work, n )
552 CALL dlacpy( 'ALL', n, jw, work, n, q( 1, kwtop ), ldq )
553 END IF
554
555 IF ( kwtop-1-istartm+1 > 0 ) THEN
556 CALL dgemm( 'N', 'N', kwtop-istartm, jw, jw, one,
557 $ a( istartm,
558 $ kwtop ), lda, zc, ldzc, zero, work,
559 $ kwtop-istartm )
560 CALL dlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
561 $ a( istartm, kwtop ), lda )
562 CALL dgemm( 'N', 'N', kwtop-istartm, jw, jw, one,
563 $ b( istartm,
564 $ kwtop ), ldb, zc, ldzc, zero, work,
565 $ kwtop-istartm )
566 CALL dlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
567 $ b( istartm, kwtop ), ldb )
568 END IF
569 IF ( ilz ) THEN
570 CALL dgemm( 'N', 'N', n, jw, jw, one, z( 1, kwtop ), ldz,
571 $ zc,
572 $ ldzc, zero, work, n )
573 CALL dlacpy( 'ALL', n, jw, work, n, z( 1, kwtop ), ldz )
574 END IF
575
576 END SUBROUTINE
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:101
subroutine dlag2(a, lda, b, ldb, safmin, scale1, scale2, wr1, wr2, wi)
DLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
Definition dlag2.f:154
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
recursive subroutine dlaqz0(wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb, alphar, alphai, beta, q, ldq, z, ldz, work, lwork, rec, info)
DLAQZ0
Definition dlaqz0.f:305
subroutine dlaqz2(ilq, ilz, k, istartm, istopm, ihi, a, lda, b, ldb, nq, qstart, q, ldq, nz, zstart, z, ldz)
DLAQZ2
Definition dlaqz2.f:173
recursive subroutine dlaqz3(ilschur, ilq, ilz, n, ilo, ihi, nw, a, lda, b, ldb, q, ldq, z, ldz, ns, nd, alphar, alphai, beta, qc, ldqc, zc, ldzc, work, lwork, rec, info)
DLAQZ3
Definition dlaqz3.f:238
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:108
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
subroutine dtgexc(wantq, wantz, n, a, lda, b, ldb, q, ldq, z, ldz, ifst, ilst, work, lwork, info)
DTGEXC
Definition dtgexc.f:218