LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ztgsyl.f
Go to the documentation of this file.
1*> \brief \b ZTGSYL
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZTGSYL + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztgsyl.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztgsyl.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztgsyl.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZTGSYL( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
22* LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK,
23* IWORK, INFO )
24*
25* .. Scalar Arguments ..
26* CHARACTER TRANS
27* INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF,
28* $ LWORK, M, N
29* DOUBLE PRECISION DIF, SCALE
30* ..
31* .. Array Arguments ..
32* INTEGER IWORK( * )
33* COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ),
34* $ D( LDD, * ), E( LDE, * ), F( LDF, * ),
35* $ WORK( * )
36* ..
37*
38*
39*> \par Purpose:
40* =============
41*>
42*> \verbatim
43*>
44*> ZTGSYL solves the generalized Sylvester equation:
45*>
46*> A * R - L * B = scale * C (1)
47*> D * R - L * E = scale * F
48*>
49*> where R and L are unknown m-by-n matrices, (A, D), (B, E) and
50*> (C, F) are given matrix pairs of size m-by-m, n-by-n and m-by-n,
51*> respectively, with complex entries. A, B, D and E are upper
52*> triangular (i.e., (A,D) and (B,E) in generalized Schur form).
53*>
54*> The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1
55*> is an output scaling factor chosen to avoid overflow.
56*>
57*> In matrix notation (1) is equivalent to solve Zx = scale*b, where Z
58*> is defined as
59*>
60*> Z = [ kron(In, A) -kron(B**H, Im) ] (2)
61*> [ kron(In, D) -kron(E**H, Im) ],
62*>
63*> Here Ix is the identity matrix of size x and X**H is the conjugate
64*> transpose of X. Kron(X, Y) is the Kronecker product between the
65*> matrices X and Y.
66*>
67*> If TRANS = 'C', y in the conjugate transposed system Z**H *y = scale*b
68*> is solved for, which is equivalent to solve for R and L in
69*>
70*> A**H * R + D**H * L = scale * C (3)
71*> R * B**H + L * E**H = scale * -F
72*>
73*> This case (TRANS = 'C') is used to compute an one-norm-based estimate
74*> of Dif[(A,D), (B,E)], the separation between the matrix pairs (A,D)
75*> and (B,E), using ZLACON.
76*>
77*> If IJOB >= 1, ZTGSYL computes a Frobenius norm-based estimate of
78*> Dif[(A,D),(B,E)]. That is, the reciprocal of a lower bound on the
79*> reciprocal of the smallest singular value of Z.
80*>
81*> This is a level-3 BLAS algorithm.
82*> \endverbatim
83*
84* Arguments:
85* ==========
86*
87*> \param[in] TRANS
88*> \verbatim
89*> TRANS is CHARACTER*1
90*> = 'N': solve the generalized sylvester equation (1).
91*> = 'C': solve the "conjugate transposed" system (3).
92*> \endverbatim
93*>
94*> \param[in] IJOB
95*> \verbatim
96*> IJOB is INTEGER
97*> Specifies what kind of functionality to be performed.
98*> =0: solve (1) only.
99*> =1: The functionality of 0 and 3.
100*> =2: The functionality of 0 and 4.
101*> =3: Only an estimate of Dif[(A,D), (B,E)] is computed.
102*> (look ahead strategy is used).
103*> =4: Only an estimate of Dif[(A,D), (B,E)] is computed.
104*> (ZGECON on sub-systems is used).
105*> Not referenced if TRANS = 'C'.
106*> \endverbatim
107*>
108*> \param[in] M
109*> \verbatim
110*> M is INTEGER
111*> The order of the matrices A and D, and the row dimension of
112*> the matrices C, F, R and L.
113*> \endverbatim
114*>
115*> \param[in] N
116*> \verbatim
117*> N is INTEGER
118*> The order of the matrices B and E, and the column dimension
119*> of the matrices C, F, R and L.
120*> \endverbatim
121*>
122*> \param[in] A
123*> \verbatim
124*> A is COMPLEX*16 array, dimension (LDA, M)
125*> The upper triangular matrix A.
126*> \endverbatim
127*>
128*> \param[in] LDA
129*> \verbatim
130*> LDA is INTEGER
131*> The leading dimension of the array A. LDA >= max(1, M).
132*> \endverbatim
133*>
134*> \param[in] B
135*> \verbatim
136*> B is COMPLEX*16 array, dimension (LDB, N)
137*> The upper triangular matrix B.
138*> \endverbatim
139*>
140*> \param[in] LDB
141*> \verbatim
142*> LDB is INTEGER
143*> The leading dimension of the array B. LDB >= max(1, N).
144*> \endverbatim
145*>
146*> \param[in,out] C
147*> \verbatim
148*> C is COMPLEX*16 array, dimension (LDC, N)
149*> On entry, C contains the right-hand-side of the first matrix
150*> equation in (1) or (3).
151*> On exit, if IJOB = 0, 1 or 2, C has been overwritten by
152*> the solution R. If IJOB = 3 or 4 and TRANS = 'N', C holds R,
153*> the solution achieved during the computation of the
154*> Dif-estimate.
155*> \endverbatim
156*>
157*> \param[in] LDC
158*> \verbatim
159*> LDC is INTEGER
160*> The leading dimension of the array C. LDC >= max(1, M).
161*> \endverbatim
162*>
163*> \param[in] D
164*> \verbatim
165*> D is COMPLEX*16 array, dimension (LDD, M)
166*> The upper triangular matrix D.
167*> \endverbatim
168*>
169*> \param[in] LDD
170*> \verbatim
171*> LDD is INTEGER
172*> The leading dimension of the array D. LDD >= max(1, M).
173*> \endverbatim
174*>
175*> \param[in] E
176*> \verbatim
177*> E is COMPLEX*16 array, dimension (LDE, N)
178*> The upper triangular matrix E.
179*> \endverbatim
180*>
181*> \param[in] LDE
182*> \verbatim
183*> LDE is INTEGER
184*> The leading dimension of the array E. LDE >= max(1, N).
185*> \endverbatim
186*>
187*> \param[in,out] F
188*> \verbatim
189*> F is COMPLEX*16 array, dimension (LDF, N)
190*> On entry, F contains the right-hand-side of the second matrix
191*> equation in (1) or (3).
192*> On exit, if IJOB = 0, 1 or 2, F has been overwritten by
193*> the solution L. If IJOB = 3 or 4 and TRANS = 'N', F holds L,
194*> the solution achieved during the computation of the
195*> Dif-estimate.
196*> \endverbatim
197*>
198*> \param[in] LDF
199*> \verbatim
200*> LDF is INTEGER
201*> The leading dimension of the array F. LDF >= max(1, M).
202*> \endverbatim
203*>
204*> \param[out] DIF
205*> \verbatim
206*> DIF is DOUBLE PRECISION
207*> On exit DIF is the reciprocal of a lower bound of the
208*> reciprocal of the Dif-function, i.e. DIF is an upper bound of
209*> Dif[(A,D), (B,E)] = sigma-min(Z), where Z as in (2).
210*> IF IJOB = 0 or TRANS = 'C', DIF is not referenced.
211*> \endverbatim
212*>
213*> \param[out] SCALE
214*> \verbatim
215*> SCALE is DOUBLE PRECISION
216*> On exit SCALE is the scaling factor in (1) or (3).
217*> If 0 < SCALE < 1, C and F hold the solutions R and L, resp.,
218*> to a slightly perturbed system but the input matrices A, B,
219*> D and E have not been changed. If SCALE = 0, R and L will
220*> hold the solutions to the homogeneous system with C = F = 0.
221*> \endverbatim
222*>
223*> \param[out] WORK
224*> \verbatim
225*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
226*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
227*> \endverbatim
228*>
229*> \param[in] LWORK
230*> \verbatim
231*> LWORK is INTEGER
232*> The dimension of the array WORK. LWORK > = 1.
233*> If IJOB = 1 or 2 and TRANS = 'N', LWORK >= max(1,2*M*N).
234*>
235*> If LWORK = -1, then a workspace query is assumed; the routine
236*> only calculates the optimal size of the WORK array, returns
237*> this value as the first entry of the WORK array, and no error
238*> message related to LWORK is issued by XERBLA.
239*> \endverbatim
240*>
241*> \param[out] IWORK
242*> \verbatim
243*> IWORK is INTEGER array, dimension (M+N+2)
244*> \endverbatim
245*>
246*> \param[out] INFO
247*> \verbatim
248*> INFO is INTEGER
249*> =0: successful exit
250*> <0: If INFO = -i, the i-th argument had an illegal value.
251*> >0: (A, D) and (B, E) have common or very close
252*> eigenvalues.
253*> \endverbatim
254*
255* Authors:
256* ========
257*
258*> \author Univ. of Tennessee
259*> \author Univ. of California Berkeley
260*> \author Univ. of Colorado Denver
261*> \author NAG Ltd.
262*
263*> \ingroup tgsyl
264*
265*> \par Contributors:
266* ==================
267*>
268*> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
269*> Umea University, S-901 87 Umea, Sweden.
270*
271*> \par References:
272* ================
273*>
274*> [1] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
275*> for Solving the Generalized Sylvester Equation and Estimating the
276*> Separation between Regular Matrix Pairs, Report UMINF - 93.23,
277*> Department of Computing Science, Umea University, S-901 87 Umea,
278*> Sweden, December 1993, Revised April 1994, Also as LAPACK Working
279*> Note 75. To appear in ACM Trans. on Math. Software, Vol 22,
280*> No 1, 1996.
281*> \n
282*> [2] B. Kagstrom, A Perturbation Analysis of the Generalized Sylvester
283*> Equation (AR - LB, DR - LE ) = (C, F), SIAM J. Matrix Anal.
284*> Appl., 15(4):1045-1060, 1994.
285*> \n
286*> [3] B. Kagstrom and L. Westin, Generalized Schur Methods with
287*> Condition Estimators for Solving the Generalized Sylvester
288*> Equation, IEEE Transactions on Automatic Control, Vol. 34, No. 7,
289*> July 1989, pp 745-751.
290*>
291* =====================================================================
292 SUBROUTINE ztgsyl( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
293 $ LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK,
294 $ IWORK, INFO )
295*
296* -- LAPACK computational routine --
297* -- LAPACK is a software package provided by Univ. of Tennessee, --
298* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
299*
300* .. Scalar Arguments ..
301 CHARACTER TRANS
302 INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF,
303 $ lwork, m, n
304 DOUBLE PRECISION DIF, SCALE
305* ..
306* .. Array Arguments ..
307 INTEGER IWORK( * )
308 COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ),
309 $ d( ldd, * ), e( lde, * ), f( ldf, * ),
310 $ work( * )
311* ..
312*
313* =====================================================================
314* Replaced various illegal calls to CCOPY by calls to CLASET.
315* Sven Hammarling, 1/5/02.
316*
317* .. Parameters ..
318 DOUBLE PRECISION ZERO, ONE
319 PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0 )
320 COMPLEX*16 CZERO
321 parameter( czero = (0.0d+0, 0.0d+0) )
322* ..
323* .. Local Scalars ..
324 LOGICAL LQUERY, NOTRAN
325 INTEGER I, IE, IFUNC, IROUND, IS, ISOLVE, J, JE, JS, K,
326 $ linfo, lwmin, mb, nb, p, pq, q
327 DOUBLE PRECISION DSCALE, DSUM, SCALE2, SCALOC
328* ..
329* .. External Functions ..
330 LOGICAL LSAME
331 INTEGER ILAENV
332 EXTERNAL lsame, ilaenv
333* ..
334* .. External Subroutines ..
335 EXTERNAL xerbla, zgemm, zlacpy, zlaset, zscal, ztgsy2
336* ..
337* .. Intrinsic Functions ..
338 INTRINSIC dble, dcmplx, max, sqrt
339* ..
340* .. Executable Statements ..
341*
342* Decode and test input parameters
343*
344 info = 0
345 notran = lsame( trans, 'N' )
346 lquery = ( lwork.EQ.-1 )
347*
348 IF( .NOT.notran .AND. .NOT.lsame( trans, 'C' ) ) THEN
349 info = -1
350 ELSE IF( notran ) THEN
351 IF( ( ijob.LT.0 ) .OR. ( ijob.GT.4 ) ) THEN
352 info = -2
353 END IF
354 END IF
355 IF( info.EQ.0 ) THEN
356 IF( m.LE.0 ) THEN
357 info = -3
358 ELSE IF( n.LE.0 ) THEN
359 info = -4
360 ELSE IF( lda.LT.max( 1, m ) ) THEN
361 info = -6
362 ELSE IF( ldb.LT.max( 1, n ) ) THEN
363 info = -8
364 ELSE IF( ldc.LT.max( 1, m ) ) THEN
365 info = -10
366 ELSE IF( ldd.LT.max( 1, m ) ) THEN
367 info = -12
368 ELSE IF( lde.LT.max( 1, n ) ) THEN
369 info = -14
370 ELSE IF( ldf.LT.max( 1, m ) ) THEN
371 info = -16
372 END IF
373 END IF
374*
375 IF( info.EQ.0 ) THEN
376 IF( notran ) THEN
377 IF( ijob.EQ.1 .OR. ijob.EQ.2 ) THEN
378 lwmin = max( 1, 2*m*n )
379 ELSE
380 lwmin = 1
381 END IF
382 ELSE
383 lwmin = 1
384 END IF
385 work( 1 ) = lwmin
386*
387 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
388 info = -20
389 END IF
390 END IF
391*
392 IF( info.NE.0 ) THEN
393 CALL xerbla( 'ZTGSYL', -info )
394 RETURN
395 ELSE IF( lquery ) THEN
396 RETURN
397 END IF
398*
399* Quick return if possible
400*
401 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
402 scale = 1
403 IF( notran ) THEN
404 IF( ijob.NE.0 ) THEN
405 dif = 0
406 END IF
407 END IF
408 RETURN
409 END IF
410*
411* Determine optimal block sizes MB and NB
412*
413 mb = ilaenv( 2, 'ZTGSYL', trans, m, n, -1, -1 )
414 nb = ilaenv( 5, 'ZTGSYL', trans, m, n, -1, -1 )
415*
416 isolve = 1
417 ifunc = 0
418 IF( notran ) THEN
419 IF( ijob.GE.3 ) THEN
420 ifunc = ijob - 2
421 CALL zlaset( 'F', m, n, czero, czero, c, ldc )
422 CALL zlaset( 'F', m, n, czero, czero, f, ldf )
423 ELSE IF( ijob.GE.1 .AND. notran ) THEN
424 isolve = 2
425 END IF
426 END IF
427*
428 IF( ( mb.LE.1 .AND. nb.LE.1 ) .OR. ( mb.GE.m .AND. nb.GE.n ) )
429 $ THEN
430*
431* Use unblocked Level 2 solver
432*
433 DO 30 iround = 1, isolve
434*
435 scale = one
436 dscale = zero
437 dsum = one
438 pq = m*n
439 CALL ztgsy2( trans, ifunc, m, n, a, lda, b, ldb, c, ldc, d,
440 $ ldd, e, lde, f, ldf, scale, dsum, dscale,
441 $ info )
442 IF( dscale.NE.zero ) THEN
443 IF( ijob.EQ.1 .OR. ijob.EQ.3 ) THEN
444 dif = sqrt( dble( 2*m*n ) ) / ( dscale*sqrt( dsum ) )
445 ELSE
446 dif = sqrt( dble( pq ) ) / ( dscale*sqrt( dsum ) )
447 END IF
448 END IF
449 IF( isolve.EQ.2 .AND. iround.EQ.1 ) THEN
450 IF( notran ) THEN
451 ifunc = ijob
452 END IF
453 scale2 = scale
454 CALL zlacpy( 'F', m, n, c, ldc, work, m )
455 CALL zlacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
456 CALL zlaset( 'F', m, n, czero, czero, c, ldc )
457 CALL zlaset( 'F', m, n, czero, czero, f, ldf )
458 ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
459 CALL zlacpy( 'F', m, n, work, m, c, ldc )
460 CALL zlacpy( 'F', m, n, work( m*n+1 ), m, f, ldf )
461 scale = scale2
462 END IF
463 30 CONTINUE
464*
465 RETURN
466*
467 END IF
468*
469* Determine block structure of A
470*
471 p = 0
472 i = 1
473 40 CONTINUE
474 IF( i.GT.m )
475 $ GO TO 50
476 p = p + 1
477 iwork( p ) = i
478 i = i + mb
479 IF( i.GE.m )
480 $ GO TO 50
481 GO TO 40
482 50 CONTINUE
483 iwork( p+1 ) = m + 1
484 IF( iwork( p ).EQ.iwork( p+1 ) )
485 $ p = p - 1
486*
487* Determine block structure of B
488*
489 q = p + 1
490 j = 1
491 60 CONTINUE
492 IF( j.GT.n )
493 $ GO TO 70
494*
495 q = q + 1
496 iwork( q ) = j
497 j = j + nb
498 IF( j.GE.n )
499 $ GO TO 70
500 GO TO 60
501*
502 70 CONTINUE
503 iwork( q+1 ) = n + 1
504 IF( iwork( q ).EQ.iwork( q+1 ) )
505 $ q = q - 1
506*
507 IF( notran ) THEN
508 DO 150 iround = 1, isolve
509*
510* Solve (I, J) - subsystem
511* A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
512* D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
513* for I = P, P - 1, ..., 1; J = 1, 2, ..., Q
514*
515 pq = 0
516 scale = one
517 dscale = zero
518 dsum = one
519 DO 130 j = p + 2, q
520 js = iwork( j )
521 je = iwork( j+1 ) - 1
522 nb = je - js + 1
523 DO 120 i = p, 1, -1
524 is = iwork( i )
525 ie = iwork( i+1 ) - 1
526 mb = ie - is + 1
527 CALL ztgsy2( trans, ifunc, mb, nb, a( is, is ), lda,
528 $ b( js, js ), ldb, c( is, js ), ldc,
529 $ d( is, is ), ldd, e( js, js ), lde,
530 $ f( is, js ), ldf, scaloc, dsum, dscale,
531 $ linfo )
532 IF( linfo.GT.0 )
533 $ info = linfo
534 pq = pq + mb*nb
535 IF( scaloc.NE.one ) THEN
536 DO 80 k = 1, js - 1
537 CALL zscal( m, dcmplx( scaloc, zero ),
538 $ c( 1, k ), 1 )
539 CALL zscal( m, dcmplx( scaloc, zero ),
540 $ f( 1, k ), 1 )
541 80 CONTINUE
542 DO 90 k = js, je
543 CALL zscal( is-1, dcmplx( scaloc, zero ),
544 $ c( 1, k ), 1 )
545 CALL zscal( is-1, dcmplx( scaloc, zero ),
546 $ f( 1, k ), 1 )
547 90 CONTINUE
548 DO 100 k = js, je
549 CALL zscal( m-ie, dcmplx( scaloc, zero ),
550 $ c( ie+1, k ), 1 )
551 CALL zscal( m-ie, dcmplx( scaloc, zero ),
552 $ f( ie+1, k ), 1 )
553 100 CONTINUE
554 DO 110 k = je + 1, n
555 CALL zscal( m, dcmplx( scaloc, zero ),
556 $ c( 1, k ), 1 )
557 CALL zscal( m, dcmplx( scaloc, zero ),
558 $ f( 1, k ), 1 )
559 110 CONTINUE
560 scale = scale*scaloc
561 END IF
562*
563* Substitute R(I,J) and L(I,J) into remaining equation.
564*
565 IF( i.GT.1 ) THEN
566 CALL zgemm( 'N', 'N', is-1, nb, mb,
567 $ dcmplx( -one, zero ), a( 1, is ), lda,
568 $ c( is, js ), ldc, dcmplx( one, zero ),
569 $ c( 1, js ), ldc )
570 CALL zgemm( 'N', 'N', is-1, nb, mb,
571 $ dcmplx( -one, zero ), d( 1, is ), ldd,
572 $ c( is, js ), ldc, dcmplx( one, zero ),
573 $ f( 1, js ), ldf )
574 END IF
575 IF( j.LT.q ) THEN
576 CALL zgemm( 'N', 'N', mb, n-je, nb,
577 $ dcmplx( one, zero ), f( is, js ), ldf,
578 $ b( js, je+1 ), ldb,
579 $ dcmplx( one, zero ), c( is, je+1 ),
580 $ ldc )
581 CALL zgemm( 'N', 'N', mb, n-je, nb,
582 $ dcmplx( one, zero ), f( is, js ), ldf,
583 $ e( js, je+1 ), lde,
584 $ dcmplx( one, zero ), f( is, je+1 ),
585 $ ldf )
586 END IF
587 120 CONTINUE
588 130 CONTINUE
589 IF( dscale.NE.zero ) THEN
590 IF( ijob.EQ.1 .OR. ijob.EQ.3 ) THEN
591 dif = sqrt( dble( 2*m*n ) ) / ( dscale*sqrt( dsum ) )
592 ELSE
593 dif = sqrt( dble( pq ) ) / ( dscale*sqrt( dsum ) )
594 END IF
595 END IF
596 IF( isolve.EQ.2 .AND. iround.EQ.1 ) THEN
597 IF( notran ) THEN
598 ifunc = ijob
599 END IF
600 scale2 = scale
601 CALL zlacpy( 'F', m, n, c, ldc, work, m )
602 CALL zlacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
603 CALL zlaset( 'F', m, n, czero, czero, c, ldc )
604 CALL zlaset( 'F', m, n, czero, czero, f, ldf )
605 ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
606 CALL zlacpy( 'F', m, n, work, m, c, ldc )
607 CALL zlacpy( 'F', m, n, work( m*n+1 ), m, f, ldf )
608 scale = scale2
609 END IF
610 150 CONTINUE
611 ELSE
612*
613* Solve transposed (I, J)-subsystem
614* A(I, I)**H * R(I, J) + D(I, I)**H * L(I, J) = C(I, J)
615* R(I, J) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
616* for I = 1,2,..., P; J = Q, Q-1,..., 1
617*
618 scale = one
619 DO 210 i = 1, p
620 is = iwork( i )
621 ie = iwork( i+1 ) - 1
622 mb = ie - is + 1
623 DO 200 j = q, p + 2, -1
624 js = iwork( j )
625 je = iwork( j+1 ) - 1
626 nb = je - js + 1
627 CALL ztgsy2( trans, ifunc, mb, nb, a( is, is ), lda,
628 $ b( js, js ), ldb, c( is, js ), ldc,
629 $ d( is, is ), ldd, e( js, js ), lde,
630 $ f( is, js ), ldf, scaloc, dsum, dscale,
631 $ linfo )
632 IF( linfo.GT.0 )
633 $ info = linfo
634 IF( scaloc.NE.one ) THEN
635 DO 160 k = 1, js - 1
636 CALL zscal( m, dcmplx( scaloc, zero ), c( 1, k ),
637 $ 1 )
638 CALL zscal( m, dcmplx( scaloc, zero ), f( 1, k ),
639 $ 1 )
640 160 CONTINUE
641 DO 170 k = js, je
642 CALL zscal( is-1, dcmplx( scaloc, zero ),
643 $ c( 1, k ), 1 )
644 CALL zscal( is-1, dcmplx( scaloc, zero ),
645 $ f( 1, k ), 1 )
646 170 CONTINUE
647 DO 180 k = js, je
648 CALL zscal( m-ie, dcmplx( scaloc, zero ),
649 $ c( ie+1, k ), 1 )
650 CALL zscal( m-ie, dcmplx( scaloc, zero ),
651 $ f( ie+1, k ), 1 )
652 180 CONTINUE
653 DO 190 k = je + 1, n
654 CALL zscal( m, dcmplx( scaloc, zero ), c( 1, k ),
655 $ 1 )
656 CALL zscal( m, dcmplx( scaloc, zero ), f( 1, k ),
657 $ 1 )
658 190 CONTINUE
659 scale = scale*scaloc
660 END IF
661*
662* Substitute R(I,J) and L(I,J) into remaining equation.
663*
664 IF( j.GT.p+2 ) THEN
665 CALL zgemm( 'N', 'C', mb, js-1, nb,
666 $ dcmplx( one, zero ), c( is, js ), ldc,
667 $ b( 1, js ), ldb, dcmplx( one, zero ),
668 $ f( is, 1 ), ldf )
669 CALL zgemm( 'N', 'C', mb, js-1, nb,
670 $ dcmplx( one, zero ), f( is, js ), ldf,
671 $ e( 1, js ), lde, dcmplx( one, zero ),
672 $ f( is, 1 ), ldf )
673 END IF
674 IF( i.LT.p ) THEN
675 CALL zgemm( 'C', 'N', m-ie, nb, mb,
676 $ dcmplx( -one, zero ), a( is, ie+1 ), lda,
677 $ c( is, js ), ldc, dcmplx( one, zero ),
678 $ c( ie+1, js ), ldc )
679 CALL zgemm( 'C', 'N', m-ie, nb, mb,
680 $ dcmplx( -one, zero ), d( is, ie+1 ), ldd,
681 $ f( is, js ), ldf, dcmplx( one, zero ),
682 $ c( ie+1, js ), ldc )
683 END IF
684 200 CONTINUE
685 210 CONTINUE
686 END IF
687*
688 work( 1 ) = lwmin
689*
690 RETURN
691*
692* End of ZTGSYL
693*
694 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:103
subroutine zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition zlaset.f:106
subroutine zscal(n, za, zx, incx)
ZSCAL
Definition zscal.f:78
subroutine ztgsy2(trans, ijob, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, rdsum, rdscal, info)
ZTGSY2 solves the generalized Sylvester equation (unblocked algorithm).
Definition ztgsy2.f:259
subroutine ztgsyl(trans, ijob, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, dif, work, lwork, iwork, info)
ZTGSYL
Definition ztgsyl.f:295