LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
ctgsy2.f
Go to the documentation of this file.
00001 *> \brief \b CTGSY2
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CTGSY2 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctgsy2.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctgsy2.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctgsy2.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CTGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
00022 *                          LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL,
00023 *                          INFO )
00024 * 
00025 *       .. Scalar Arguments ..
00026 *       CHARACTER          TRANS
00027 *       INTEGER            IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N
00028 *       REAL               RDSCAL, RDSUM, SCALE
00029 *       ..
00030 *       .. Array Arguments ..
00031 *       COMPLEX            A( LDA, * ), B( LDB, * ), C( LDC, * ),
00032 *      $                   D( LDD, * ), E( LDE, * ), F( LDF, * )
00033 *       ..
00034 *  
00035 *
00036 *> \par Purpose:
00037 *  =============
00038 *>
00039 *> \verbatim
00040 *>
00041 *> CTGSY2 solves the generalized Sylvester equation
00042 *>
00043 *>             A * R - L * B = scale *  C               (1)
00044 *>             D * R - L * E = scale * F
00045 *>
00046 *> using Level 1 and 2 BLAS, where R and L are unknown M-by-N matrices,
00047 *> (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M,
00048 *> N-by-N and M-by-N, respectively. A, B, D and E are upper triangular
00049 *> (i.e., (A,D) and (B,E) in generalized Schur form).
00050 *>
00051 *> The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an output
00052 *> scaling factor chosen to avoid overflow.
00053 *>
00054 *> In matrix notation solving equation (1) corresponds to solve
00055 *> Zx = scale * b, where Z is defined as
00056 *>
00057 *>        Z = [ kron(In, A)  -kron(B**H, Im) ]             (2)
00058 *>            [ kron(In, D)  -kron(E**H, Im) ],
00059 *>
00060 *> Ik is the identity matrix of size k and X**H is the transpose of X.
00061 *> kron(X, Y) is the Kronecker product between the matrices X and Y.
00062 *>
00063 *> If TRANS = 'C', y in the conjugate transposed system Z**H*y = scale*b
00064 *> is solved for, which is equivalent to solve for R and L in
00065 *>
00066 *>             A**H * R  + D**H * L   = scale * C           (3)
00067 *>             R  * B**H + L  * E**H  = scale * -F
00068 *>
00069 *> This case is used to compute an estimate of Dif[(A, D), (B, E)] =
00070 *> = sigma_min(Z) using reverse communicaton with CLACON.
00071 *>
00072 *> CTGSY2 also (IJOB >= 1) contributes to the computation in CTGSYL
00073 *> of an upper bound on the separation between to matrix pairs. Then
00074 *> the input (A, D), (B, E) are sub-pencils of two matrix pairs in
00075 *> CTGSYL.
00076 *> \endverbatim
00077 *
00078 *  Arguments:
00079 *  ==========
00080 *
00081 *> \param[in] TRANS
00082 *> \verbatim
00083 *>          TRANS is CHARACTER*1
00084 *>          = 'N', solve the generalized Sylvester equation (1).
00085 *>          = 'T': solve the 'transposed' system (3).
00086 *> \endverbatim
00087 *>
00088 *> \param[in] IJOB
00089 *> \verbatim
00090 *>          IJOB is INTEGER
00091 *>          Specifies what kind of functionality to be performed.
00092 *>          =0: solve (1) only.
00093 *>          =1: A contribution from this subsystem to a Frobenius
00094 *>              norm-based estimate of the separation between two matrix
00095 *>              pairs is computed. (look ahead strategy is used).
00096 *>          =2: A contribution from this subsystem to a Frobenius
00097 *>              norm-based estimate of the separation between two matrix
00098 *>              pairs is computed. (SGECON on sub-systems is used.)
00099 *>          Not referenced if TRANS = 'T'.
00100 *> \endverbatim
00101 *>
00102 *> \param[in] M
00103 *> \verbatim
00104 *>          M is INTEGER
00105 *>          On entry, M specifies the order of A and D, and the row
00106 *>          dimension of C, F, R and L.
00107 *> \endverbatim
00108 *>
00109 *> \param[in] N
00110 *> \verbatim
00111 *>          N is INTEGER
00112 *>          On entry, N specifies the order of B and E, and the column
00113 *>          dimension of C, F, R and L.
00114 *> \endverbatim
00115 *>
00116 *> \param[in] A
00117 *> \verbatim
00118 *>          A is COMPLEX array, dimension (LDA, M)
00119 *>          On entry, A contains an upper triangular matrix.
00120 *> \endverbatim
00121 *>
00122 *> \param[in] LDA
00123 *> \verbatim
00124 *>          LDA is INTEGER
00125 *>          The leading dimension of the matrix A. LDA >= max(1, M).
00126 *> \endverbatim
00127 *>
00128 *> \param[in] B
00129 *> \verbatim
00130 *>          B is COMPLEX array, dimension (LDB, N)
00131 *>          On entry, B contains an upper triangular matrix.
00132 *> \endverbatim
00133 *>
00134 *> \param[in] LDB
00135 *> \verbatim
00136 *>          LDB is INTEGER
00137 *>          The leading dimension of the matrix B. LDB >= max(1, N).
00138 *> \endverbatim
00139 *>
00140 *> \param[in,out] C
00141 *> \verbatim
00142 *>          C is COMPLEX array, dimension (LDC, N)
00143 *>          On entry, C contains the right-hand-side of the first matrix
00144 *>          equation in (1).
00145 *>          On exit, if IJOB = 0, C has been overwritten by the solution
00146 *>          R.
00147 *> \endverbatim
00148 *>
00149 *> \param[in] LDC
00150 *> \verbatim
00151 *>          LDC is INTEGER
00152 *>          The leading dimension of the matrix C. LDC >= max(1, M).
00153 *> \endverbatim
00154 *>
00155 *> \param[in] D
00156 *> \verbatim
00157 *>          D is COMPLEX array, dimension (LDD, M)
00158 *>          On entry, D contains an upper triangular matrix.
00159 *> \endverbatim
00160 *>
00161 *> \param[in] LDD
00162 *> \verbatim
00163 *>          LDD is INTEGER
00164 *>          The leading dimension of the matrix D. LDD >= max(1, M).
00165 *> \endverbatim
00166 *>
00167 *> \param[in] E
00168 *> \verbatim
00169 *>          E is COMPLEX array, dimension (LDE, N)
00170 *>          On entry, E contains an upper triangular matrix.
00171 *> \endverbatim
00172 *>
00173 *> \param[in] LDE
00174 *> \verbatim
00175 *>          LDE is INTEGER
00176 *>          The leading dimension of the matrix E. LDE >= max(1, N).
00177 *> \endverbatim
00178 *>
00179 *> \param[in,out] F
00180 *> \verbatim
00181 *>          F is COMPLEX array, dimension (LDF, N)
00182 *>          On entry, F contains the right-hand-side of the second matrix
00183 *>          equation in (1).
00184 *>          On exit, if IJOB = 0, F has been overwritten by the solution
00185 *>          L.
00186 *> \endverbatim
00187 *>
00188 *> \param[in] LDF
00189 *> \verbatim
00190 *>          LDF is INTEGER
00191 *>          The leading dimension of the matrix F. LDF >= max(1, M).
00192 *> \endverbatim
00193 *>
00194 *> \param[out] SCALE
00195 *> \verbatim
00196 *>          SCALE is REAL
00197 *>          On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions
00198 *>          R and L (C and F on entry) will hold the solutions to a
00199 *>          slightly perturbed system but the input matrices A, B, D and
00200 *>          E have not been changed. If SCALE = 0, R and L will hold the
00201 *>          solutions to the homogeneous system with C = F = 0.
00202 *>          Normally, SCALE = 1.
00203 *> \endverbatim
00204 *>
00205 *> \param[in,out] RDSUM
00206 *> \verbatim
00207 *>          RDSUM is REAL
00208 *>          On entry, the sum of squares of computed contributions to
00209 *>          the Dif-estimate under computation by CTGSYL, where the
00210 *>          scaling factor RDSCAL (see below) has been factored out.
00211 *>          On exit, the corresponding sum of squares updated with the
00212 *>          contributions from the current sub-system.
00213 *>          If TRANS = 'T' RDSUM is not touched.
00214 *>          NOTE: RDSUM only makes sense when CTGSY2 is called by
00215 *>          CTGSYL.
00216 *> \endverbatim
00217 *>
00218 *> \param[in,out] RDSCAL
00219 *> \verbatim
00220 *>          RDSCAL is REAL
00221 *>          On entry, scaling factor used to prevent overflow in RDSUM.
00222 *>          On exit, RDSCAL is updated w.r.t. the current contributions
00223 *>          in RDSUM.
00224 *>          If TRANS = 'T', RDSCAL is not touched.
00225 *>          NOTE: RDSCAL only makes sense when CTGSY2 is called by
00226 *>          CTGSYL.
00227 *> \endverbatim
00228 *>
00229 *> \param[out] INFO
00230 *> \verbatim
00231 *>          INFO is INTEGER
00232 *>          On exit, if INFO is set to
00233 *>            =0: Successful exit
00234 *>            <0: If INFO = -i, input argument number i is illegal.
00235 *>            >0: The matrix pairs (A, D) and (B, E) have common or very
00236 *>                close eigenvalues.
00237 *> \endverbatim
00238 *
00239 *  Authors:
00240 *  ========
00241 *
00242 *> \author Univ. of Tennessee 
00243 *> \author Univ. of California Berkeley 
00244 *> \author Univ. of Colorado Denver 
00245 *> \author NAG Ltd. 
00246 *
00247 *> \date November 2011
00248 *
00249 *> \ingroup complexSYauxiliary
00250 *
00251 *> \par Contributors:
00252 *  ==================
00253 *>
00254 *>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
00255 *>     Umea University, S-901 87 Umea, Sweden.
00256 *
00257 *  =====================================================================
00258       SUBROUTINE CTGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
00259      $                   LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL,
00260      $                   INFO )
00261 *
00262 *  -- LAPACK auxiliary routine (version 3.4.0) --
00263 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00264 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00265 *     November 2011
00266 *
00267 *     .. Scalar Arguments ..
00268       CHARACTER          TRANS
00269       INTEGER            IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N
00270       REAL               RDSCAL, RDSUM, SCALE
00271 *     ..
00272 *     .. Array Arguments ..
00273       COMPLEX            A( LDA, * ), B( LDB, * ), C( LDC, * ),
00274      $                   D( LDD, * ), E( LDE, * ), F( LDF, * )
00275 *     ..
00276 *
00277 *  =====================================================================
00278 *
00279 *     .. Parameters ..
00280       REAL               ZERO, ONE
00281       INTEGER            LDZ
00282       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, LDZ = 2 )
00283 *     ..
00284 *     .. Local Scalars ..
00285       LOGICAL            NOTRAN
00286       INTEGER            I, IERR, J, K
00287       REAL               SCALOC
00288       COMPLEX            ALPHA
00289 *     ..
00290 *     .. Local Arrays ..
00291       INTEGER            IPIV( LDZ ), JPIV( LDZ )
00292       COMPLEX            RHS( LDZ ), Z( LDZ, LDZ )
00293 *     ..
00294 *     .. External Functions ..
00295       LOGICAL            LSAME
00296       EXTERNAL           LSAME
00297 *     ..
00298 *     .. External Subroutines ..
00299       EXTERNAL           CAXPY, CGESC2, CGETC2, CSCAL, CLATDF, XERBLA
00300 *     ..
00301 *     .. Intrinsic Functions ..
00302       INTRINSIC          CMPLX, CONJG, MAX
00303 *     ..
00304 *     .. Executable Statements ..
00305 *
00306 *     Decode and test input parameters
00307 *
00308       INFO = 0
00309       IERR = 0
00310       NOTRAN = LSAME( TRANS, 'N' )
00311       IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
00312          INFO = -1
00313       ELSE IF( NOTRAN ) THEN
00314          IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.2 ) ) THEN
00315             INFO = -2
00316          END IF
00317       END IF
00318       IF( INFO.EQ.0 ) THEN
00319          IF( M.LE.0 ) THEN
00320             INFO = -3
00321          ELSE IF( N.LE.0 ) THEN
00322             INFO = -4
00323          ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00324             INFO = -5
00325          ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00326             INFO = -8
00327          ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
00328             INFO = -10
00329          ELSE IF( LDD.LT.MAX( 1, M ) ) THEN
00330             INFO = -12
00331          ELSE IF( LDE.LT.MAX( 1, N ) ) THEN
00332             INFO = -14
00333          ELSE IF( LDF.LT.MAX( 1, M ) ) THEN
00334             INFO = -16
00335          END IF
00336       END IF
00337       IF( INFO.NE.0 ) THEN
00338          CALL XERBLA( 'CTGSY2', -INFO )
00339          RETURN
00340       END IF
00341 *
00342       IF( NOTRAN ) THEN
00343 *
00344 *        Solve (I, J) - system
00345 *           A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
00346 *           D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
00347 *        for I = M, M - 1, ..., 1; J = 1, 2, ..., N
00348 *
00349          SCALE = ONE
00350          SCALOC = ONE
00351          DO 30 J = 1, N
00352             DO 20 I = M, 1, -1
00353 *
00354 *              Build 2 by 2 system
00355 *
00356                Z( 1, 1 ) = A( I, I )
00357                Z( 2, 1 ) = D( I, I )
00358                Z( 1, 2 ) = -B( J, J )
00359                Z( 2, 2 ) = -E( J, J )
00360 *
00361 *              Set up right hand side(s)
00362 *
00363                RHS( 1 ) = C( I, J )
00364                RHS( 2 ) = F( I, J )
00365 *
00366 *              Solve Z * x = RHS
00367 *
00368                CALL CGETC2( LDZ, Z, LDZ, IPIV, JPIV, IERR )
00369                IF( IERR.GT.0 )
00370      $            INFO = IERR
00371                IF( IJOB.EQ.0 ) THEN
00372                   CALL CGESC2( LDZ, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
00373                   IF( SCALOC.NE.ONE ) THEN
00374                      DO 10 K = 1, N
00375                         CALL CSCAL( M, CMPLX( SCALOC, ZERO ), C( 1, K ),
00376      $                              1 )
00377                         CALL CSCAL( M, CMPLX( SCALOC, ZERO ), F( 1, K ),
00378      $                              1 )
00379    10                CONTINUE
00380                      SCALE = SCALE*SCALOC
00381                   END IF
00382                ELSE
00383                   CALL CLATDF( IJOB, LDZ, Z, LDZ, RHS, RDSUM, RDSCAL,
00384      $                         IPIV, JPIV )
00385                END IF
00386 *
00387 *              Unpack solution vector(s)
00388 *
00389                C( I, J ) = RHS( 1 )
00390                F( I, J ) = RHS( 2 )
00391 *
00392 *              Substitute R(I, J) and L(I, J) into remaining equation.
00393 *
00394                IF( I.GT.1 ) THEN
00395                   ALPHA = -RHS( 1 )
00396                   CALL CAXPY( I-1, ALPHA, A( 1, I ), 1, C( 1, J ), 1 )
00397                   CALL CAXPY( I-1, ALPHA, D( 1, I ), 1, F( 1, J ), 1 )
00398                END IF
00399                IF( J.LT.N ) THEN
00400                   CALL CAXPY( N-J, RHS( 2 ), B( J, J+1 ), LDB,
00401      $                        C( I, J+1 ), LDC )
00402                   CALL CAXPY( N-J, RHS( 2 ), E( J, J+1 ), LDE,
00403      $                        F( I, J+1 ), LDF )
00404                END IF
00405 *
00406    20       CONTINUE
00407    30    CONTINUE
00408       ELSE
00409 *
00410 *        Solve transposed (I, J) - system:
00411 *           A(I, I)**H * R(I, J) + D(I, I)**H * L(J, J) = C(I, J)
00412 *           R(I, I) * B(J, J) + L(I, J) * E(J, J)   = -F(I, J)
00413 *        for I = 1, 2, ..., M, J = N, N - 1, ..., 1
00414 *
00415          SCALE = ONE
00416          SCALOC = ONE
00417          DO 80 I = 1, M
00418             DO 70 J = N, 1, -1
00419 *
00420 *              Build 2 by 2 system Z**H
00421 *
00422                Z( 1, 1 ) = CONJG( A( I, I ) )
00423                Z( 2, 1 ) = -CONJG( B( J, J ) )
00424                Z( 1, 2 ) = CONJG( D( I, I ) )
00425                Z( 2, 2 ) = -CONJG( E( J, J ) )
00426 *
00427 *
00428 *              Set up right hand side(s)
00429 *
00430                RHS( 1 ) = C( I, J )
00431                RHS( 2 ) = F( I, J )
00432 *
00433 *              Solve Z**H * x = RHS
00434 *
00435                CALL CGETC2( LDZ, Z, LDZ, IPIV, JPIV, IERR )
00436                IF( IERR.GT.0 )
00437      $            INFO = IERR
00438                CALL CGESC2( LDZ, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
00439                IF( SCALOC.NE.ONE ) THEN
00440                   DO 40 K = 1, N
00441                      CALL CSCAL( M, CMPLX( SCALOC, ZERO ), C( 1, K ),
00442      $                           1 )
00443                      CALL CSCAL( M, CMPLX( SCALOC, ZERO ), F( 1, K ),
00444      $                           1 )
00445    40             CONTINUE
00446                   SCALE = SCALE*SCALOC
00447                END IF
00448 *
00449 *              Unpack solution vector(s)
00450 *
00451                C( I, J ) = RHS( 1 )
00452                F( I, J ) = RHS( 2 )
00453 *
00454 *              Substitute R(I, J) and L(I, J) into remaining equation.
00455 *
00456                DO 50 K = 1, J - 1
00457                   F( I, K ) = F( I, K ) + RHS( 1 )*CONJG( B( K, J ) ) +
00458      $                        RHS( 2 )*CONJG( E( K, J ) )
00459    50          CONTINUE
00460                DO 60 K = I + 1, M
00461                   C( K, J ) = C( K, J ) - CONJG( A( I, K ) )*RHS( 1 ) -
00462      $                        CONJG( D( I, K ) )*RHS( 2 )
00463    60          CONTINUE
00464 *
00465    70       CONTINUE
00466    80    CONTINUE
00467       END IF
00468       RETURN
00469 *
00470 *     End of CTGSY2
00471 *
00472       END
 All Files Functions