1: *> \brief \b DLAQZ4
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLAQZ4 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqz4.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqz4.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqz4.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLAQZ4( 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: * DOUBLE PRECISION, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ),
32: * $ Q( LDQ, * ), Z( LDZ, * ), QC( LDQC, * ), ZC( LDZC, * ),
33: * $ WORK( * ), SR( * ), SI( * ), SS( * )
34: *
35: * INTEGER, INTENT( OUT ) :: INFO
36: * ..
37: *
38: *
39: *> \par Purpose:
40: * =============
41: *>
42: *> \verbatim
43: *>
44: *> DLAQZ4 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] SR
96: *> \verbatim
97: *> SR is DOUBLE PRECISION array. SR contains
98: *> the real parts of the shifts to use.
99: *> \endverbatim
100: *>
101: *> \param[in] SI
102: *> \verbatim
103: *> SI is DOUBLE PRECISION array. SI contains
104: *> the imaginary parts of the shifts to use.
105: *> \endverbatim
106: *>
107: *> \param[in] SS
108: *> \verbatim
109: *> SS is DOUBLE PRECISION array. SS contains
110: *> the scale of the shifts to use.
111: *> \endverbatim
112: *>
113: *> \param[in,out] A
114: *> \verbatim
115: *> A is DOUBLE PRECISION array, dimension (LDA, N)
116: *> \endverbatim
117: *>
118: *> \param[in] LDA
119: *> \verbatim
120: *> LDA is INTEGER
121: *> The leading dimension of the array A. LDA >= max( 1, N ).
122: *> \endverbatim
123: *>
124: *> \param[in,out] B
125: *> \verbatim
126: *> B is DOUBLE PRECISION array, dimension (LDB, N)
127: *> \endverbatim
128: *>
129: *> \param[in] LDB
130: *> \verbatim
131: *> LDB is INTEGER
132: *> The leading dimension of the array B. LDB >= max( 1, N ).
133: *> \endverbatim
134: *>
135: *> \param[in,out] Q
136: *> \verbatim
137: *> Q is DOUBLE PRECISION array, dimension (LDQ, N)
138: *> \endverbatim
139: *>
140: *> \param[in] LDQ
141: *> \verbatim
142: *> LDQ is INTEGER
143: *> \endverbatim
144: *>
145: *> \param[in,out] Z
146: *> \verbatim
147: *> Z is DOUBLE PRECISION array, dimension (LDZ, N)
148: *> \endverbatim
149: *>
150: *> \param[in] LDZ
151: *> \verbatim
152: *> LDZ is INTEGER
153: *> \endverbatim
154: *>
155: *> \param[in,out] QC
156: *> \verbatim
157: *> QC is DOUBLE PRECISION array, dimension (LDQC, NBLOCK_DESIRED)
158: *> \endverbatim
159: *>
160: *> \param[in] LDQC
161: *> \verbatim
162: *> LDQC is INTEGER
163: *> \endverbatim
164: *>
165: *> \param[in,out] ZC
166: *> \verbatim
167: *> ZC is DOUBLE PRECISION array, dimension (LDZC, NBLOCK_DESIRED)
168: *> \endverbatim
169: *>
170: *> \param[in] LDZC
171: *> \verbatim
172: *> LDZ is INTEGER
173: *> \endverbatim
174: *>
175: *> \param[out] WORK
176: *> \verbatim
177: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
178: *> On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
179: *> \endverbatim
180: *>
181: *> \param[in] LWORK
182: *> \verbatim
183: *> LWORK is INTEGER
184: *> The dimension of the array WORK. LWORK >= max(1,N).
185: *>
186: *> If LWORK = -1, then a workspace query is assumed; the routine
187: *> only calculates the optimal size of the WORK array, returns
188: *> this value as the first entry of the WORK array, and no error
189: *> message related to LWORK is issued by XERBLA.
190: *> \endverbatim
191: *>
192: *> \param[out] INFO
193: *> \verbatim
194: *> INFO is INTEGER
195: *> = 0: successful exit
196: *> < 0: if INFO = -i, the i-th argument had an illegal value
197: *> \endverbatim
198: *
199: * Authors:
200: * ========
201: *
202: *> \author Thijs Steel, KU Leuven
203: *
204: *> \date May 2020
205: *
206: *> \ingroup doubleGEcomputational
207: *>
208: * =====================================================================
209: SUBROUTINE DLAQZ4( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NSHIFTS,
210: $ NBLOCK_DESIRED, SR, SI, SS, A, LDA, B, LDB, Q,
211: $ LDQ, Z, LDZ, QC, LDQC, ZC, LDZC, WORK, LWORK,
212: $ INFO )
213: IMPLICIT NONE
214:
215: * Function arguments
216: LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
217: INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK,
218: $ NSHIFTS, NBLOCK_DESIRED, LDQC, LDZC
219:
220: DOUBLE PRECISION, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ),
221: $ Q( LDQ, * ), Z( LDZ, * ), QC( LDQC, * ),
222: $ ZC( LDZC, * ), WORK( * ), SR( * ), SI( * ),
223: $ SS( * )
224:
225: INTEGER, INTENT( OUT ) :: INFO
226:
227: * Parameters
228: DOUBLE PRECISION :: ZERO, ONE, HALF
229: PARAMETER( ZERO = 0.0D0, ONE = 1.0D0, HALF = 0.5D0 )
230:
231: * Local scalars
232: INTEGER :: I, J, NS, ISTARTM, ISTOPM, SHEIGHT, SWIDTH, K, NP,
233: $ ISTARTB, ISTOPB, ISHIFT, NBLOCK, NPOS
234: DOUBLE PRECISION :: TEMP, V( 3 ), C1, S1, C2, S2, SWAP
235: *
236: * External functions
237: EXTERNAL :: XERBLA, DGEMM, DLAQZ1, DLAQZ2, DLASET, DLARTG, DROT,
238: $ DLACPY
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( 'DLAQZ4', -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 DLASET( 'FULL', NS+1, NS+1, ZERO, ONE, QC, LDQC )
314: CALL DLASET( 'FULL', NS, NS, ZERO, ONE, ZC, LDZC )
315:
316: DO I = 1, NS, 2
317: * Introduce the shift
318: CALL DLAQZ1( 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 DLARTG( TEMP, V( 3 ), C1, S1, V( 2 ) )
323: CALL DLARTG( V( 1 ), V( 2 ), C2, S2, TEMP )
324:
325: CALL DROT( NS, A( ILO+1, ILO ), LDA, A( ILO+2, ILO ), LDA, C1,
326: $ S1 )
327: CALL DROT( NS, A( ILO, ILO ), LDA, A( ILO+1, ILO ), LDA, C2,
328: $ S2 )
329: CALL DROT( NS, B( ILO+1, ILO ), LDB, B( ILO+2, ILO ), LDB, C1,
330: $ S1 )
331: CALL DROT( NS, B( ILO, ILO ), LDB, B( ILO+1, ILO ), LDB, C2,
332: $ S2 )
333: CALL DROT( NS+1, QC( 1, 2 ), 1, QC( 1, 3 ), 1, C1, S1 )
334: CALL DROT( 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 DLAQZ2( .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 DGEMM( 'T', 'N', SHEIGHT, SWIDTH, SHEIGHT, ONE, QC, LDQC,
355: $ A( ILO, ILO+NS ), LDA, ZERO, WORK, SHEIGHT )
356: CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, A( ILO,
357: $ ILO+NS ), LDA )
358: CALL DGEMM( 'T', 'N', SHEIGHT, SWIDTH, SHEIGHT, ONE, QC, LDQC,
359: $ B( ILO, ILO+NS ), LDB, ZERO, WORK, SHEIGHT )
360: CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, B( ILO,
361: $ ILO+NS ), LDB )
362: END IF
363: IF ( ILQ ) THEN
364: CALL DGEMM( 'N', 'N', N, SHEIGHT, SHEIGHT, ONE, Q( 1, ILO ),
365: $ LDQ, QC, LDQC, ZERO, WORK, N )
366: CALL DLACPY( '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 DGEMM( 'N', 'N', SHEIGHT, SWIDTH, SWIDTH, ONE, A( ISTARTM,
375: $ ILO ), LDA, ZC, LDZC, ZERO, WORK, SHEIGHT )
376: CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, A( ISTARTM,
377: $ ILO ), LDA )
378: CALL DGEMM( 'N', 'N', SHEIGHT, SWIDTH, SWIDTH, ONE, B( ISTARTM,
379: $ ILO ), LDB, ZC, LDZC, ZERO, WORK, SHEIGHT )
380: CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, B( ISTARTM,
381: $ ILO ), LDB )
382: END IF
383: IF ( ILZ ) THEN
384: CALL DGEMM( 'N', 'N', N, SWIDTH, SWIDTH, ONE, Z( 1, ILO ), LDZ,
385: $ ZC, LDZC, ZERO, WORK, N )
386: CALL DLACPY( '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 DLASET( 'FULL', NS+NP, NS+NP, ZERO, ONE, QC, LDQC )
404: CALL DLASET( '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 DLAQZ2( .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 DGEMM( 'T', 'N', SHEIGHT, SWIDTH, SHEIGHT, ONE, QC,
427: $ LDQC, A( K+1, K+NS+NP ), LDA, ZERO, WORK,
428: $ SHEIGHT )
429: CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, A( K+1,
430: $ K+NS+NP ), LDA )
431: CALL DGEMM( 'T', 'N', SHEIGHT, SWIDTH, SHEIGHT, ONE, QC,
432: $ LDQC, B( K+1, K+NS+NP ), LDB, ZERO, WORK,
433: $ SHEIGHT )
434: CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, B( K+1,
435: $ K+NS+NP ), LDB )
436: END IF
437: IF ( ILQ ) THEN
438: CALL DGEMM( 'N', 'N', N, NBLOCK, NBLOCK, ONE, Q( 1, K+1 ),
439: $ LDQ, QC, LDQC, ZERO, WORK, N )
440: CALL DLACPY( '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 DGEMM( 'N', 'N', SHEIGHT, SWIDTH, SWIDTH, ONE,
449: $ A( ISTARTM, K ), LDA, ZC, LDZC, ZERO, WORK,
450: $ SHEIGHT )
451: CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT,
452: $ A( ISTARTM, K ), LDA )
453: CALL DGEMM( 'N', 'N', SHEIGHT, SWIDTH, SWIDTH, ONE,
454: $ B( ISTARTM, K ), LDB, ZC, LDZC, ZERO, WORK,
455: $ SHEIGHT )
456: CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT,
457: $ B( ISTARTM, K ), LDB )
458: END IF
459: IF ( ILZ ) THEN
460: CALL DGEMM( 'N', 'N', N, NBLOCK, NBLOCK, ONE, Z( 1, K ),
461: $ LDZ, ZC, LDZC, ZERO, WORK, N )
462: CALL DLACPY( '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 DLASET( 'FULL', NS, NS, ZERO, ONE, QC, LDQC )
473: CALL DLASET( '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 DLAQZ2( .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 DGEMM( 'T', 'N', SHEIGHT, SWIDTH, SHEIGHT, ONE, QC, LDQC,
498: $ A( IHI-NS+1, IHI+1 ), LDA, ZERO, WORK, SHEIGHT )
499: CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT,
500: $ A( IHI-NS+1, IHI+1 ), LDA )
501: CALL DGEMM( 'T', 'N', SHEIGHT, SWIDTH, SHEIGHT, ONE, QC, LDQC,
502: $ B( IHI-NS+1, IHI+1 ), LDB, ZERO, WORK, SHEIGHT )
503: CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT,
504: $ B( IHI-NS+1, IHI+1 ), LDB )
505: END IF
506: IF ( ILQ ) THEN
507: CALL DGEMM( 'N', 'N', N, NS, NS, ONE, Q( 1, IHI-NS+1 ), LDQ,
508: $ QC, LDQC, ZERO, WORK, N )
509: CALL DLACPY( '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 DGEMM( 'N', 'N', SHEIGHT, SWIDTH, SWIDTH, ONE, A( ISTARTM,
518: $ IHI-NS ), LDA, ZC, LDZC, ZERO, WORK, SHEIGHT )
519: CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, A( ISTARTM,
520: $ IHI-NS ), LDA )
521: CALL DGEMM( 'N', 'N', SHEIGHT, SWIDTH, SWIDTH, ONE, B( ISTARTM,
522: $ IHI-NS ), LDB, ZC, LDZC, ZERO, WORK, SHEIGHT )
523: CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, B( ISTARTM,
524: $ IHI-NS ), LDB )
525: END IF
526: IF ( ILZ ) THEN
527: CALL DGEMM( 'N', 'N', N, NS+1, NS+1, ONE, Z( 1, IHI-NS ), LDZ,
528: $ ZC, LDZC, ZERO, WORK, N )
529: CALL DLACPY( 'ALL', N, NS+1, WORK, N, Z( 1, IHI-NS ), LDZ )
530: END IF
531:
532: END SUBROUTINE
CVSweb interface <joel.bertrand@systella.fr>