LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
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 homogenious 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 *> \date December 2016
264 *
265 *> \ingroup complex16SYcomputational
266 *
267 *> \par Contributors:
268 * ==================
269 *>
270 *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
271 *> Umea University, S-901 87 Umea, Sweden.
272 *
273 *> \par References:
274 * ================
275 *>
276 *> [1] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
277 *> for Solving the Generalized Sylvester Equation and Estimating the
278 *> Separation between Regular Matrix Pairs, Report UMINF - 93.23,
279 *> Department of Computing Science, Umea University, S-901 87 Umea,
280 *> Sweden, December 1993, Revised April 1994, Also as LAPACK Working
281 *> Note 75. To appear in ACM Trans. on Math. Software, Vol 22,
282 *> No 1, 1996.
283 *> \n
284 *> [2] B. Kagstrom, A Perturbation Analysis of the Generalized Sylvester
285 *> Equation (AR - LB, DR - LE ) = (C, F), SIAM J. Matrix Anal.
286 *> Appl., 15(4):1045-1060, 1994.
287 *> \n
288 *> [3] B. Kagstrom and L. Westin, Generalized Schur Methods with
289 *> Condition Estimators for Solving the Generalized Sylvester
290 *> Equation, IEEE Transactions on Automatic Control, Vol. 34, No. 7,
291 *> July 1989, pp 745-751.
292 *>
293 * =====================================================================
294  SUBROUTINE ztgsyl( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
295  $ LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK,
296  $ IWORK, INFO )
297 *
298 * -- LAPACK computational routine (version 3.7.0) --
299 * -- LAPACK is a software package provided by Univ. of Tennessee, --
300 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
301 * December 2016
302 *
303 * .. Scalar Arguments ..
304  CHARACTER TRANS
305  INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF,
306  $ lwork, m, n
307  DOUBLE PRECISION DIF, SCALE
308 * ..
309 * .. Array Arguments ..
310  INTEGER IWORK( * )
311  COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ),
312  $ d( ldd, * ), e( lde, * ), f( ldf, * ),
313  $ work( * )
314 * ..
315 *
316 * =====================================================================
317 * Replaced various illegal calls to CCOPY by calls to CLASET.
318 * Sven Hammarling, 1/5/02.
319 *
320 * .. Parameters ..
321  DOUBLE PRECISION ZERO, ONE
322  PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0 )
323  COMPLEX*16 CZERO
324  parameter( czero = (0.0d+0, 0.0d+0) )
325 * ..
326 * .. Local Scalars ..
327  LOGICAL LQUERY, NOTRAN
328  INTEGER I, IE, IFUNC, IROUND, IS, ISOLVE, J, JE, JS, K,
329  $ linfo, lwmin, mb, nb, p, pq, q
330  DOUBLE PRECISION DSCALE, DSUM, SCALE2, SCALOC
331 * ..
332 * .. External Functions ..
333  LOGICAL LSAME
334  INTEGER ILAENV
335  EXTERNAL lsame, ilaenv
336 * ..
337 * .. External Subroutines ..
338  EXTERNAL xerbla, zgemm, zlacpy, zlaset, zscal, ztgsy2
339 * ..
340 * .. Intrinsic Functions ..
341  INTRINSIC dble, dcmplx, max, sqrt
342 * ..
343 * .. Executable Statements ..
344 *
345 * Decode and test input parameters
346 *
347  info = 0
348  notran = lsame( trans, 'N' )
349  lquery = ( lwork.EQ.-1 )
350 *
351  IF( .NOT.notran .AND. .NOT.lsame( trans, 'C' ) ) THEN
352  info = -1
353  ELSE IF( notran ) THEN
354  IF( ( ijob.LT.0 ) .OR. ( ijob.GT.4 ) ) THEN
355  info = -2
356  END IF
357  END IF
358  IF( info.EQ.0 ) THEN
359  IF( m.LE.0 ) THEN
360  info = -3
361  ELSE IF( n.LE.0 ) THEN
362  info = -4
363  ELSE IF( lda.LT.max( 1, m ) ) THEN
364  info = -6
365  ELSE IF( ldb.LT.max( 1, n ) ) THEN
366  info = -8
367  ELSE IF( ldc.LT.max( 1, m ) ) THEN
368  info = -10
369  ELSE IF( ldd.LT.max( 1, m ) ) THEN
370  info = -12
371  ELSE IF( lde.LT.max( 1, n ) ) THEN
372  info = -14
373  ELSE IF( ldf.LT.max( 1, m ) ) THEN
374  info = -16
375  END IF
376  END IF
377 *
378  IF( info.EQ.0 ) THEN
379  IF( notran ) THEN
380  IF( ijob.EQ.1 .OR. ijob.EQ.2 ) THEN
381  lwmin = max( 1, 2*m*n )
382  ELSE
383  lwmin = 1
384  END IF
385  ELSE
386  lwmin = 1
387  END IF
388  work( 1 ) = lwmin
389 *
390  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
391  info = -20
392  END IF
393  END IF
394 *
395  IF( info.NE.0 ) THEN
396  CALL xerbla( 'ZTGSYL', -info )
397  RETURN
398  ELSE IF( lquery ) THEN
399  RETURN
400  END IF
401 *
402 * Quick return if possible
403 *
404  IF( m.EQ.0 .OR. n.EQ.0 ) THEN
405  scale = 1
406  IF( notran ) THEN
407  IF( ijob.NE.0 ) THEN
408  dif = 0
409  END IF
410  END IF
411  RETURN
412  END IF
413 *
414 * Determine optimal block sizes MB and NB
415 *
416  mb = ilaenv( 2, 'ZTGSYL', trans, m, n, -1, -1 )
417  nb = ilaenv( 5, 'ZTGSYL', trans, m, n, -1, -1 )
418 *
419  isolve = 1
420  ifunc = 0
421  IF( notran ) THEN
422  IF( ijob.GE.3 ) THEN
423  ifunc = ijob - 2
424  CALL zlaset( 'F', m, n, czero, czero, c, ldc )
425  CALL zlaset( 'F', m, n, czero, czero, f, ldf )
426  ELSE IF( ijob.GE.1 .AND. notran ) THEN
427  isolve = 2
428  END IF
429  END IF
430 *
431  IF( ( mb.LE.1 .AND. nb.LE.1 ) .OR. ( mb.GE.m .AND. nb.GE.n ) )
432  $ THEN
433 *
434 * Use unblocked Level 2 solver
435 *
436  DO 30 iround = 1, isolve
437 *
438  scale = one
439  dscale = zero
440  dsum = one
441  pq = m*n
442  CALL ztgsy2( trans, ifunc, m, n, a, lda, b, ldb, c, ldc, d,
443  $ ldd, e, lde, f, ldf, scale, dsum, dscale,
444  $ info )
445  IF( dscale.NE.zero ) THEN
446  IF( ijob.EQ.1 .OR. ijob.EQ.3 ) THEN
447  dif = sqrt( dble( 2*m*n ) ) / ( dscale*sqrt( dsum ) )
448  ELSE
449  dif = sqrt( dble( pq ) ) / ( dscale*sqrt( dsum ) )
450  END IF
451  END IF
452  IF( isolve.EQ.2 .AND. iround.EQ.1 ) THEN
453  IF( notran ) THEN
454  ifunc = ijob
455  END IF
456  scale2 = scale
457  CALL zlacpy( 'F', m, n, c, ldc, work, m )
458  CALL zlacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
459  CALL zlaset( 'F', m, n, czero, czero, c, ldc )
460  CALL zlaset( 'F', m, n, czero, czero, f, ldf )
461  ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
462  CALL zlacpy( 'F', m, n, work, m, c, ldc )
463  CALL zlacpy( 'F', m, n, work( m*n+1 ), m, f, ldf )
464  scale = scale2
465  END IF
466  30 CONTINUE
467 *
468  RETURN
469 *
470  END IF
471 *
472 * Determine block structure of A
473 *
474  p = 0
475  i = 1
476  40 CONTINUE
477  IF( i.GT.m )
478  $ GO TO 50
479  p = p + 1
480  iwork( p ) = i
481  i = i + mb
482  IF( i.GE.m )
483  $ GO TO 50
484  GO TO 40
485  50 CONTINUE
486  iwork( p+1 ) = m + 1
487  IF( iwork( p ).EQ.iwork( p+1 ) )
488  $ p = p - 1
489 *
490 * Determine block structure of B
491 *
492  q = p + 1
493  j = 1
494  60 CONTINUE
495  IF( j.GT.n )
496  $ GO TO 70
497 *
498  q = q + 1
499  iwork( q ) = j
500  j = j + nb
501  IF( j.GE.n )
502  $ GO TO 70
503  GO TO 60
504 *
505  70 CONTINUE
506  iwork( q+1 ) = n + 1
507  IF( iwork( q ).EQ.iwork( q+1 ) )
508  $ q = q - 1
509 *
510  IF( notran ) THEN
511  DO 150 iround = 1, isolve
512 *
513 * Solve (I, J) - subsystem
514 * A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
515 * D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
516 * for I = P, P - 1, ..., 1; J = 1, 2, ..., Q
517 *
518  pq = 0
519  scale = one
520  dscale = zero
521  dsum = one
522  DO 130 j = p + 2, q
523  js = iwork( j )
524  je = iwork( j+1 ) - 1
525  nb = je - js + 1
526  DO 120 i = p, 1, -1
527  is = iwork( i )
528  ie = iwork( i+1 ) - 1
529  mb = ie - is + 1
530  CALL ztgsy2( trans, ifunc, mb, nb, a( is, is ), lda,
531  $ b( js, js ), ldb, c( is, js ), ldc,
532  $ d( is, is ), ldd, e( js, js ), lde,
533  $ f( is, js ), ldf, scaloc, dsum, dscale,
534  $ linfo )
535  IF( linfo.GT.0 )
536  $ info = linfo
537  pq = pq + mb*nb
538  IF( scaloc.NE.one ) THEN
539  DO 80 k = 1, js - 1
540  CALL zscal( m, dcmplx( scaloc, zero ),
541  $ c( 1, k ), 1 )
542  CALL zscal( m, dcmplx( scaloc, zero ),
543  $ f( 1, k ), 1 )
544  80 CONTINUE
545  DO 90 k = js, je
546  CALL zscal( is-1, dcmplx( scaloc, zero ),
547  $ c( 1, k ), 1 )
548  CALL zscal( is-1, dcmplx( scaloc, zero ),
549  $ f( 1, k ), 1 )
550  90 CONTINUE
551  DO 100 k = js, je
552  CALL zscal( m-ie, dcmplx( scaloc, zero ),
553  $ c( ie+1, k ), 1 )
554  CALL zscal( m-ie, dcmplx( scaloc, zero ),
555  $ f( ie+1, k ), 1 )
556  100 CONTINUE
557  DO 110 k = je + 1, n
558  CALL zscal( m, dcmplx( scaloc, zero ),
559  $ c( 1, k ), 1 )
560  CALL zscal( m, dcmplx( scaloc, zero ),
561  $ f( 1, k ), 1 )
562  110 CONTINUE
563  scale = scale*scaloc
564  END IF
565 *
566 * Substitute R(I,J) and L(I,J) into remaining equation.
567 *
568  IF( i.GT.1 ) THEN
569  CALL zgemm( 'N', 'N', is-1, nb, mb,
570  $ dcmplx( -one, zero ), a( 1, is ), lda,
571  $ c( is, js ), ldc, dcmplx( one, zero ),
572  $ c( 1, js ), ldc )
573  CALL zgemm( 'N', 'N', is-1, nb, mb,
574  $ dcmplx( -one, zero ), d( 1, is ), ldd,
575  $ c( is, js ), ldc, dcmplx( one, zero ),
576  $ f( 1, js ), ldf )
577  END IF
578  IF( j.LT.q ) THEN
579  CALL zgemm( 'N', 'N', mb, n-je, nb,
580  $ dcmplx( one, zero ), f( is, js ), ldf,
581  $ b( js, je+1 ), ldb,
582  $ dcmplx( one, zero ), c( is, je+1 ),
583  $ ldc )
584  CALL zgemm( 'N', 'N', mb, n-je, nb,
585  $ dcmplx( one, zero ), f( is, js ), ldf,
586  $ e( js, je+1 ), lde,
587  $ dcmplx( one, zero ), f( is, je+1 ),
588  $ ldf )
589  END IF
590  120 CONTINUE
591  130 CONTINUE
592  IF( dscale.NE.zero ) THEN
593  IF( ijob.EQ.1 .OR. ijob.EQ.3 ) THEN
594  dif = sqrt( dble( 2*m*n ) ) / ( dscale*sqrt( dsum ) )
595  ELSE
596  dif = sqrt( dble( pq ) ) / ( dscale*sqrt( dsum ) )
597  END IF
598  END IF
599  IF( isolve.EQ.2 .AND. iround.EQ.1 ) THEN
600  IF( notran ) THEN
601  ifunc = ijob
602  END IF
603  scale2 = scale
604  CALL zlacpy( 'F', m, n, c, ldc, work, m )
605  CALL zlacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
606  CALL zlaset( 'F', m, n, czero, czero, c, ldc )
607  CALL zlaset( 'F', m, n, czero, czero, f, ldf )
608  ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
609  CALL zlacpy( 'F', m, n, work, m, c, ldc )
610  CALL zlacpy( 'F', m, n, work( m*n+1 ), m, f, ldf )
611  scale = scale2
612  END IF
613  150 CONTINUE
614  ELSE
615 *
616 * Solve transposed (I, J)-subsystem
617 * A(I, I)**H * R(I, J) + D(I, I)**H * L(I, J) = C(I, J)
618 * R(I, J) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
619 * for I = 1,2,..., P; J = Q, Q-1,..., 1
620 *
621  scale = one
622  DO 210 i = 1, p
623  is = iwork( i )
624  ie = iwork( i+1 ) - 1
625  mb = ie - is + 1
626  DO 200 j = q, p + 2, -1
627  js = iwork( j )
628  je = iwork( j+1 ) - 1
629  nb = je - js + 1
630  CALL ztgsy2( trans, ifunc, mb, nb, a( is, is ), lda,
631  $ b( js, js ), ldb, c( is, js ), ldc,
632  $ d( is, is ), ldd, e( js, js ), lde,
633  $ f( is, js ), ldf, scaloc, dsum, dscale,
634  $ linfo )
635  IF( linfo.GT.0 )
636  $ info = linfo
637  IF( scaloc.NE.one ) THEN
638  DO 160 k = 1, js - 1
639  CALL zscal( m, dcmplx( scaloc, zero ), c( 1, k ),
640  $ 1 )
641  CALL zscal( m, dcmplx( scaloc, zero ), f( 1, k ),
642  $ 1 )
643  160 CONTINUE
644  DO 170 k = js, je
645  CALL zscal( is-1, dcmplx( scaloc, zero ),
646  $ c( 1, k ), 1 )
647  CALL zscal( is-1, dcmplx( scaloc, zero ),
648  $ f( 1, k ), 1 )
649  170 CONTINUE
650  DO 180 k = js, je
651  CALL zscal( m-ie, dcmplx( scaloc, zero ),
652  $ c( ie+1, k ), 1 )
653  CALL zscal( m-ie, dcmplx( scaloc, zero ),
654  $ f( ie+1, k ), 1 )
655  180 CONTINUE
656  DO 190 k = je + 1, n
657  CALL zscal( m, dcmplx( scaloc, zero ), c( 1, k ),
658  $ 1 )
659  CALL zscal( m, dcmplx( scaloc, zero ), f( 1, k ),
660  $ 1 )
661  190 CONTINUE
662  scale = scale*scaloc
663  END IF
664 *
665 * Substitute R(I,J) and L(I,J) into remaining equation.
666 *
667  IF( j.GT.p+2 ) THEN
668  CALL zgemm( 'N', 'C', mb, js-1, nb,
669  $ dcmplx( one, zero ), c( is, js ), ldc,
670  $ b( 1, js ), ldb, dcmplx( one, zero ),
671  $ f( is, 1 ), ldf )
672  CALL zgemm( 'N', 'C', mb, js-1, nb,
673  $ dcmplx( one, zero ), f( is, js ), ldf,
674  $ e( 1, js ), lde, dcmplx( one, zero ),
675  $ f( is, 1 ), ldf )
676  END IF
677  IF( i.LT.p ) THEN
678  CALL zgemm( 'C', 'N', m-ie, nb, mb,
679  $ dcmplx( -one, zero ), a( is, ie+1 ), lda,
680  $ c( is, js ), ldc, dcmplx( one, zero ),
681  $ c( ie+1, js ), ldc )
682  CALL zgemm( 'C', 'N', m-ie, nb, mb,
683  $ dcmplx( -one, zero ), d( is, ie+1 ), ldd,
684  $ f( is, js ), ldf, dcmplx( one, zero ),
685  $ c( ie+1, js ), ldc )
686  END IF
687  200 CONTINUE
688  210 CONTINUE
689  END IF
690 *
691  work( 1 ) = lwmin
692 *
693  RETURN
694 *
695 * End of ZTGSYL
696 *
697  END
ztgsy2
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:261
zgemm
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
ztgsyl
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:297
zlaset
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:108
zlacpy
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
zscal
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:80