LAPACK  3.8.0 LAPACK: Linear Algebra PACKage
ztgsy2.f
Go to the documentation of this file.
1 *> \brief \b ZTGSY2 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
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztgsy2.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztgsy2.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztgsy2.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZTGSY2( 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 * DOUBLE PRECISION RDSCAL, RDSUM, SCALE
29 * ..
30 * .. Array Arguments ..
31 * COMPLEX*16 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 *> ZTGSY2 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 conjuguate 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 ZLACON.
71 *>
72 *> ZTGSY2 also (IJOB >= 1) contributes to the computation in ZTGSYL
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 *> ZTGSYL.
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. (DGECON 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*16 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*16 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*16 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*16 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*16 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*16 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 DOUBLE PRECISION
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 DOUBLE PRECISION
208 *> On entry, the sum of squares of computed contributions to
209 *> the Dif-estimate under computation by ZTGSYL, 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 ZTGSY2 is called by
215 *> ZTGSYL.
216 *> \endverbatim
217 *>
218 *> \param[in,out] RDSCAL
219 *> \verbatim
220 *> RDSCAL is DOUBLE PRECISION
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 ZTGSY2 is called by
226 *> ZTGSYL.
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 December 2016
248 *
249 *> \ingroup complex16SYauxiliary
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 ztgsy2( 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.7.0) --
263 * -- LAPACK is a software package provided by Univ. of Tennessee, --
264 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
265 * December 2016
266 *
267 * .. Scalar Arguments ..
268  CHARACTER TRANS
269  INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N
270  DOUBLE PRECISION RDSCAL, RDSUM, SCALE
271 * ..
272 * .. Array Arguments ..
273  COMPLEX*16 A( lda, * ), B( ldb, * ), C( ldc, * ),
274  \$ d( ldd, * ), e( lde, * ), f( ldf, * )
275 * ..
276 *
277 * =====================================================================
278 *
279 * .. Parameters ..
280  DOUBLE PRECISION ZERO, ONE
281  INTEGER LDZ
282  parameter( zero = 0.0d+0, one = 1.0d+0, ldz = 2 )
283 * ..
284 * .. Local Scalars ..
285  LOGICAL NOTRAN
286  INTEGER I, IERR, J, K
287  DOUBLE PRECISION SCALOC
288  COMPLEX*16 ALPHA
289 * ..
290 * .. Local Arrays ..
291  INTEGER IPIV( ldz ), JPIV( ldz )
292  COMPLEX*16 RHS( ldz ), Z( ldz, ldz )
293 * ..
294 * .. External Functions ..
295  LOGICAL LSAME
296  EXTERNAL lsame
297 * ..
298 * .. External Subroutines ..
299  EXTERNAL xerbla, zaxpy, zgesc2, zgetc2, zlatdf, zscal
300 * ..
301 * .. Intrinsic Functions ..
302  INTRINSIC dcmplx, dconjg, 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 = -6
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( 'ZTGSY2', -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 zgetc2( ldz, z, ldz, ipiv, jpiv, ierr )
369  IF( ierr.GT.0 )
370  \$ info = ierr
371  IF( ijob.EQ.0 ) THEN
372  CALL zgesc2( ldz, z, ldz, rhs, ipiv, jpiv, scaloc )
373  IF( scaloc.NE.one ) THEN
374  DO 10 k = 1, n
375  CALL zscal( m, dcmplx( scaloc, zero ),
376  \$ c( 1, k ), 1 )
377  CALL zscal( m, dcmplx( scaloc, zero ),
378  \$ f( 1, k ), 1 )
379  10 CONTINUE
380  scale = scale*scaloc
381  END IF
382  ELSE
383  CALL zlatdf( 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 zaxpy( i-1, alpha, a( 1, i ), 1, c( 1, j ), 1 )
397  CALL zaxpy( i-1, alpha, d( 1, i ), 1, f( 1, j ), 1 )
398  END IF
399  IF( j.LT.n ) THEN
400  CALL zaxpy( n-j, rhs( 2 ), b( j, j+1 ), ldb,
401  \$ c( i, j+1 ), ldc )
402  CALL zaxpy( 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 ) = dconjg( a( i, i ) )
423  z( 2, 1 ) = -dconjg( b( j, j ) )
424  z( 1, 2 ) = dconjg( d( i, i ) )
425  z( 2, 2 ) = -dconjg( 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 zgetc2( ldz, z, ldz, ipiv, jpiv, ierr )
436  IF( ierr.GT.0 )
437  \$ info = ierr
438  CALL zgesc2( ldz, z, ldz, rhs, ipiv, jpiv, scaloc )
439  IF( scaloc.NE.one ) THEN
440  DO 40 k = 1, n
441  CALL zscal( m, dcmplx( scaloc, zero ), c( 1, k ),
442  \$ 1 )
443  CALL zscal( m, dcmplx( 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 )*dconjg( b( k, j ) ) +
458  \$ rhs( 2 )*dconjg( e( k, j ) )
459  50 CONTINUE
460  DO 60 k = i + 1, m
461  c( k, j ) = c( k, j ) - dconjg( a( i, k ) )*rhs( 1 ) -
462  \$ dconjg( d( i, k ) )*rhs( 2 )
463  60 CONTINUE
464 *
465  70 CONTINUE
466  80 CONTINUE
467  END IF
468  RETURN
469 *
470 * End of ZTGSY2
471 *
472  END
subroutine zlatdf(IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV, JPIV)
ZLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution ...
Definition: zlatdf.f:171
subroutine zgetc2(N, A, LDA, IPIV, JPIV, INFO)
ZGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix...
Definition: zgetc2.f:113
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zgesc2(N, A, LDA, RHS, IPIV, JPIV, SCALE)
ZGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed...
Definition: zgesc2.f:117
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
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:90
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:80