LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
ctgsyl.f
Go to the documentation of this file.
1 *> \brief \b CTGSYL
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CTGSYL + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctgsyl.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctgsyl.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctgsyl.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CTGSYL( 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 * REAL DIF, SCALE
30 * ..
31 * .. Array Arguments ..
32 * INTEGER IWORK( * )
33 * COMPLEX 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 *> CTGSYL 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 CLACON.
76 *>
77 *> If IJOB >= 1, CTGSYL 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 *> (CGECON 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 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 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 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 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 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 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 REAL
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 REAL
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 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 November 2011
264 *
265 *> \ingroup complexSYcomputational
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 ctgsyl( 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.4.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 * November 2011
302 *
303 * .. Scalar Arguments ..
304  CHARACTER trans
305  INTEGER ijob, info, lda, ldb, ldc, ldd, lde, ldf,
306  $ lwork, m, n
307  REAL dif, scale
308 * ..
309 * .. Array Arguments ..
310  INTEGER iwork( * )
311  COMPLEX 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  REAL zero, one
322  parameter( zero = 0.0e+0, one = 1.0e+0 )
323  COMPLEX czero
324  parameter( czero = (0.0e+0, 0.0e+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  REAL dscale, dsum, scale2, scaloc
331 * ..
332 * .. External Functions ..
333  LOGICAL lsame
334  INTEGER ilaenv
335  EXTERNAL lsame, ilaenv
336 * ..
337 * .. External Subroutines ..
338  EXTERNAL cgemm, clacpy, claset, cscal, ctgsy2, xerbla
339 * ..
340 * .. Intrinsic Functions ..
341  INTRINSIC cmplx, max, REAL, 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( 'CTGSYL', -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, 'CTGSYL', trans, m, n, -1, -1 )
417  nb = ilaenv( 5, 'CTGSYL', 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 claset( 'F', m, n, czero, czero, c, ldc )
425  CALL claset( '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 ctgsy2( 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( REAL( 2*M*N ) ) / ( dscale*sqrt( dsum ) )
448  ELSE
449  dif = sqrt( REAL( 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 clacpy( 'F', m, n, c, ldc, work, m )
458  CALL clacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
459  CALL claset( 'F', m, n, czero, czero, c, ldc )
460  CALL claset( 'F', m, n, czero, czero, f, ldf )
461  ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
462  CALL clacpy( 'F', m, n, work, m, c, ldc )
463  CALL clacpy( '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 ctgsy2( 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 cscal( m, cmplx( scaloc, zero ), c( 1, k ),
541  $ 1 )
542  CALL cscal( m, cmplx( scaloc, zero ), f( 1, k ),
543  $ 1 )
544  80 CONTINUE
545  DO 90 k = js, je
546  CALL cscal( is-1, cmplx( scaloc, zero ),
547  $ c( 1, k ), 1 )
548  CALL cscal( is-1, cmplx( scaloc, zero ),
549  $ f( 1, k ), 1 )
550  90 CONTINUE
551  DO 100 k = js, je
552  CALL cscal( m-ie, cmplx( scaloc, zero ),
553  $ c( ie+1, k ), 1 )
554  CALL cscal( m-ie, cmplx( scaloc, zero ),
555  $ f( ie+1, k ), 1 )
556  100 CONTINUE
557  DO 110 k = je + 1, n
558  CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
559  $ 1 )
560  CALL cscal( m, cmplx( scaloc, zero ), f( 1, k ),
561  $ 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 cgemm( 'N', 'N', is-1, nb, mb,
570  $ cmplx( -one, zero ), a( 1, is ), lda,
571  $ c( is, js ), ldc, cmplx( one, zero ),
572  $ c( 1, js ), ldc )
573  CALL cgemm( 'N', 'N', is-1, nb, mb,
574  $ cmplx( -one, zero ), d( 1, is ), ldd,
575  $ c( is, js ), ldc, cmplx( one, zero ),
576  $ f( 1, js ), ldf )
577  END IF
578  IF( j.LT.q ) THEN
579  CALL cgemm( 'N', 'N', mb, n-je, nb,
580  $ cmplx( one, zero ), f( is, js ), ldf,
581  $ b( js, je+1 ), ldb, cmplx( one, zero ),
582  $ c( is, je+1 ), ldc )
583  CALL cgemm( 'N', 'N', mb, n-je, nb,
584  $ cmplx( one, zero ), f( is, js ), ldf,
585  $ e( js, je+1 ), lde, cmplx( one, zero ),
586  $ f( is, je+1 ), ldf )
587  END IF
588  120 CONTINUE
589  130 CONTINUE
590  IF( dscale.NE.zero ) THEN
591  IF( ijob.EQ.1 .OR. ijob.EQ.3 ) THEN
592  dif = sqrt( REAL( 2*M*N ) ) / ( dscale*sqrt( dsum ) )
593  ELSE
594  dif = sqrt( REAL( PQ ) ) / ( dscale*sqrt( dsum ) )
595  END IF
596  END IF
597  IF( isolve.EQ.2 .AND. iround.EQ.1 ) THEN
598  IF( notran ) THEN
599  ifunc = ijob
600  END IF
601  scale2 = scale
602  CALL clacpy( 'F', m, n, c, ldc, work, m )
603  CALL clacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
604  CALL claset( 'F', m, n, czero, czero, c, ldc )
605  CALL claset( 'F', m, n, czero, czero, f, ldf )
606  ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
607  CALL clacpy( 'F', m, n, work, m, c, ldc )
608  CALL clacpy( 'F', m, n, work( m*n+1 ), m, f, ldf )
609  scale = scale2
610  END IF
611  150 CONTINUE
612  ELSE
613 *
614 * Solve transposed (I, J)-subsystem
615 * A(I, I)**H * R(I, J) + D(I, I)**H * L(I, J) = C(I, J)
616 * R(I, J) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
617 * for I = 1,2,..., P; J = Q, Q-1,..., 1
618 *
619  scale = one
620  DO 210 i = 1, p
621  is = iwork( i )
622  ie = iwork( i+1 ) - 1
623  mb = ie - is + 1
624  DO 200 j = q, p + 2, -1
625  js = iwork( j )
626  je = iwork( j+1 ) - 1
627  nb = je - js + 1
628  CALL ctgsy2( trans, ifunc, mb, nb, a( is, is ), lda,
629  $ b( js, js ), ldb, c( is, js ), ldc,
630  $ d( is, is ), ldd, e( js, js ), lde,
631  $ f( is, js ), ldf, scaloc, dsum, dscale,
632  $ linfo )
633  IF( linfo.GT.0 )
634  $ info = linfo
635  IF( scaloc.NE.one ) THEN
636  DO 160 k = 1, js - 1
637  CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
638  $ 1 )
639  CALL cscal( m, cmplx( scaloc, zero ), f( 1, k ),
640  $ 1 )
641  160 CONTINUE
642  DO 170 k = js, je
643  CALL cscal( is-1, cmplx( scaloc, zero ), c( 1, k ),
644  $ 1 )
645  CALL cscal( is-1, cmplx( scaloc, zero ), f( 1, k ),
646  $ 1 )
647  170 CONTINUE
648  DO 180 k = js, je
649  CALL cscal( m-ie, cmplx( scaloc, zero ),
650  $ c( ie+1, k ), 1 )
651  CALL cscal( m-ie, cmplx( scaloc, zero ),
652  $ f( ie+1, k ), 1 )
653  180 CONTINUE
654  DO 190 k = je + 1, n
655  CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
656  $ 1 )
657  CALL cscal( m, cmplx( scaloc, zero ), f( 1, k ),
658  $ 1 )
659  190 CONTINUE
660  scale = scale*scaloc
661  END IF
662 *
663 * Substitute R(I,J) and L(I,J) into remaining equation.
664 *
665  IF( j.GT.p+2 ) THEN
666  CALL cgemm( 'N', 'C', mb, js-1, nb,
667  $ cmplx( one, zero ), c( is, js ), ldc,
668  $ b( 1, js ), ldb, cmplx( one, zero ),
669  $ f( is, 1 ), ldf )
670  CALL cgemm( 'N', 'C', mb, js-1, nb,
671  $ cmplx( one, zero ), f( is, js ), ldf,
672  $ e( 1, js ), lde, cmplx( one, zero ),
673  $ f( is, 1 ), ldf )
674  END IF
675  IF( i.LT.p ) THEN
676  CALL cgemm( 'C', 'N', m-ie, nb, mb,
677  $ cmplx( -one, zero ), a( is, ie+1 ), lda,
678  $ c( is, js ), ldc, cmplx( one, zero ),
679  $ c( ie+1, js ), ldc )
680  CALL cgemm( 'C', 'N', m-ie, nb, mb,
681  $ cmplx( -one, zero ), d( is, ie+1 ), ldd,
682  $ f( is, js ), ldf, cmplx( one, zero ),
683  $ c( ie+1, js ), ldc )
684  END IF
685  200 CONTINUE
686  210 CONTINUE
687  END IF
688 *
689  work( 1 ) = lwmin
690 *
691  RETURN
692 *
693 * End of CTGSYL
694 *
695  END