LAPACK 3.11.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 doubleGEcomputational
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
240 info = 0
241 IF ( nblock_desired .LT. nshifts+1 ) THEN
242 info = -8
243 END IF
244 IF ( lwork .EQ.-1 ) THEN
245* workspace query, quick return
246 work( 1 ) = n*nblock_desired
247 RETURN
248 ELSE IF ( lwork .LT. n*nblock_desired ) THEN
249 info = -25
250 END IF
251
252 IF( info.NE.0 ) THEN
253 CALL xerbla( 'SLAQZ4', -info )
254 RETURN
255 END IF
256
257* Executable statements
258
259 IF ( nshifts .LT. 2 ) THEN
260 RETURN
261 END IF
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* Shuffle shifts into pairs of real shifts and pairs
276* of complex conjugate shifts assuming complex
277* conjugate shifts are already adjacent to one
278* another
279
280 DO i = 1, nshifts-2, 2
281 IF( si( i ).NE.-si( i+1 ) ) THEN
282*
283 swap = sr( i )
284 sr( i ) = sr( i+1 )
285 sr( i+1 ) = sr( i+2 )
286 sr( i+2 ) = swap
287
288 swap = si( i )
289 si( i ) = si( i+1 )
290 si( i+1 ) = si( i+2 )
291 si( i+2 ) = swap
292
293 swap = ss( i )
294 ss( i ) = ss( i+1 )
295 ss( i+1 ) = ss( i+2 )
296 ss( i+2 ) = swap
297 END IF
298 END DO
299
300* NSHFTS is supposed to be even, but if it is odd,
301* then simply reduce it by one. The shuffle above
302* ensures that the dropped shift is real and that
303* the remaining shifts are paired.
304
305 ns = nshifts-mod( nshifts, 2 )
306 npos = max( nblock_desired-ns, 1 )
307
308* The following block introduces the shifts and chases
309* them down one by one just enough to make space for
310* the other shifts. The near-the-diagonal block is
311* of size (ns+1) x ns.
312
313 CALL slaset( 'FULL', ns+1, ns+1, zero, one, qc, ldqc )
314 CALL slaset( 'FULL', ns, ns, zero, one, zc, ldzc )
315
316 DO i = 1, ns, 2
317* Introduce the shift
318 CALL slaqz1( a( ilo, ilo ), lda, b( ilo, ilo ), ldb, 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, c1,
326 $ s1 )
327 CALL srot( ns, a( ilo, ilo ), lda, a( ilo+1, ilo ), lda, c2,
328 $ s2 )
329 CALL srot( ns, b( ilo+1, ilo ), ldb, b( ilo+2, ilo ), ldb, c1,
330 $ s1 )
331 CALL srot( ns, b( ilo, ilo ), ldb, b( ilo+1, ilo ), ldb, c2,
332 $ s2 )
333 CALL srot( ns+1, qc( 1, 2 ), 1, qc( 1, 3 ), 1, c1, s1 )
334 CALL srot( ns+1, qc( 1, 1 ), 1, qc( 1, 2 ), 1, c2, s2 )
335
336* Chase the shift down
337 DO j = 1, ns-1-i
338
339 CALL slaqz2( .true., .true., j, 1, ns, ihi-ilo+1, a( ilo,
340 $ ilo ), lda, b( ilo, ilo ), ldb, ns+1, 1, qc,
341 $ ldqc, ns, 1, zc, ldzc )
342
343 END DO
344
345 END DO
346
347* Update the rest of the pencil
348
349* Update A(ilo:ilo+ns,ilo+ns:istopm) and B(ilo:ilo+ns,ilo+ns:istopm)
350* from the left with Qc(1:ns+1,1:ns+1)'
351 sheight = ns+1
352 swidth = istopm-( ilo+ns )+1
353 IF ( swidth > 0 ) THEN
354 CALL sgemm( 'T', 'N', sheight, swidth, sheight, one, qc, ldqc,
355 $ a( ilo, ilo+ns ), lda, zero, work, sheight )
356 CALL slacpy( 'ALL', sheight, swidth, work, sheight, a( ilo,
357 $ ilo+ns ), lda )
358 CALL sgemm( 'T', 'N', sheight, swidth, sheight, one, qc, ldqc,
359 $ b( ilo, ilo+ns ), ldb, zero, work, sheight )
360 CALL slacpy( 'ALL', sheight, swidth, work, sheight, b( ilo,
361 $ ilo+ns ), ldb )
362 END IF
363 IF ( ilq ) THEN
364 CALL sgemm( 'N', 'N', n, sheight, sheight, one, q( 1, ilo ),
365 $ ldq, qc, ldqc, zero, work, n )
366 CALL slacpy( 'ALL', n, sheight, work, n, q( 1, ilo ), ldq )
367 END IF
368
369* Update A(istartm:ilo-1,ilo:ilo+ns-1) and B(istartm:ilo-1,ilo:ilo+ns-1)
370* from the right with Zc(1:ns,1:ns)
371 sheight = ilo-1-istartm+1
372 swidth = ns
373 IF ( sheight > 0 ) THEN
374 CALL sgemm( 'N', 'N', sheight, swidth, swidth, one, a( istartm,
375 $ ilo ), lda, zc, ldzc, zero, work, sheight )
376 CALL slacpy( 'ALL', sheight, swidth, work, sheight, a( istartm,
377 $ ilo ), lda )
378 CALL sgemm( 'N', 'N', sheight, swidth, swidth, one, b( istartm,
379 $ ilo ), ldb, zc, ldzc, zero, work, sheight )
380 CALL slacpy( 'ALL', sheight, swidth, work, sheight, b( istartm,
381 $ ilo ), ldb )
382 END IF
383 IF ( ilz ) THEN
384 CALL sgemm( 'N', 'N', n, swidth, swidth, one, z( 1, ilo ), ldz,
385 $ zc, ldzc, zero, work, n )
386 CALL slacpy( 'ALL', n, swidth, work, n, z( 1, ilo ), ldz )
387 END IF
388
389* The following block chases the shifts down to the bottom
390* right block. If possible, a shift is moved down npos
391* positions at a time
392
393 k = ilo
394 DO WHILE ( k < ihi-ns )
395 np = min( ihi-ns-k, npos )
396* Size of the near-the-diagonal block
397 nblock = ns+np
398* istartb points to the first row we will be updating
399 istartb = k+1
400* istopb points to the last column we will be updating
401 istopb = k+nblock-1
402
403 CALL slaset( 'FULL', ns+np, ns+np, zero, one, qc, ldqc )
404 CALL slaset( 'FULL', ns+np, ns+np, zero, one, zc, ldzc )
405
406* Near the diagonal shift chase
407 DO i = ns-1, 0, -2
408 DO j = 0, np-1
409* Move down the block with index k+i+j-1, updating
410* the (ns+np x ns+np) block:
411* (k:k+ns+np,k:k+ns+np-1)
412 CALL slaqz2( .true., .true., k+i+j-1, istartb, istopb,
413 $ ihi, a, lda, b, ldb, nblock, k+1, qc, ldqc,
414 $ nblock, k, zc, ldzc )
415 END DO
416 END DO
417
418* Update rest of the pencil
419
420* Update A(k+1:k+ns+np, k+ns+np:istopm) and
421* B(k+1:k+ns+np, k+ns+np:istopm)
422* from the left with Qc(1:ns+np,1:ns+np)'
423 sheight = ns+np
424 swidth = istopm-( k+ns+np )+1
425 IF ( swidth > 0 ) THEN
426 CALL sgemm( 'T', 'N', sheight, swidth, sheight, one, qc,
427 $ ldqc, a( k+1, k+ns+np ), lda, zero, work,
428 $ sheight )
429 CALL slacpy( 'ALL', sheight, swidth, work, sheight, a( k+1,
430 $ k+ns+np ), lda )
431 CALL sgemm( 'T', 'N', sheight, swidth, sheight, one, qc,
432 $ ldqc, b( k+1, k+ns+np ), ldb, zero, work,
433 $ sheight )
434 CALL slacpy( 'ALL', sheight, swidth, work, sheight, b( k+1,
435 $ k+ns+np ), ldb )
436 END IF
437 IF ( ilq ) THEN
438 CALL sgemm( 'N', 'N', n, nblock, nblock, one, q( 1, k+1 ),
439 $ ldq, qc, ldqc, zero, work, n )
440 CALL slacpy( 'ALL', n, nblock, work, n, q( 1, k+1 ), ldq )
441 END IF
442
443* Update A(istartm:k,k:k+ns+npos-1) and B(istartm:k,k:k+ns+npos-1)
444* from the right with Zc(1:ns+np,1:ns+np)
445 sheight = k-istartm+1
446 swidth = nblock
447 IF ( sheight > 0 ) THEN
448 CALL sgemm( 'N', 'N', sheight, swidth, swidth, one,
449 $ a( istartm, k ), lda, zc, ldzc, zero, work,
450 $ sheight )
451 CALL slacpy( 'ALL', sheight, swidth, work, sheight,
452 $ a( istartm, k ), lda )
453 CALL sgemm( 'N', 'N', sheight, swidth, swidth, one,
454 $ b( istartm, k ), ldb, zc, ldzc, zero, work,
455 $ sheight )
456 CALL slacpy( 'ALL', sheight, swidth, work, sheight,
457 $ b( istartm, k ), ldb )
458 END IF
459 IF ( ilz ) THEN
460 CALL sgemm( 'N', 'N', n, nblock, nblock, one, z( 1, k ),
461 $ ldz, zc, ldzc, zero, work, n )
462 CALL slacpy( 'ALL', n, nblock, work, n, z( 1, k ), ldz )
463 END IF
464
465 k = k+np
466
467 END DO
468
469* The following block removes the shifts from the bottom right corner
470* one by one. Updates are initially applied to A(ihi-ns+1:ihi,ihi-ns:ihi).
471
472 CALL slaset( 'FULL', ns, ns, zero, one, qc, ldqc )
473 CALL slaset( 'FULL', ns+1, ns+1, zero, one, zc, ldzc )
474
475* istartb points to the first row we will be updating
476 istartb = ihi-ns+1
477* istopb points to the last column we will be updating
478 istopb = ihi
479
480 DO i = 1, ns, 2
481* Chase the shift down to the bottom right corner
482 DO ishift = ihi-i-1, ihi-2
483 CALL slaqz2( .true., .true., ishift, istartb, istopb, ihi,
484 $ a, lda, b, ldb, ns, ihi-ns+1, qc, ldqc, ns+1,
485 $ ihi-ns, zc, ldzc )
486 END DO
487
488 END DO
489
490* Update rest of the pencil
491
492* Update A(ihi-ns+1:ihi, ihi+1:istopm)
493* from the left with Qc(1:ns,1:ns)'
494 sheight = ns
495 swidth = istopm-( ihi+1 )+1
496 IF ( swidth > 0 ) THEN
497 CALL sgemm( 'T', 'N', sheight, swidth, sheight, one, qc, ldqc,
498 $ a( ihi-ns+1, ihi+1 ), lda, zero, work, sheight )
499 CALL slacpy( 'ALL', sheight, swidth, work, sheight,
500 $ a( ihi-ns+1, ihi+1 ), lda )
501 CALL sgemm( 'T', 'N', sheight, swidth, sheight, one, qc, ldqc,
502 $ b( ihi-ns+1, ihi+1 ), ldb, zero, work, sheight )
503 CALL slacpy( 'ALL', sheight, swidth, work, sheight,
504 $ b( ihi-ns+1, ihi+1 ), ldb )
505 END IF
506 IF ( ilq ) THEN
507 CALL sgemm( 'N', 'N', n, ns, ns, one, q( 1, ihi-ns+1 ), ldq,
508 $ qc, ldqc, zero, work, n )
509 CALL slacpy( 'ALL', n, ns, work, n, q( 1, ihi-ns+1 ), ldq )
510 END IF
511
512* Update A(istartm:ihi-ns,ihi-ns:ihi)
513* from the right with Zc(1:ns+1,1:ns+1)
514 sheight = ihi-ns-istartm+1
515 swidth = ns+1
516 IF ( sheight > 0 ) THEN
517 CALL sgemm( 'N', 'N', sheight, swidth, swidth, one, a( istartm,
518 $ ihi-ns ), lda, zc, ldzc, zero, work, sheight )
519 CALL slacpy( 'ALL', sheight, swidth, work, sheight, a( istartm,
520 $ ihi-ns ), lda )
521 CALL sgemm( 'N', 'N', sheight, swidth, swidth, one, b( istartm,
522 $ ihi-ns ), ldb, zc, ldzc, zero, work, sheight )
523 CALL slacpy( 'ALL', sheight, swidth, work, sheight, b( istartm,
524 $ ihi-ns ), ldb )
525 END IF
526 IF ( ilz ) THEN
527 CALL sgemm( 'N', 'N', n, ns+1, ns+1, one, z( 1, ihi-ns ), ldz, zc,
528 $ ldzc, zero, work, n )
529 CALL slacpy( 'ALL', n, ns+1, work, n, z( 1, ihi-ns ), ldz )
530 END IF
531
532 END SUBROUTINE
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
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
slartg
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f90:111
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
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
slaqz1
subroutine slaqz1(A, LDA, B, LDB, SR1, SR2, SI, BETA1, BETA2, V)
SLAQZ1
Definition: slaqz1.f:127
srot
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:92
sgemm
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187