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