SRC/ztgsy2.f

   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
   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 2015
 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.6.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 *     November 2015
 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