LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
ctgsy2.f
Go to the documentation of this file.
1 *> \brief \b CTGSY2 solves the generalized Sylvester equation (unblocked algorithm).
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CTGSY2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctgsy2.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctgsy2.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctgsy2.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CTGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
22 * LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL,
23 * INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER TRANS
27 * INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N
28 * REAL RDSCAL, RDSUM, SCALE
29 * ..
30 * .. Array Arguments ..
31 * COMPLEX A( LDA, * ), B( LDB, * ), C( LDC, * ),
32 * $ D( LDD, * ), E( LDE, * ), F( LDF, * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> CTGSY2 solves the generalized Sylvester equation
42 *>
43 *> A * R - L * B = scale * C (1)
44 *> D * R - L * E = scale * F
45 *>
46 *> using Level 1 and 2 BLAS, where R and L are unknown M-by-N matrices,
47 *> (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M,
48 *> N-by-N and M-by-N, respectively. A, B, D and E are upper triangular
49 *> (i.e., (A,D) and (B,E) in generalized Schur form).
50 *>
51 *> The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an output
52 *> scaling factor chosen to avoid overflow.
53 *>
54 *> In matrix notation solving equation (1) corresponds to solve
55 *> Zx = scale * b, where Z is defined as
56 *>
57 *> Z = [ kron(In, A) -kron(B**H, Im) ] (2)
58 *> [ kron(In, D) -kron(E**H, Im) ],
59 *>
60 *> Ik is the identity matrix of size k and X**H is the transpose of X.
61 *> kron(X, Y) is the Kronecker product between the matrices X and Y.
62 *>
63 *> If TRANS = 'C', y in the conjugate transposed system Z**H*y = scale*b
64 *> is solved for, which is equivalent to solve for R and L in
65 *>
66 *> A**H * R + D**H * L = scale * C (3)
67 *> R * B**H + L * E**H = scale * -F
68 *>
69 *> This case is used to compute an estimate of Dif[(A, D), (B, E)] =
70 *> = sigma_min(Z) using reverse communicaton with CLACON.
71 *>
72 *> CTGSY2 also (IJOB >= 1) contributes to the computation in CTGSYL
73 *> of an upper bound on the separation between to matrix pairs. Then
74 *> the input (A, D), (B, E) are sub-pencils of two matrix pairs in
75 *> CTGSYL.
76 *> \endverbatim
77 *
78 * Arguments:
79 * ==========
80 *
81 *> \param[in] TRANS
82 *> \verbatim
83 *> TRANS is CHARACTER*1
84 *> = 'N', solve the generalized Sylvester equation (1).
85 *> = 'T': solve the 'transposed' system (3).
86 *> \endverbatim
87 *>
88 *> \param[in] IJOB
89 *> \verbatim
90 *> IJOB is INTEGER
91 *> Specifies what kind of functionality to be performed.
92 *> =0: solve (1) only.
93 *> =1: A contribution from this subsystem to a Frobenius
94 *> norm-based estimate of the separation between two matrix
95 *> pairs is computed. (look ahead strategy is used).
96 *> =2: A contribution from this subsystem to a Frobenius
97 *> norm-based estimate of the separation between two matrix
98 *> pairs is computed. (SGECON on sub-systems is used.)
99 *> Not referenced if TRANS = 'T'.
100 *> \endverbatim
101 *>
102 *> \param[in] M
103 *> \verbatim
104 *> M is INTEGER
105 *> On entry, M specifies the order of A and D, and the row
106 *> dimension of C, F, R and L.
107 *> \endverbatim
108 *>
109 *> \param[in] N
110 *> \verbatim
111 *> N is INTEGER
112 *> On entry, N specifies the order of B and E, and the column
113 *> dimension of C, F, R and L.
114 *> \endverbatim
115 *>
116 *> \param[in] A
117 *> \verbatim
118 *> A is COMPLEX array, dimension (LDA, M)
119 *> On entry, A contains an upper triangular matrix.
120 *> \endverbatim
121 *>
122 *> \param[in] LDA
123 *> \verbatim
124 *> LDA is INTEGER
125 *> The leading dimension of the matrix A. LDA >= max(1, M).
126 *> \endverbatim
127 *>
128 *> \param[in] B
129 *> \verbatim
130 *> B is COMPLEX array, dimension (LDB, N)
131 *> On entry, B contains an upper triangular matrix.
132 *> \endverbatim
133 *>
134 *> \param[in] LDB
135 *> \verbatim
136 *> LDB is INTEGER
137 *> The leading dimension of the matrix B. LDB >= max(1, N).
138 *> \endverbatim
139 *>
140 *> \param[in,out] C
141 *> \verbatim
142 *> C is COMPLEX array, dimension (LDC, N)
143 *> On entry, C contains the right-hand-side of the first matrix
144 *> equation in (1).
145 *> On exit, if IJOB = 0, C has been overwritten by the solution
146 *> R.
147 *> \endverbatim
148 *>
149 *> \param[in] LDC
150 *> \verbatim
151 *> LDC is INTEGER
152 *> The leading dimension of the matrix C. LDC >= max(1, M).
153 *> \endverbatim
154 *>
155 *> \param[in] D
156 *> \verbatim
157 *> D is COMPLEX array, dimension (LDD, M)
158 *> On entry, D contains an upper triangular matrix.
159 *> \endverbatim
160 *>
161 *> \param[in] LDD
162 *> \verbatim
163 *> LDD is INTEGER
164 *> The leading dimension of the matrix D. LDD >= max(1, M).
165 *> \endverbatim
166 *>
167 *> \param[in] E
168 *> \verbatim
169 *> E is COMPLEX array, dimension (LDE, N)
170 *> On entry, E contains an upper triangular matrix.
171 *> \endverbatim
172 *>
173 *> \param[in] LDE
174 *> \verbatim
175 *> LDE is INTEGER
176 *> The leading dimension of the matrix E. LDE >= max(1, N).
177 *> \endverbatim
178 *>
179 *> \param[in,out] F
180 *> \verbatim
181 *> F is COMPLEX array, dimension (LDF, N)
182 *> On entry, F contains the right-hand-side of the second matrix
183 *> equation in (1).
184 *> On exit, if IJOB = 0, F has been overwritten by the solution
185 *> L.
186 *> \endverbatim
187 *>
188 *> \param[in] LDF
189 *> \verbatim
190 *> LDF is INTEGER
191 *> The leading dimension of the matrix F. LDF >= max(1, M).
192 *> \endverbatim
193 *>
194 *> \param[out] SCALE
195 *> \verbatim
196 *> SCALE is REAL
197 *> On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions
198 *> R and L (C and F on entry) will hold the solutions to a
199 *> slightly perturbed system but the input matrices A, B, D and
200 *> E have not been changed. If SCALE = 0, R and L will hold the
201 *> solutions to the homogeneous system with C = F = 0.
202 *> Normally, SCALE = 1.
203 *> \endverbatim
204 *>
205 *> \param[in,out] RDSUM
206 *> \verbatim
207 *> RDSUM is REAL
208 *> On entry, the sum of squares of computed contributions to
209 *> the Dif-estimate under computation by CTGSYL, where the
210 *> scaling factor RDSCAL (see below) has been factored out.
211 *> On exit, the corresponding sum of squares updated with the
212 *> contributions from the current sub-system.
213 *> If TRANS = 'T' RDSUM is not touched.
214 *> NOTE: RDSUM only makes sense when CTGSY2 is called by
215 *> CTGSYL.
216 *> \endverbatim
217 *>
218 *> \param[in,out] RDSCAL
219 *> \verbatim
220 *> RDSCAL is REAL
221 *> On entry, scaling factor used to prevent overflow in RDSUM.
222 *> On exit, RDSCAL is updated w.r.t. the current contributions
223 *> in RDSUM.
224 *> If TRANS = 'T', RDSCAL is not touched.
225 *> NOTE: RDSCAL only makes sense when CTGSY2 is called by
226 *> CTGSYL.
227 *> \endverbatim
228 *>
229 *> \param[out] INFO
230 *> \verbatim
231 *> INFO is INTEGER
232 *> On exit, if INFO is set to
233 *> =0: Successful exit
234 *> <0: If INFO = -i, input argument number i is illegal.
235 *> >0: The matrix pairs (A, D) and (B, E) have common or very
236 *> close eigenvalues.
237 *> \endverbatim
238 *
239 * Authors:
240 * ========
241 *
242 *> \author Univ. of Tennessee
243 *> \author Univ. of California Berkeley
244 *> \author Univ. of Colorado Denver
245 *> \author NAG Ltd.
246 *
247 *> \date September 2012
248 *
249 *> \ingroup complexSYauxiliary
250 *
251 *> \par Contributors:
252 * ==================
253 *>
254 *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
255 *> Umea University, S-901 87 Umea, Sweden.
256 *
257 * =====================================================================
258  SUBROUTINE ctgsy2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
259  $ ldd, e, lde, f, ldf, scale, rdsum, rdscal,
260  $ info )
261 *
262 * -- LAPACK auxiliary routine (version 3.4.2) --
263 * -- LAPACK is a software package provided by Univ. of Tennessee, --
264 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
265 * September 2012
266 *
267 * .. Scalar Arguments ..
268  CHARACTER trans
269  INTEGER ijob, info, lda, ldb, ldc, ldd, lde, ldf, m, n
270  REAL rdscal, rdsum, scale
271 * ..
272 * .. Array Arguments ..
273  COMPLEX a( lda, * ), b( ldb, * ), c( ldc, * ),
274  $ d( ldd, * ), e( lde, * ), f( ldf, * )
275 * ..
276 *
277 * =====================================================================
278 *
279 * .. Parameters ..
280  REAL zero, one
281  INTEGER ldz
282  parameter( zero = 0.0e+0, one = 1.0e+0, ldz = 2 )
283 * ..
284 * .. Local Scalars ..
285  LOGICAL notran
286  INTEGER i, ierr, j, k
287  REAL scaloc
288  COMPLEX alpha
289 * ..
290 * .. Local Arrays ..
291  INTEGER ipiv( ldz ), jpiv( ldz )
292  COMPLEX rhs( ldz ), z( ldz, ldz )
293 * ..
294 * .. External Functions ..
295  LOGICAL lsame
296  EXTERNAL lsame
297 * ..
298 * .. External Subroutines ..
299  EXTERNAL caxpy, cgesc2, cgetc2, cscal, clatdf, xerbla
300 * ..
301 * .. Intrinsic Functions ..
302  INTRINSIC cmplx, conjg, max
303 * ..
304 * .. Executable Statements ..
305 *
306 * Decode and test input parameters
307 *
308  info = 0
309  ierr = 0
310  notran = lsame( trans, 'N' )
311  IF( .NOT.notran .AND. .NOT.lsame( trans, 'C' ) ) THEN
312  info = -1
313  ELSE IF( notran ) THEN
314  IF( ( ijob.LT.0 ) .OR. ( ijob.GT.2 ) ) THEN
315  info = -2
316  END IF
317  END IF
318  IF( info.EQ.0 ) THEN
319  IF( m.LE.0 ) THEN
320  info = -3
321  ELSE IF( n.LE.0 ) THEN
322  info = -4
323  ELSE IF( lda.LT.max( 1, m ) ) THEN
324  info = -5
325  ELSE IF( ldb.LT.max( 1, n ) ) THEN
326  info = -8
327  ELSE IF( ldc.LT.max( 1, m ) ) THEN
328  info = -10
329  ELSE IF( ldd.LT.max( 1, m ) ) THEN
330  info = -12
331  ELSE IF( lde.LT.max( 1, n ) ) THEN
332  info = -14
333  ELSE IF( ldf.LT.max( 1, m ) ) THEN
334  info = -16
335  END IF
336  END IF
337  IF( info.NE.0 ) THEN
338  CALL xerbla( 'CTGSY2', -info )
339  return
340  END IF
341 *
342  IF( notran ) THEN
343 *
344 * Solve (I, J) - system
345 * A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
346 * D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
347 * for I = M, M - 1, ..., 1; J = 1, 2, ..., N
348 *
349  scale = one
350  scaloc = one
351  DO 30 j = 1, n
352  DO 20 i = m, 1, -1
353 *
354 * Build 2 by 2 system
355 *
356  z( 1, 1 ) = a( i, i )
357  z( 2, 1 ) = d( i, i )
358  z( 1, 2 ) = -b( j, j )
359  z( 2, 2 ) = -e( j, j )
360 *
361 * Set up right hand side(s)
362 *
363  rhs( 1 ) = c( i, j )
364  rhs( 2 ) = f( i, j )
365 *
366 * Solve Z * x = RHS
367 *
368  CALL cgetc2( ldz, z, ldz, ipiv, jpiv, ierr )
369  IF( ierr.GT.0 )
370  $ info = ierr
371  IF( ijob.EQ.0 ) THEN
372  CALL cgesc2( ldz, z, ldz, rhs, ipiv, jpiv, scaloc )
373  IF( scaloc.NE.one ) THEN
374  DO 10 k = 1, n
375  CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
376  $ 1 )
377  CALL cscal( m, cmplx( scaloc, zero ), f( 1, k ),
378  $ 1 )
379  10 continue
380  scale = scale*scaloc
381  END IF
382  ELSE
383  CALL clatdf( ijob, ldz, z, ldz, rhs, rdsum, rdscal,
384  $ ipiv, jpiv )
385  END IF
386 *
387 * Unpack solution vector(s)
388 *
389  c( i, j ) = rhs( 1 )
390  f( i, j ) = rhs( 2 )
391 *
392 * Substitute R(I, J) and L(I, J) into remaining equation.
393 *
394  IF( i.GT.1 ) THEN
395  alpha = -rhs( 1 )
396  CALL caxpy( i-1, alpha, a( 1, i ), 1, c( 1, j ), 1 )
397  CALL caxpy( i-1, alpha, d( 1, i ), 1, f( 1, j ), 1 )
398  END IF
399  IF( j.LT.n ) THEN
400  CALL caxpy( n-j, rhs( 2 ), b( j, j+1 ), ldb,
401  $ c( i, j+1 ), ldc )
402  CALL caxpy( n-j, rhs( 2 ), e( j, j+1 ), lde,
403  $ f( i, j+1 ), ldf )
404  END IF
405 *
406  20 continue
407  30 continue
408  ELSE
409 *
410 * Solve transposed (I, J) - system:
411 * A(I, I)**H * R(I, J) + D(I, I)**H * L(J, J) = C(I, J)
412 * R(I, I) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
413 * for I = 1, 2, ..., M, J = N, N - 1, ..., 1
414 *
415  scale = one
416  scaloc = one
417  DO 80 i = 1, m
418  DO 70 j = n, 1, -1
419 *
420 * Build 2 by 2 system Z**H
421 *
422  z( 1, 1 ) = conjg( a( i, i ) )
423  z( 2, 1 ) = -conjg( b( j, j ) )
424  z( 1, 2 ) = conjg( d( i, i ) )
425  z( 2, 2 ) = -conjg( e( j, j ) )
426 *
427 *
428 * Set up right hand side(s)
429 *
430  rhs( 1 ) = c( i, j )
431  rhs( 2 ) = f( i, j )
432 *
433 * Solve Z**H * x = RHS
434 *
435  CALL cgetc2( ldz, z, ldz, ipiv, jpiv, ierr )
436  IF( ierr.GT.0 )
437  $ info = ierr
438  CALL cgesc2( ldz, z, ldz, rhs, ipiv, jpiv, scaloc )
439  IF( scaloc.NE.one ) THEN
440  DO 40 k = 1, n
441  CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
442  $ 1 )
443  CALL cscal( m, cmplx( scaloc, zero ), f( 1, k ),
444  $ 1 )
445  40 continue
446  scale = scale*scaloc
447  END IF
448 *
449 * Unpack solution vector(s)
450 *
451  c( i, j ) = rhs( 1 )
452  f( i, j ) = rhs( 2 )
453 *
454 * Substitute R(I, J) and L(I, J) into remaining equation.
455 *
456  DO 50 k = 1, j - 1
457  f( i, k ) = f( i, k ) + rhs( 1 )*conjg( b( k, j ) ) +
458  $ rhs( 2 )*conjg( e( k, j ) )
459  50 continue
460  DO 60 k = i + 1, m
461  c( k, j ) = c( k, j ) - conjg( a( i, k ) )*rhs( 1 ) -
462  $ conjg( d( i, k ) )*rhs( 2 )
463  60 continue
464 *
465  70 continue
466  80 continue
467  END IF
468  return
469 *
470 * End of CTGSY2
471 *
472  END