Annotation of rpl/lapack/lapack/ztgsy2.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b ZTGSY2
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZTGSY2 + dependencies
! 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 November 2011
! 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: * =====================================================================
1.1 bertrand 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: *
1.9 ! bertrand 262: * -- LAPACK auxiliary routine (version 3.4.0) --
1.1 bertrand 263: * -- LAPACK is a software package provided by Univ. of Tennessee, --
264: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 265: * November 2011
1.1 bertrand 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 = -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( '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:
1.8 bertrand 411: * A(I, I)**H * R(I, J) + D(I, I)**H * L(J, J) = C(I, J)
1.1 bertrand 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: *
1.8 bertrand 420: * Build 2 by 2 system Z**H
1.1 bertrand 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: *
1.8 bertrand 433: * Solve Z**H * x = RHS
1.1 bertrand 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ztgsy2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack