OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
zlaqz3.f
Go to the documentation of this file.
1*> \brief \b ZLAQZ3
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZLAQZ3 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ZLAQZ3.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ZLAQZ3.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ZLAQZ3.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZLAQZ3( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NSHIFTS,
22* $ NBLOCK_DESIRED, ALPHA, BETA, A, LDA, B, LDB, Q, LDQ, Z, LDZ,
23* $ QC, LDQC, ZC, LDZC, WORK, LWORK, INFO )
24* IMPLICIT NONE
25*
26* Function arguments
27* LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
28* INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK,
29* $ NSHIFTS, NBLOCK_DESIRED, LDQC, LDZC
30*
31* COMPLEX*16, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ,
32* $ * ), Z( LDZ, * ), QC( LDQC, * ), ZC( LDZC, * ), WORK( * ),
33* $ ALPHA( * ), BETA( * )
34*
35* INTEGER, INTENT( OUT ) :: INFO
36* ..
37*
38*
39*> \par Purpose:
40* =============
41*>
42*> \verbatim
43*>
44*> ZLAQZ3 Executes a single multishift QZ sweep
45*> \endverbatim
46*
47* Arguments:
48* ==========
49*
50*> \param[in] ILSCHUR
51*> \verbatim
52*> ILSCHUR is LOGICAL
53*> Determines whether or not to update the full Schur form
54*> \endverbatim
55*>
56*> \param[in] ILQ
57*> \verbatim
58*> ILQ is LOGICAL
59*> Determines whether or not to update the matrix Q
60*> \endverbatim
61*>
62*> \param[in] ILZ
63*> \verbatim
64*> ILZ is LOGICAL
65*> Determines whether or not to update the matrix Z
66*> \endverbatim
67*>
68*> \param[in] N
69*> \verbatim
70*> N is INTEGER
71*> The order of the matrices A, B, Q, and Z. N >= 0.
72*> \endverbatim
73*>
74*> \param[in] ILO
75*> \verbatim
76*> ILO is INTEGER
77*> \endverbatim
78*>
79*> \param[in] IHI
80*> \verbatim
81*> IHI is INTEGER
82*> \endverbatim
83*>
84*> \param[in] NSHIFTS
85*> \verbatim
86*> NSHIFTS is INTEGER
87*> The desired number of shifts to use
88*> \endverbatim
89*>
90*> \param[in] NBLOCK_DESIRED
91*> \verbatim
92*> NBLOCK_DESIRED is INTEGER
93*> The desired size of the computational windows
94*> \endverbatim
95*>
96*> \param[in] ALPHA
97*> \verbatim
98*> ALPHA is COMPLEX*16 array. SR contains
99*> the alpha parts of the shifts to use.
100*> \endverbatim
101*>
102*> \param[in] BETA
103*> \verbatim
104*> BETA is COMPLEX*16 array. SS contains
105*> the scale of the shifts to use.
106*> \endverbatim
107*>
108*> \param[in,out] A
109*> \verbatim
110*> A is COMPLEX*16 array, dimension (LDA, N)
111*> \endverbatim
112*>
113*> \param[in] LDA
114*> \verbatim
115*> LDA is INTEGER
116*> The leading dimension of the array A. LDA >= max( 1, N ).
117*> \endverbatim
118*>
119*> \param[in,out] B
120*> \verbatim
121*> B is COMPLEX*16 array, dimension (LDB, N)
122*> \endverbatim
123*>
124*> \param[in] LDB
125*> \verbatim
126*> LDB is INTEGER
127*> The leading dimension of the array B. LDB >= max( 1, N ).
128*> \endverbatim
129*>
130*> \param[in,out] Q
131*> \verbatim
132*> Q is COMPLEX*16 array, dimension (LDQ, N)
133*> \endverbatim
134*>
135*> \param[in] LDQ
136*> \verbatim
137*> LDQ is INTEGER
138*> \endverbatim
139*>
140*> \param[in,out] Z
141*> \verbatim
142*> Z is COMPLEX*16 array, dimension (LDZ, N)
143*> \endverbatim
144*>
145*> \param[in] LDZ
146*> \verbatim
147*> LDZ is INTEGER
148*> \endverbatim
149*>
150*> \param[in,out] QC
151*> \verbatim
152*> QC is COMPLEX*16 array, dimension (LDQC, NBLOCK_DESIRED)
153*> \endverbatim
154*>
155*> \param[in] LDQC
156*> \verbatim
157*> LDQC is INTEGER
158*> \endverbatim
159*>
160*> \param[in,out] ZC
161*> \verbatim
162*> ZC is COMPLEX*16 array, dimension (LDZC, NBLOCK_DESIRED)
163*> \endverbatim
164*>
165*> \param[in] LDZC
166*> \verbatim
167*> LDZ is INTEGER
168*> \endverbatim
169*>
170*> \param[out] WORK
171*> \verbatim
172*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
173*> On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
174*> \endverbatim
175*>
176*> \param[in] LWORK
177*> \verbatim
178*> LWORK is INTEGER
179*> The dimension of the array WORK. LWORK >= max(1,N).
180*>
181*> If LWORK = -1, then a workspace query is assumed; the routine
182*> only calculates the optimal size of the WORK array, returns
183*> this value as the first entry of the WORK array, and no error
184*> message related to LWORK is issued by XERBLA.
185*> \endverbatim
186*>
187*> \param[out] INFO
188*> \verbatim
189*> INFO is INTEGER
190*> = 0: successful exit
191*> < 0: if INFO = -i, the i-th argument had an illegal value
192*> \endverbatim
193*
194* Authors:
195* ========
196*
197*> \author Thijs Steel, KU Leuven
198*
199*> \date May 2020
200*
201*> \ingroup complex16GEcomputational
202*>
203* =====================================================================
204 SUBROUTINE zlaqz3( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NSHIFTS,
205 $ NBLOCK_DESIRED, ALPHA, BETA, A, LDA, B, LDB,
206 $ Q, LDQ, Z, LDZ, QC, LDQC, ZC, LDZC, WORK,
207 $ LWORK, INFO )
208 IMPLICIT NONE
209
210* Function arguments
211 LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
212 INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK,
213 $ NSHIFTS, NBLOCK_DESIRED, LDQC, LDZC
214
215 COMPLEX*16, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ,
216 $ * ), Z( LDZ, * ), QC( LDQC, * ), ZC( LDZC, * ), WORK( * ),
217 $ ALPHA( * ), BETA( * )
218
219 INTEGER, INTENT( OUT ) :: INFO
220
221* Parameters
222 COMPLEX*16 CZERO, CONE
223 PARAMETER ( CZERO = ( 0.0d+0, 0.0d+0 ), cone = ( 1.0d+0,
224 $ 0.0d+0 ) )
225 DOUBLE PRECISION :: ZERO, ONE, HALF
226 parameter( zero = 0.0d0, one = 1.0d0, half = 0.5d0 )
227
228* Local scalars
229 INTEGER :: I, J, NS, ISTARTM, ISTOPM, SHEIGHT, SWIDTH, K, NP,
230 $ ISTARTB, ISTOPB, ISHIFT, NBLOCK, NPOS
231 DOUBLE PRECISION :: SAFMIN, SAFMAX, C, SCALE
232 COMPLEX*16 :: TEMP, TEMP2, TEMP3, S
233
234* External Functions
235 EXTERNAL :: xerbla, dlabad, zlaset, zlartg, zrot, zlaqz1, zgemm,
236 $ zlacpy
237 DOUBLE PRECISION, EXTERNAL :: DLAMCH
238
239 info = 0
240 IF ( nblock_desired .LT. nshifts+1 ) THEN
241 info = -8
242 END IF
243 IF ( lwork .EQ.-1 ) THEN
244* workspace query, quick return
245 work( 1 ) = n*nblock_desired
246 RETURN
247 ELSE IF ( lwork .LT. n*nblock_desired ) THEN
248 info = -25
249 END IF
250
251 IF( info.NE.0 ) THEN
252 CALL xerbla( 'ZLAQZ3', -info )
253 RETURN
254 END IF
255
256*
257* Executable statements
258*
259
260* Get machine constants
261 safmin = dlamch( 'SAFE MINIMUM' )
262 safmax = one/safmin
263 CALL dlabad( safmin, safmax )
264
265 IF ( ilo .GE. ihi ) THEN
266 RETURN
267 END IF
268
269 IF ( ilschur ) THEN
270 istartm = 1
271 istopm = n
272 ELSE
273 istartm = ilo
274 istopm = ihi
275 END IF
276
277 ns = nshifts
278 npos = max( nblock_desired-ns, 1 )
279
280
281* The following block introduces the shifts and chases
282* them down one by one just enough to make space for
283* the other shifts. The near-the-diagonal block is
284* of size (ns+1) x ns.
285
286 CALL zlaset( 'FULL', ns+1, ns+1, czero, cone, qc, ldqc )
287 CALL zlaset( 'FULL', ns, ns, czero, cone, zc, ldzc )
288
289 DO i = 1, ns
290* Introduce the shift
291 scale = sqrt( abs( alpha( i ) ) ) * sqrt( abs( beta( i ) ) )
292 IF( scale .GE. safmin .AND. scale .LE. safmax ) THEN
293 alpha( i ) = alpha( i )/scale
294 beta( i ) = beta( i )/scale
295 END IF
296
297 temp2 = beta( i )*a( ilo, ilo )-alpha( i )*b( ilo, ilo )
298 temp3 = beta( i )*a( ilo+1, ilo )
299
300 IF ( abs( temp2 ) .GT. safmax .OR.
301 $ abs( temp3 ) .GT. safmax ) THEN
302 temp2 = cone
303 temp3 = czero
304 END IF
305
306 CALL zlartg( temp2, temp3, c, s, temp )
307 CALL zrot( ns, a( ilo, ilo ), lda, a( ilo+1, ilo ), lda, c,
308 $ s )
309 CALL zrot( ns, b( ilo, ilo ), ldb, b( ilo+1, ilo ), ldb, c,
310 $ s )
311 CALL zrot( ns+1, qc( 1, 1 ), 1, qc( 1, 2 ), 1, c,
312 $ dconjg( s ) )
313
314* Chase the shift down
315 DO j = 1, ns-i
316
317 CALL zlaqz1( .true., .true., j, 1, ns, ihi-ilo+1, a( ilo,
318 $ ilo ), lda, b( ilo, ilo ), ldb, ns+1, 1, qc,
319 $ ldqc, ns, 1, zc, ldzc )
320
321 END DO
322
323 END DO
324
325* Update the rest of the pencil
326
327* Update A(ilo:ilo+ns,ilo+ns:istopm) and B(ilo:ilo+ns,ilo+ns:istopm)
328* from the left with Qc(1:ns+1,1:ns+1)'
329 sheight = ns+1
330 swidth = istopm-( ilo+ns )+1
331 IF ( swidth > 0 ) THEN
332 CALL zgemm( 'C', 'N', sheight, swidth, sheight, cone, qc, ldqc,
333 $ a( ilo, ilo+ns ), lda, czero, work, sheight )
334 CALL zlacpy( 'ALL', sheight, swidth, work, sheight, a( ilo,
335 $ ilo+ns ), lda )
336 CALL zgemm( 'C', 'N', sheight, swidth, sheight, cone, qc, ldqc,
337 $ b( ilo, ilo+ns ), ldb, czero, work, sheight )
338 CALL zlacpy( 'ALL', sheight, swidth, work, sheight, b( ilo,
339 $ ilo+ns ), ldb )
340 END IF
341 IF ( ilq ) THEN
342 CALL zgemm( 'N', 'N', n, sheight, sheight, cone, q( 1, ilo ),
343 $ ldq, qc, ldqc, czero, work, n )
344 CALL zlacpy( 'ALL', n, sheight, work, n, q( 1, ilo ), ldq )
345 END IF
346
347* Update A(istartm:ilo-1,ilo:ilo+ns-1) and B(istartm:ilo-1,ilo:ilo+ns-1)
348* from the right with Zc(1:ns,1:ns)
349 sheight = ilo-1-istartm+1
350 swidth = ns
351 IF ( sheight > 0 ) THEN
352 CALL zgemm( 'N', 'N', sheight, swidth, swidth, cone,
353 $ a( istartm, ilo ), lda, zc, ldzc, czero, work,
354 $ sheight )
355 CALL zlacpy( 'ALL', sheight, swidth, work, sheight, a( istartm,
356 $ ilo ), lda )
357 CALL zgemm( 'N', 'N', sheight, swidth, swidth, cone,
358 $ b( istartm, ilo ), ldb, zc, ldzc, czero, work,
359 $ sheight )
360 CALL zlacpy( 'ALL', sheight, swidth, work, sheight, b( istartm,
361 $ ilo ), ldb )
362 END IF
363 IF ( ilz ) THEN
364 CALL zgemm( 'N', 'N', n, swidth, swidth, cone, z( 1, ilo ),
365 $ ldz, zc, ldzc, czero, work, n )
366 CALL zlacpy( 'ALL', n, swidth, work, n, z( 1, ilo ), ldz )
367 END IF
368
369* The following block chases the shifts down to the bottom
370* right block. If possible, a shift is moved down npos
371* positions at a time
372
373 k = ilo
374 DO WHILE ( k < ihi-ns )
375 np = min( ihi-ns-k, npos )
376* Size of the near-the-diagonal block
377 nblock = ns+np
378* istartb points to the first row we will be updating
379 istartb = k+1
380* istopb points to the last column we will be updating
381 istopb = k+nblock-1
382
383 CALL zlaset( 'FULL', ns+np, ns+np, czero, cone, qc, ldqc )
384 CALL zlaset( 'FULL', ns+np, ns+np, czero, cone, zc, ldzc )
385
386* Near the diagonal shift chase
387 DO i = ns-1, 0, -1
388 DO j = 0, np-1
389* Move down the block with index k+i+j, updating
390* the (ns+np x ns+np) block:
391* (k:k+ns+np,k:k+ns+np-1)
392 CALL zlaqz1( .true., .true., k+i+j, istartb, istopb, ihi,
393 $ a, lda, b, ldb, nblock, k+1, qc, ldqc,
394 $ nblock, k, zc, ldzc )
395 END DO
396 END DO
397
398* Update rest of the pencil
399
400* Update A(k+1:k+ns+np, k+ns+np:istopm) and
401* B(k+1:k+ns+np, k+ns+np:istopm)
402* from the left with Qc(1:ns+np,1:ns+np)'
403 sheight = ns+np
404 swidth = istopm-( k+ns+np )+1
405 IF ( swidth > 0 ) THEN
406 CALL zgemm( 'C', 'N', sheight, swidth, sheight, cone, qc,
407 $ ldqc, a( k+1, k+ns+np ), lda, czero, work,
408 $ sheight )
409 CALL zlacpy( 'ALL', sheight, swidth, work, sheight, a( k+1,
410 $ k+ns+np ), lda )
411 CALL zgemm( 'C', 'N', sheight, swidth, sheight, cone, qc,
412 $ ldqc, b( k+1, k+ns+np ), ldb, czero, work,
413 $ sheight )
414 CALL zlacpy( 'ALL', sheight, swidth, work, sheight, b( k+1,
415 $ k+ns+np ), ldb )
416 END IF
417 IF ( ilq ) THEN
418 CALL zgemm( 'N', 'N', n, nblock, nblock, cone, q( 1, k+1 ),
419 $ ldq, qc, ldqc, czero, work, n )
420 CALL zlacpy( 'ALL', n, nblock, work, n, q( 1, k+1 ), ldq )
421 END IF
422
423* Update A(istartm:k,k:k+ns+npos-1) and B(istartm:k,k:k+ns+npos-1)
424* from the right with Zc(1:ns+np,1:ns+np)
425 sheight = k-istartm+1
426 swidth = nblock
427 IF ( sheight > 0 ) THEN
428 CALL zgemm( 'N', 'N', sheight, swidth, swidth, cone,
429 $ a( istartm, k ), lda, zc, ldzc, czero, work,
430 $ sheight )
431 CALL zlacpy( 'ALL', sheight, swidth, work, sheight,
432 $ a( istartm, k ), lda )
433 CALL zgemm( 'N', 'N', sheight, swidth, swidth, cone,
434 $ b( istartm, k ), ldb, zc, ldzc, czero, work,
435 $ sheight )
436 CALL zlacpy( 'ALL', sheight, swidth, work, sheight,
437 $ b( istartm, k ), ldb )
438 END IF
439 IF ( ilz ) THEN
440 CALL zgemm( 'N', 'N', n, nblock, nblock, cone, z( 1, k ),
441 $ ldz, zc, ldzc, czero, work, n )
442 CALL zlacpy( 'ALL', n, nblock, work, n, z( 1, k ), ldz )
443 END IF
444
445 k = k+np
446
447 END DO
448
449* The following block removes the shifts from the bottom right corner
450* one by one. Updates are initially applied to A(ihi-ns+1:ihi,ihi-ns:ihi).
451
452 CALL zlaset( 'full', NS, NS, CZERO, CONE, QC, LDQC )
453 CALL ZLASET( 'full', NS+1, NS+1, CZERO, CONE, ZC, LDZC )
454
455* istartb points to the first row we will be updating
456 ISTARTB = IHI-NS+1
457* istopb points to the last column we will be updating
458 ISTOPB = IHI
459
460 DO I = 1, NS
461* Chase the shift down to the bottom right corner
462 DO ISHIFT = IHI-I, IHI-1
463 CALL ZLAQZ1( .TRUE., .TRUE., ISHIFT, ISTARTB, ISTOPB, IHI,
464 $ A, LDA, B, LDB, NS, IHI-NS+1, QC, LDQC, NS+1,
465 $ IHI-NS, ZC, LDZC )
466 END DO
467
468 END DO
469
470* Update rest of the pencil
471
472* Update A(ihi-ns+1:ihi, ihi+1:istopm)
473* from the left with Qc(1:ns,1:ns)'
474 SHEIGHT = NS
475 SWIDTH = ISTOPM-( IHI+1 )+1
476 IF ( SWIDTH > 0 ) THEN
477 CALL ZGEMM( 'c', 'n', SHEIGHT, SWIDTH, SHEIGHT, CONE, QC, LDQC,
478 $ A( IHI-NS+1, IHI+1 ), LDA, CZERO, WORK, SHEIGHT )
479 CALL ZLACPY( 'all', SHEIGHT, SWIDTH, WORK, SHEIGHT,
480 $ A( IHI-NS+1, IHI+1 ), LDA )
481 CALL ZGEMM( 'c', 'n', SHEIGHT, SWIDTH, SHEIGHT, CONE, QC, LDQC,
482 $ B( IHI-NS+1, IHI+1 ), LDB, CZERO, WORK, SHEIGHT )
483 CALL ZLACPY( 'all', SHEIGHT, SWIDTH, WORK, SHEIGHT,
484 $ B( IHI-NS+1, IHI+1 ), LDB )
485 END IF
486 IF ( ILQ ) THEN
487 CALL ZGEMM( 'n', 'n', n, ns, ns, cone, q( 1, ihi-ns+1 ), ldq,
488 $ qc, ldqc, czero, work, n )
489 CALL zlacpy( 'ALL', n, ns, work, n, q( 1, ihi-ns+1 ), ldq )
490 END IF
491
492* Update A(istartm:ihi-ns,ihi-ns:ihi)
493* from the right with Zc(1:ns+1,1:ns+1)
494 sheight = ihi-ns-istartm+1
495 swidth = ns+1
496 IF ( sheight > 0 ) THEN
497 CALL zgemm( 'N', 'N', sheight, swidth, swidth, cone,
498 $ a( istartm, ihi-ns ), lda, zc, ldzc, czero, work,
499 $ sheight )
500 CALL zlacpy( 'ALL', sheight, swidth, work, sheight, a( istartm,
501 $ ihi-ns ), lda )
502 CALL zgemm( 'N', 'N', sheight, swidth, swidth, cone,
503 $ b( istartm, ihi-ns ), ldb, zc, ldzc, czero, work,
504 $ sheight )
505 CALL zlacpy( 'ALL', sheight, swidth, work, sheight, b( istartm,
506 $ ihi-ns ), ldb )
507 END IF
508 IF ( ilz ) THEN
509 CALL zgemm( 'N', 'N', n, ns+1, ns+1, cone, z( 1, ihi-ns ), ldz,
510 $ zc, ldzc, czero, work, n )
511 CALL zlacpy( 'ALL', n, ns+1, work, n, z( 1, ihi-ns ), ldz )
512 END IF
513
514 END SUBROUTINE
zlartg
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
Definition zlartg.f90:118
dlabad
subroutine dlabad(small, large)
DLABAD
Definition dlabad.f:74
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
zlaqz1
subroutine zlaqz1(ilq, ilz, k, istartm, istopm, ihi, a, lda, b, ldb, nq, qstart, q, ldq, nz, zstart, z, ldz)
ZLAQZ1
Definition zlaqz1.f:173
zlaqz3
subroutine zlaqz3(ilschur, ilq, ilz, n, ilo, ihi, nshifts, nblock_desired, alpha, beta, a, lda, b, ldb, q, ldq, z, ldz, qc, ldqc, zc, ldzc, work, lwork, info)
ZLAQZ3
Definition zlaqz3.f:208
zlacpy
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:103
zrot
subroutine zrot(n, cx, incx, cy, incy, c, s)
ZROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
Definition zrot.f:103
zlaset
subroutine zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition zlaset.f:106
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:187
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21