1*> \brief \b SLAQZ3
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SLAQZ3 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaqz3.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaqz3.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaqz3.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*      SUBROUTINE SLAQZ3( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NW, A, LDA, B,
22*     $    LDB, Q, LDQ, Z, LDZ, NS, ND, ALPHAR, ALPHAI, BETA, QC, LDQC,
23*     $    ZC, LDZC, WORK, LWORK, REC, INFO )
24*      IMPLICIT NONE
25*
26*      Arguments
27*      LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
28*      INTEGER, INTENT( IN ) :: N, ILO, IHI, NW, LDA, LDB, LDQ, LDZ,
29*     $    LDQC, LDZC, LWORK, REC
30*
31*      REAL, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
32*     $    Z( LDZ, * ), ALPHAR( * ), ALPHAI( * ), BETA( * )
33*      INTEGER, INTENT( OUT ) :: NS, ND, INFO
34*      REAL :: QC( LDQC, * ), ZC( LDZC, * ), WORK( * )
35*       ..
36*
37*
38*> \par Purpose:
39*  =============
40*>
41*> \verbatim
42*>
43*> SLAQZ3 performs AED
44*> \endverbatim
45*
46*  Arguments:
47*  ==========
48*
49*> \param[in] ILSCHUR
50*> \verbatim
51*>          ILSCHUR is LOGICAL
52*>              Determines whether or not to update the full Schur form
53*> \endverbatim
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*>          ILO and IHI mark the rows and columns of (A,B) which
81*>          are to be normalized
82*> \endverbatim
83*>
84*> \param[in] NW
85*> \verbatim
86*>          NW is INTEGER
87*>          The desired size of the deflation window.
88*> \endverbatim
89*>
90*> \param[in,out] A
91*> \verbatim
92*>          A is REAL array, dimension (LDA, N)
93*> \endverbatim
94*>
95*> \param[in] LDA
96*> \verbatim
97*>          LDA is INTEGER
98*>          The leading dimension of the array A.  LDA >= max( 1, N ).
99*> \endverbatim
100*>
101*> \param[in,out] B
102*> \verbatim
103*>          B is REAL array, dimension (LDB, N)
104*> \endverbatim
105*>
106*> \param[in] LDB
107*> \verbatim
108*>          LDB is INTEGER
109*>          The leading dimension of the array B.  LDB >= max( 1, N ).
110*> \endverbatim
111*>
112*> \param[in,out] Q
113*> \verbatim
114*>          Q is REAL array, dimension (LDQ, N)
115*> \endverbatim
116*>
117*> \param[in] LDQ
118*> \verbatim
119*>          LDQ is INTEGER
120*> \endverbatim
121*>
122*> \param[in,out] Z
123*> \verbatim
124*>          Z is REAL array, dimension (LDZ, N)
125*> \endverbatim
126*>
127*> \param[in] LDZ
128*> \verbatim
129*>          LDZ is INTEGER
130*> \endverbatim
131*>
132*> \param[out] NS
133*> \verbatim
134*>          NS is INTEGER
135*>          The number of unconverged eigenvalues available to
136*>          use as shifts.
137*> \endverbatim
138*>
139*> \param[out] ND
140*> \verbatim
141*>          ND is INTEGER
142*>          The number of converged eigenvalues found.
143*> \endverbatim
144*>
145*> \param[out] ALPHAR
146*> \verbatim
147*>          ALPHAR is REAL array, dimension (N)
148*>          The real parts of each scalar alpha defining an eigenvalue
149*>          of GNEP.
150*> \endverbatim
151*>
152*> \param[out] ALPHAI
153*> \verbatim
154*>          ALPHAI is REAL array, dimension (N)
155*>          The imaginary parts of each scalar alpha defining an
156*>          eigenvalue of GNEP.
157*>          If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
158*>          positive, then the j-th and (j+1)-st eigenvalues are a
159*>          complex conjugate pair, with ALPHAI(j+1) = -ALPHAI(j).
160*> \endverbatim
161*>
162*> \param[out] BETA
163*> \verbatim
164*>          BETA is REAL array, dimension (N)
165*>          The scalars beta that define the eigenvalues of GNEP.
166*>          Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
167*>          beta = BETA(j) represent the j-th eigenvalue of the matrix
168*>          pair (A,B), in one of the forms lambda = alpha/beta or
169*>          mu = beta/alpha.  Since either lambda or mu may overflow,
170*>          they should not, in general, be computed.
171*> \endverbatim
172*>
173*> \param[in,out] QC
174*> \verbatim
175*>          QC is REAL array, dimension (LDQC, NW)
176*> \endverbatim
177*>
178*> \param[in] LDQC
179*> \verbatim
180*>          LDQC is INTEGER
181*> \endverbatim
182*>
183*> \param[in,out] ZC
184*> \verbatim
185*>          ZC is REAL array, dimension (LDZC, NW)
186*> \endverbatim
187*>
188*> \param[in] LDZC
189*> \verbatim
190*>          LDZ is INTEGER
191*> \endverbatim
192*>
193*> \param[out] WORK
194*> \verbatim
195*>          WORK is REAL array, dimension (MAX(1,LWORK))
196*>          On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
197*> \endverbatim
198*>
199*> \param[in] LWORK
200*> \verbatim
201*>          LWORK is INTEGER
202*>          The dimension of the array WORK.  LWORK >= max(1,N).
203*>
204*>          If LWORK = -1, then a workspace query is assumed; the routine
205*>          only calculates the optimal size of the WORK array, returns
206*>          this value as the first entry of the WORK array, and no error
207*>          message related to LWORK is issued by XERBLA.
208*> \endverbatim
209*>
210*> \param[in] REC
211*> \verbatim
212*>          REC is INTEGER
213*>             REC indicates the current recursion level. Should be set
214*>             to 0 on first call.
215*> \endverbatim
216*>
217*> \param[out] INFO
218*> \verbatim
219*>          INFO is INTEGER
220*>          = 0: successful exit
221*>          < 0: if INFO = -i, the i-th argument had an illegal value
222*> \endverbatim
223*
224*  Authors:
225*  ========
226*
227*> \author Thijs Steel, KU Leuven
228*
229*> \date May 2020
230*
231*> \ingroup doubleGEcomputational
232*>
233*  =====================================================================
234      RECURSIVE SUBROUTINE SLAQZ3( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NW,
235     $                             A, LDA, B, LDB, Q, LDQ, Z, LDZ, NS,
236     $                             ND, ALPHAR, ALPHAI, BETA, QC, LDQC,
237     $                             ZC, LDZC, WORK, LWORK, REC, INFO )
238      IMPLICIT NONE
239
240*     Arguments
241      LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
242      INTEGER, INTENT( IN ) :: N, ILO, IHI, NW, LDA, LDB, LDQ, LDZ,
243     $         LDQC, LDZC, LWORK, REC
244
245      REAL, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
246     $   Z( LDZ, * ), ALPHAR( * ), ALPHAI( * ), BETA( * )
247      INTEGER, INTENT( OUT ) :: NS, ND, INFO
248      REAL :: QC( LDQC, * ), ZC( LDZC, * ), WORK( * )
249
250*     Parameters
251      REAL :: ZERO, ONE, HALF
252      PARAMETER( ZERO = 0.0, ONE = 1.0, HALF = 0.5 )
253
254*     Local Scalars
255      LOGICAL :: BULGE
256      INTEGER :: JW, KWTOP, KWBOT, ISTOPM, ISTARTM, K, K2, STGEXC_INFO,
257     $           IFST, ILST, LWORKREQ, QZ_SMALL_INFO
258      REAL :: S, SMLNUM, ULP, SAFMIN, SAFMAX, C1, S1, TEMP
259
260*     External Functions
261      EXTERNAL :: XERBLA, STGEXC, SLABAD, SLAQZ0, SLACPY, SLASET,
262     $            SLAQZ2, SROT, SLARTG, SLAG2, SGEMM
263      REAL, EXTERNAL :: SLAMCH
264
265      INFO = 0
266
267*     Set up deflation window
268      JW = MIN( NW, IHI-ILO+1 )
269      KWTOP = IHI-JW+1
270      IF ( KWTOP .EQ. ILO ) THEN
271         S = ZERO
272      ELSE
273         S = A( KWTOP, KWTOP-1 )
274      END IF
275
276*     Determine required workspace
277      IFST = 1
278      ILST = JW
279      CALL STGEXC( .TRUE., .TRUE., JW, A, LDA, B, LDB, QC, LDQC, ZC,
280     $             LDZC, IFST, ILST, WORK, -1, STGEXC_INFO )
281      LWORKREQ = INT( WORK( 1 ) )
282      CALL SLAQZ0( 'S', 'V', 'V', JW, 1, JW, A( KWTOP, KWTOP ), LDA,
283     $             B( KWTOP, KWTOP ), LDB, ALPHAR, ALPHAI, BETA, QC,
284     $             LDQC, ZC, LDZC, WORK, -1, REC+1, QZ_SMALL_INFO )
285      LWORKREQ = MAX( LWORKREQ, INT( WORK( 1 ) )+2*JW**2 )
286      LWORKREQ = MAX( LWORKREQ, N*NW, 2*NW**2+N )
287      IF ( LWORK .EQ.-1 ) THEN
288*        workspace query, quick return
289         WORK( 1 ) = LWORKREQ
290         RETURN
291      ELSE IF ( LWORK .LT. LWORKREQ ) THEN
292         INFO = -26
293      END IF
294
295      IF( INFO.NE.0 ) THEN
296         CALL XERBLA( 'SLAQZ3', -INFO )
297         RETURN
298      END IF
299
300*     Get machine constants
301      SAFMIN = SLAMCH( 'SAFE MINIMUM' )
302      SAFMAX = ONE/SAFMIN
303      CALL SLABAD( SAFMIN, SAFMAX )
304      ULP = SLAMCH( 'PRECISION' )
305      SMLNUM = SAFMIN*( REAL( N )/ULP )
306
307      IF ( IHI .EQ. KWTOP ) THEN
308*        1 by 1 deflation window, just try a regular deflation
309         ALPHAR( KWTOP ) = A( KWTOP, KWTOP )
310         ALPHAI( KWTOP ) = ZERO
311         BETA( KWTOP ) = B( KWTOP, KWTOP )
312         NS = 1
313         ND = 0
314         IF ( ABS( S ) .LE. MAX( SMLNUM, ULP*ABS( A( KWTOP,
315     $      KWTOP ) ) ) ) THEN
316            NS = 0
317            ND = 1
318            IF ( KWTOP .GT. ILO ) THEN
319               A( KWTOP, KWTOP-1 ) = ZERO
320            END IF
321         END IF
322      END IF
323
324
325*     Store window in case of convergence failure
326      CALL SLACPY( 'ALL', JW, JW, A( KWTOP, KWTOP ), LDA, WORK, JW )
327      CALL SLACPY( 'ALL', JW, JW, B( KWTOP, KWTOP ), LDB, WORK( JW**2+
328     $             1 ), JW )
329
330*     Transform window to real schur form
331      CALL SLASET( 'FULL', JW, JW, ZERO, ONE, QC, LDQC )
332      CALL SLASET( 'FULL', JW, JW, ZERO, ONE, ZC, LDZC )
333      CALL SLAQZ0( 'S', 'V', 'V', JW, 1, JW, A( KWTOP, KWTOP ), LDA,
334     $             B( KWTOP, KWTOP ), LDB, ALPHAR, ALPHAI, BETA, QC,
335     $             LDQC, ZC, LDZC, WORK( 2*JW**2+1 ), LWORK-2*JW**2,
336     $             REC+1, QZ_SMALL_INFO )
337
338      IF( QZ_SMALL_INFO .NE. 0 ) THEN
339*        Convergence failure, restore the window and exit
340         ND = 0
341         NS = JW-QZ_SMALL_INFO
342         CALL SLACPY( 'ALL', JW, JW, WORK, JW, A( KWTOP, KWTOP ), LDA )
343         CALL SLACPY( 'ALL', JW, JW, WORK( JW**2+1 ), JW, B( KWTOP,
344     $                KWTOP ), LDB )
345         RETURN
346      END IF
347
348*     Deflation detection loop
349      IF ( KWTOP .EQ. ILO .OR. S .EQ. ZERO ) THEN
350         KWBOT = KWTOP-1
351      ELSE
352         KWBOT = IHI
353         K = 1
354         K2 = 1
355         DO WHILE ( K .LE. JW )
356            BULGE = .FALSE.
357            IF ( KWBOT-KWTOP+1 .GE. 2 ) THEN
358               BULGE = A( KWBOT, KWBOT-1 ) .NE. ZERO
359            END IF
360            IF ( BULGE ) THEN
361
362*              Try to deflate complex conjugate eigenvalue pair
363               TEMP = ABS( A( KWBOT, KWBOT ) )+SQRT( ABS( A( KWBOT,
364     $            KWBOT-1 ) ) )*SQRT( ABS( A( KWBOT-1, KWBOT ) ) )
365               IF( TEMP .EQ. ZERO )THEN
366                  TEMP = ABS( S )
367               END IF
368               IF ( MAX( ABS( S*QC( 1, KWBOT-KWTOP ) ), ABS( S*QC( 1,
369     $            KWBOT-KWTOP+1 ) ) ) .LE. MAX( SMLNUM,
370     $            ULP*TEMP ) ) THEN
371*                 Deflatable
372                  KWBOT = KWBOT-2
373               ELSE
374*                 Not deflatable, move out of the way
375                  IFST = KWBOT-KWTOP+1
376                  ILST = K2
377                  CALL STGEXC( .TRUE., .TRUE., JW, A( KWTOP, KWTOP ),
378     $                         LDA, B( KWTOP, KWTOP ), LDB, QC, LDQC,
379     $                         ZC, LDZC, IFST, ILST, WORK, LWORK,
380     $                         STGEXC_INFO )
381                  K2 = K2+2
382               END IF
383               K = K+2
384            ELSE
385
386*              Try to deflate real eigenvalue
387               TEMP = ABS( A( KWBOT, KWBOT ) )
388               IF( TEMP .EQ. ZERO ) THEN
389                  TEMP = ABS( S )
390               END IF
391               IF ( ( ABS( S*QC( 1, KWBOT-KWTOP+1 ) ) ) .LE. MAX( ULP*
392     $            TEMP, SMLNUM ) ) THEN
393*                 Deflatable
394                  KWBOT = KWBOT-1
395               ELSE
396*                 Not deflatable, move out of the way
397                  IFST = KWBOT-KWTOP+1
398                  ILST = K2
399                  CALL STGEXC( .TRUE., .TRUE., JW, A( KWTOP, KWTOP ),
400     $                         LDA, B( KWTOP, KWTOP ), LDB, QC, LDQC,
401     $                         ZC, LDZC, IFST, ILST, WORK, LWORK,
402     $                         STGEXC_INFO )
403                  K2 = K2+1
404               END IF
405
406               K = K+1
407
408            END IF
409         END DO
410      END IF
411
412*     Store eigenvalues
413      ND = IHI-KWBOT
414      NS = JW-ND
415      K = KWTOP
416      DO WHILE ( K .LE. IHI )
417         BULGE = .FALSE.
418         IF ( K .LT. IHI ) THEN
419            IF ( A( K+1, K ) .NE. ZERO ) THEN
420               BULGE = .TRUE.
421            END IF
422         END IF
423         IF ( BULGE ) THEN
424*           2x2 eigenvalue block
425            CALL SLAG2( A( K, K ), LDA, B( K, K ), LDB, SAFMIN,
426     $                  BETA( K ), BETA( K+1 ), ALPHAR( K ),
427     $                  ALPHAR( K+1 ), ALPHAI( K ) )
428            ALPHAI( K+1 ) = -ALPHAI( K )
429            K = K+2
430         ELSE
431*           1x1 eigenvalue block
432            ALPHAR( K ) = A( K, K )
433            ALPHAI( K ) = ZERO
434            BETA( K ) = B( K, K )
435            K = K+1
436         END IF
437      END DO
438
439      IF ( KWTOP .NE. ILO .AND. S .NE. ZERO ) THEN
440*        Reflect spike back, this will create optimally packed bulges
441         A( KWTOP:KWBOT, KWTOP-1 ) = A( KWTOP, KWTOP-1 )*QC( 1,
442     $      1:JW-ND )
443         DO K = KWBOT-1, KWTOP, -1
444            CALL SLARTG( A( K, KWTOP-1 ), A( K+1, KWTOP-1 ), C1, S1,
445     $                   TEMP )
446            A( K, KWTOP-1 ) = TEMP
447            A( K+1, KWTOP-1 ) = ZERO
448            K2 = MAX( KWTOP, K-1 )
449            CALL SROT( IHI-K2+1, A( K, K2 ), LDA, A( K+1, K2 ), LDA, C1,
450     $                 S1 )
451            CALL SROT( IHI-( K-1 )+1, B( K, K-1 ), LDB, B( K+1, K-1 ),
452     $                 LDB, C1, S1 )
453            CALL SROT( JW, QC( 1, K-KWTOP+1 ), 1, QC( 1, K+1-KWTOP+1 ),
454     $                 1, C1, S1 )
455         END DO
456
457*        Chase bulges down
458         ISTARTM = KWTOP
459         ISTOPM = IHI
460         K = KWBOT-1
461         DO WHILE ( K .GE. KWTOP )
462            IF ( ( K .GE. KWTOP+1 ) .AND. A( K+1, K-1 ) .NE. ZERO ) THEN
463
464*              Move double pole block down and remove it
465               DO K2 = K-1, KWBOT-2
466                  CALL SLAQZ2( .TRUE., .TRUE., K2, KWTOP, KWTOP+JW-1,
467     $                         KWBOT, A, LDA, B, LDB, JW, KWTOP, QC,
468     $                         LDQC, JW, KWTOP, ZC, LDZC )
469               END DO
470
471               K = K-2
472            ELSE
473
474*              k points to single shift
475               DO K2 = K, KWBOT-2
476
477*                 Move shift down
478                  CALL SLARTG( B( K2+1, K2+1 ), B( K2+1, K2 ), C1, S1,
479     $                         TEMP )
480                  B( K2+1, K2+1 ) = TEMP
481                  B( K2+1, K2 ) = ZERO
482                  CALL SROT( K2+2-ISTARTM+1, A( ISTARTM, K2+1 ), 1,
483     $                       A( ISTARTM, K2 ), 1, C1, S1 )
484                  CALL SROT( K2-ISTARTM+1, B( ISTARTM, K2+1 ), 1,
485     $                       B( ISTARTM, K2 ), 1, C1, S1 )
486                  CALL SROT( JW, ZC( 1, K2+1-KWTOP+1 ), 1, ZC( 1,
487     $                       K2-KWTOP+1 ), 1, C1, S1 )
488
489                  CALL SLARTG( A( K2+1, K2 ), A( K2+2, K2 ), C1, S1,
490     $                         TEMP )
491                  A( K2+1, K2 ) = TEMP
492                  A( K2+2, K2 ) = ZERO
493                  CALL SROT( ISTOPM-K2, A( K2+1, K2+1 ), LDA, A( K2+2,
494     $                       K2+1 ), LDA, C1, S1 )
495                  CALL SROT( ISTOPM-K2, B( K2+1, K2+1 ), LDB, B( K2+2,
496     $                       K2+1 ), LDB, C1, S1 )
497                  CALL SROT( JW, QC( 1, K2+1-KWTOP+1 ), 1, QC( 1,
498     $                       K2+2-KWTOP+1 ), 1, C1, S1 )
499
500               END DO
501
502*              Remove the shift
503               CALL SLARTG( B( KWBOT, KWBOT ), B( KWBOT, KWBOT-1 ), C1,
504     $                      S1, TEMP )
505               B( KWBOT, KWBOT ) = TEMP
506               B( KWBOT, KWBOT-1 ) = ZERO
507               CALL SROT( KWBOT-ISTARTM, B( ISTARTM, KWBOT ), 1,
508     $                    B( ISTARTM, KWBOT-1 ), 1, C1, S1 )
509               CALL SROT( KWBOT-ISTARTM+1, A( ISTARTM, KWBOT ), 1,
510     $                    A( ISTARTM, KWBOT-1 ), 1, C1, S1 )
511               CALL SROT( JW, ZC( 1, KWBOT-KWTOP+1 ), 1, ZC( 1,
512     $                    KWBOT-1-KWTOP+1 ), 1, C1, S1 )
513
514               K = K-1
515            END IF
516         END DO
517
518      END IF
519
520*     Apply Qc and Zc to rest of the matrix
521      IF ( ILSCHUR ) THEN
522         ISTARTM = 1
523         ISTOPM = N
524      ELSE
525         ISTARTM = ILO
526         ISTOPM = IHI
527      END IF
528
529      IF ( ISTOPM-IHI > 0 ) THEN
530         CALL SGEMM( 'T', 'N', JW, ISTOPM-IHI, JW, ONE, QC, LDQC,
531     $               A( KWTOP, IHI+1 ), LDA, ZERO, WORK, JW )
532         CALL SLACPY( 'ALL', JW, ISTOPM-IHI, WORK, JW, A( KWTOP,
533     $                IHI+1 ), LDA )
534         CALL SGEMM( 'T', 'N', JW, ISTOPM-IHI, JW, ONE, QC, LDQC,
535     $               B( KWTOP, IHI+1 ), LDB, ZERO, WORK, JW )
536         CALL SLACPY( 'ALL', JW, ISTOPM-IHI, WORK, JW, B( KWTOP,
537     $                IHI+1 ), LDB )
538      END IF
539      IF ( ILQ ) THEN
540         CALL SGEMM( 'N', 'N', N, JW, JW, ONE, Q( 1, KWTOP ), LDQ, QC,
541     $               LDQC, ZERO, WORK, N )
542         CALL SLACPY( 'ALL', N, JW, WORK, N, Q( 1, KWTOP ), LDQ )
543      END IF
544
545      IF ( KWTOP-1-ISTARTM+1 > 0 ) THEN
546         CALL SGEMM( 'N', 'N', KWTOP-ISTARTM, JW, JW, ONE, A( ISTARTM,
547     $               KWTOP ), LDA, ZC, LDZC, ZERO, WORK,
548     $               KWTOP-ISTARTM )
549         CALL SLACPY( 'ALL', KWTOP-ISTARTM, JW, WORK, KWTOP-ISTARTM,
550     $                A( ISTARTM, KWTOP ), LDA )
551         CALL SGEMM( 'N', 'N', KWTOP-ISTARTM, JW, JW, ONE, B( ISTARTM,
552     $               KWTOP ), LDB, ZC, LDZC, ZERO, WORK,
553     $               KWTOP-ISTARTM )
554         CALL SLACPY( 'ALL', KWTOP-ISTARTM, JW, WORK, KWTOP-ISTARTM,
555     $                B( ISTARTM, KWTOP ), LDB )
556      END IF
557      IF ( ILZ ) THEN
558         CALL SGEMM( 'N', 'N', N, JW, JW, ONE, Z( 1, KWTOP ), LDZ, ZC,
559     $               LDZC, ZERO, WORK, N )
560         CALL SLACPY( 'ALL', N, JW, WORK, N, Z( 1, KWTOP ), LDZ )
561      END IF
562
563      END SUBROUTINE
564