LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
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
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: claset.f:108
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189
subroutine ctgsy2(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, INFO)
CTGSY2 solves the generalized Sylvester equation (unblocked algorithm).
Definition: ctgsy2.f:261
subroutine ctgsyl(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, IWORK, INFO)
CTGSYL
Definition: ctgsyl.f:297
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
logical function lde(RI, RJ, LR)
Definition: dblat2.f:2945
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:54