LAPACK  3.10.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 homogeneous 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 *> \ingroup complexSYcomputational
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 ctgsyl( 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  REAL DIF, SCALE
305 * ..
306 * .. Array Arguments ..
307  INTEGER IWORK( * )
308  COMPLEX 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  REAL ZERO, ONE
319  PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0 )
320  COMPLEX CZERO
321  parameter( czero = (0.0e+0, 0.0e+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  REAL DSCALE, DSUM, SCALE2, SCALOC
328 * ..
329 * .. External Functions ..
330  LOGICAL LSAME
331  INTEGER ILAENV
332  EXTERNAL lsame, ilaenv
333 * ..
334 * .. External Subroutines ..
335  EXTERNAL cgemm, clacpy, claset, cscal, ctgsy2, xerbla
336 * ..
337 * .. Intrinsic Functions ..
338  INTRINSIC cmplx, max, real, 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( 'CTGSYL', -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, 'CTGSYL', trans, m, n, -1, -1 )
414  nb = ilaenv( 5, 'CTGSYL', 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 claset( 'F', m, n, czero, czero, c, ldc )
422  CALL claset( '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 ctgsy2( 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( real( 2*m*n ) ) / ( dscale*sqrt( dsum ) )
445  ELSE
446  dif = sqrt( real( 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 clacpy( 'F', m, n, c, ldc, work, m )
455  CALL clacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
456  CALL claset( 'F', m, n, czero, czero, c, ldc )
457  CALL claset( 'F', m, n, czero, czero, f, ldf )
458  ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
459  CALL clacpy( 'F', m, n, work, m, c, ldc )
460  CALL clacpy( '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 ctgsy2( 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 cscal( m, cmplx( scaloc, zero ), c( 1, k ),
538  $ 1 )
539  CALL cscal( m, cmplx( scaloc, zero ), f( 1, k ),
540  $ 1 )
541  80 CONTINUE
542  DO 90 k = js, je
543  CALL cscal( is-1, cmplx( scaloc, zero ),
544  $ c( 1, k ), 1 )
545  CALL cscal( is-1, cmplx( scaloc, zero ),
546  $ f( 1, k ), 1 )
547  90 CONTINUE
548  DO 100 k = js, je
549  CALL cscal( m-ie, cmplx( scaloc, zero ),
550  $ c( ie+1, k ), 1 )
551  CALL cscal( m-ie, cmplx( scaloc, zero ),
552  $ f( ie+1, k ), 1 )
553  100 CONTINUE
554  DO 110 k = je + 1, n
555  CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
556  $ 1 )
557  CALL cscal( m, cmplx( scaloc, zero ), f( 1, k ),
558  $ 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 cgemm( 'N', 'N', is-1, nb, mb,
567  $ cmplx( -one, zero ), a( 1, is ), lda,
568  $ c( is, js ), ldc, cmplx( one, zero ),
569  $ c( 1, js ), ldc )
570  CALL cgemm( 'N', 'N', is-1, nb, mb,
571  $ cmplx( -one, zero ), d( 1, is ), ldd,
572  $ c( is, js ), ldc, cmplx( one, zero ),
573  $ f( 1, js ), ldf )
574  END IF
575  IF( j.LT.q ) THEN
576  CALL cgemm( 'N', 'N', mb, n-je, nb,
577  $ cmplx( one, zero ), f( is, js ), ldf,
578  $ b( js, je+1 ), ldb, cmplx( one, zero ),
579  $ c( is, je+1 ), ldc )
580  CALL cgemm( 'N', 'N', mb, n-je, nb,
581  $ cmplx( one, zero ), f( is, js ), ldf,
582  $ e( js, je+1 ), lde, cmplx( one, zero ),
583  $ f( is, je+1 ), ldf )
584  END IF
585  120 CONTINUE
586  130 CONTINUE
587  IF( dscale.NE.zero ) THEN
588  IF( ijob.EQ.1 .OR. ijob.EQ.3 ) THEN
589  dif = sqrt( real( 2*m*n ) ) / ( dscale*sqrt( dsum ) )
590  ELSE
591  dif = sqrt( real( pq ) ) / ( dscale*sqrt( dsum ) )
592  END IF
593  END IF
594  IF( isolve.EQ.2 .AND. iround.EQ.1 ) THEN
595  IF( notran ) THEN
596  ifunc = ijob
597  END IF
598  scale2 = scale
599  CALL clacpy( 'F', m, n, c, ldc, work, m )
600  CALL clacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
601  CALL claset( 'F', m, n, czero, czero, c, ldc )
602  CALL claset( 'F', m, n, czero, czero, f, ldf )
603  ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
604  CALL clacpy( 'F', m, n, work, m, c, ldc )
605  CALL clacpy( 'F', m, n, work( m*n+1 ), m, f, ldf )
606  scale = scale2
607  END IF
608  150 CONTINUE
609  ELSE
610 *
611 * Solve transposed (I, J)-subsystem
612 * A(I, I)**H * R(I, J) + D(I, I)**H * L(I, J) = C(I, J)
613 * R(I, J) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
614 * for I = 1,2,..., P; J = Q, Q-1,..., 1
615 *
616  scale = one
617  DO 210 i = 1, p
618  is = iwork( i )
619  ie = iwork( i+1 ) - 1
620  mb = ie - is + 1
621  DO 200 j = q, p + 2, -1
622  js = iwork( j )
623  je = iwork( j+1 ) - 1
624  nb = je - js + 1
625  CALL ctgsy2( trans, ifunc, mb, nb, a( is, is ), lda,
626  $ b( js, js ), ldb, c( is, js ), ldc,
627  $ d( is, is ), ldd, e( js, js ), lde,
628  $ f( is, js ), ldf, scaloc, dsum, dscale,
629  $ linfo )
630  IF( linfo.GT.0 )
631  $ info = linfo
632  IF( scaloc.NE.one ) THEN
633  DO 160 k = 1, js - 1
634  CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
635  $ 1 )
636  CALL cscal( m, cmplx( scaloc, zero ), f( 1, k ),
637  $ 1 )
638  160 CONTINUE
639  DO 170 k = js, je
640  CALL cscal( is-1, cmplx( scaloc, zero ), c( 1, k ),
641  $ 1 )
642  CALL cscal( is-1, cmplx( scaloc, zero ), f( 1, k ),
643  $ 1 )
644  170 CONTINUE
645  DO 180 k = js, je
646  CALL cscal( m-ie, cmplx( scaloc, zero ),
647  $ c( ie+1, k ), 1 )
648  CALL cscal( m-ie, cmplx( scaloc, zero ),
649  $ f( ie+1, k ), 1 )
650  180 CONTINUE
651  DO 190 k = je + 1, n
652  CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
653  $ 1 )
654  CALL cscal( m, cmplx( scaloc, zero ), f( 1, k ),
655  $ 1 )
656  190 CONTINUE
657  scale = scale*scaloc
658  END IF
659 *
660 * Substitute R(I,J) and L(I,J) into remaining equation.
661 *
662  IF( j.GT.p+2 ) THEN
663  CALL cgemm( 'N', 'C', mb, js-1, nb,
664  $ cmplx( one, zero ), c( is, js ), ldc,
665  $ b( 1, js ), ldb, cmplx( one, zero ),
666  $ f( is, 1 ), ldf )
667  CALL cgemm( 'N', 'C', mb, js-1, nb,
668  $ cmplx( one, zero ), f( is, js ), ldf,
669  $ e( 1, js ), lde, cmplx( one, zero ),
670  $ f( is, 1 ), ldf )
671  END IF
672  IF( i.LT.p ) THEN
673  CALL cgemm( 'C', 'N', m-ie, nb, mb,
674  $ cmplx( -one, zero ), a( is, ie+1 ), lda,
675  $ c( is, js ), ldc, cmplx( one, zero ),
676  $ c( ie+1, js ), ldc )
677  CALL cgemm( 'C', 'N', m-ie, nb, mb,
678  $ cmplx( -one, zero ), d( is, ie+1 ), ldd,
679  $ f( is, js ), ldf, cmplx( one, zero ),
680  $ c( ie+1, js ), ldc )
681  END IF
682  200 CONTINUE
683  210 CONTINUE
684  END IF
685 *
686  work( 1 ) = lwmin
687 *
688  RETURN
689 *
690 * End of CTGSYL
691 *
692  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:78
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:187
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:106
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:103
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:259
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:295