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