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