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