LAPACK: SRC/dtgsy2.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DTGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
00002      $                   LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL,
00003      $                   IWORK, PQ, INFO )
00004 *
00005 *  -- LAPACK auxiliary routine (version 3.2) --
00006 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00007 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00008 *     January 2007
00009 *
00010 *     .. Scalar Arguments ..
00011       CHARACTER          TRANS
00012       INTEGER            IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N,
00013      $                   PQ
00014       DOUBLE PRECISION   RDSCAL, RDSUM, SCALE
00015 *     ..
00016 *     .. Array Arguments ..
00017       INTEGER            IWORK( * )
00018       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( LDC, * ),
00019      $                   D( LDD, * ), E( LDE, * ), F( LDF, * )
00020 *     ..
00021 *
00022 *  Purpose
00023 *  =======
00024 *
00025 *  DTGSY2 solves the generalized Sylvester equation:
00026 *
00027 *              A * R - L * B = scale * C                (1)
00028 *              D * R - L * E = scale * F,
00029 *
00030 *  using Level 1 and 2 BLAS. where R and L are unknown M-by-N matrices,
00031 *  (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M,
00032 *  N-by-N and M-by-N, respectively, with real entries. (A, D) and (B, E)
00033 *  must be in generalized Schur canonical form, i.e. A, B are upper
00034 *  quasi triangular and D, E are upper triangular. The solution (R, L)
00035 *  overwrites (C, F). 0 <= SCALE <= 1 is an output scaling factor
00036 *  chosen to avoid overflow.
00037 *
00038 *  In matrix notation solving equation (1) corresponds to solve
00039 *  Z*x = scale*b, where Z is defined as
00040 *
00041 *         Z = [ kron(In, A)  -kron(B', Im) ]             (2)
00042 *             [ kron(In, D)  -kron(E', Im) ],
00043 *
00044 *  Ik is the identity matrix of size k and X' is the transpose of X.
00045 *  kron(X, Y) is the Kronecker product between the matrices X and Y.
00046 *  In the process of solving (1), we solve a number of such systems
00047 *  where Dim(In), Dim(In) = 1 or 2.
00048 *
00049 *  If TRANS = 'T', solve the transposed system Z'*y = scale*b for y,
00050 *  which is equivalent to solve for R and L in
00051 *
00052 *              A' * R  + D' * L   = scale *  C           (3)
00053 *              R  * B' + L  * E'  = scale * -F
00054 *
00055 *  This case is used to compute an estimate of Dif[(A, D), (B, E)] =
00056 *  sigma_min(Z) using reverse communicaton with DLACON.
00057 *
00058 *  DTGSY2 also (IJOB >= 1) contributes to the computation in DTGSYL
00059 *  of an upper bound on the separation between to matrix pairs. Then
00060 *  the input (A, D), (B, E) are sub-pencils of the matrix pair in
00061 *  DTGSYL. See DTGSYL for details.
00062 *
00063 *  Arguments
00064 *  =========
00065 *
00066 *  TRANS   (input) CHARACTER*1
00067 *          = 'N', solve the generalized Sylvester equation (1).
00068 *          = 'T': solve the 'transposed' system (3).
00069 *
00070 *  IJOB    (input) INTEGER
00071 *          Specifies what kind of functionality to be performed.
00072 *          = 0: solve (1) only.
00073 *          = 1: A contribution from this subsystem to a Frobenius
00074 *               norm-based estimate of the separation between two matrix
00075 *               pairs is computed. (look ahead strategy is used).
00076 *          = 2: A contribution from this subsystem to a Frobenius
00077 *               norm-based estimate of the separation between two matrix
00078 *               pairs is computed. (DGECON on sub-systems is used.)
00079 *          Not referenced if TRANS = 'T'.
00080 *
00081 *  M       (input) INTEGER
00082 *          On entry, M specifies the order of A and D, and the row
00083 *          dimension of C, F, R and L.
00084 *
00085 *  N       (input) INTEGER
00086 *          On entry, N specifies the order of B and E, and the column
00087 *          dimension of C, F, R and L.
00088 *
00089 *  A       (input) DOUBLE PRECISION array, dimension (LDA, M)
00090 *          On entry, A contains an upper quasi triangular matrix.
00091 *
00092 *  LDA     (input) INTEGER
00093 *          The leading dimension of the matrix A. LDA >= max(1, M).
00094 *
00095 *  B       (input) DOUBLE PRECISION array, dimension (LDB, N)
00096 *          On entry, B contains an upper quasi triangular matrix.
00097 *
00098 *  LDB     (input) INTEGER
00099 *          The leading dimension of the matrix B. LDB >= max(1, N).
00100 *
00101 *  C       (input/output) DOUBLE PRECISION array, dimension (LDC, N)
00102 *          On entry, C contains the right-hand-side of the first matrix
00103 *          equation in (1).
00104 *          On exit, if IJOB = 0, C has been overwritten by the
00105 *          solution R.
00106 *
00107 *  LDC     (input) INTEGER
00108 *          The leading dimension of the matrix C. LDC >= max(1, M).
00109 *
00110 *  D       (input) DOUBLE PRECISION array, dimension (LDD, M)
00111 *          On entry, D contains an upper triangular matrix.
00112 *
00113 *  LDD     (input) INTEGER
00114 *          The leading dimension of the matrix D. LDD >= max(1, M).
00115 *
00116 *  E       (input) DOUBLE PRECISION array, dimension (LDE, N)
00117 *          On entry, E contains an upper triangular matrix.
00118 *
00119 *  LDE     (input) INTEGER
00120 *          The leading dimension of the matrix E. LDE >= max(1, N).
00121 *
00122 *  F       (input/output) DOUBLE PRECISION array, dimension (LDF, N)
00123 *          On entry, F contains the right-hand-side of the second matrix
00124 *          equation in (1).
00125 *          On exit, if IJOB = 0, F has been overwritten by the
00126 *          solution L.
00127 *
00128 *  LDF     (input) INTEGER
00129 *          The leading dimension of the matrix F. LDF >= max(1, M).
00130 *
00131 *  SCALE   (output) DOUBLE PRECISION
00132 *          On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions
00133 *          R and L (C and F on entry) will hold the solutions to a
00134 *          slightly perturbed system but the input matrices A, B, D and
00135 *          E have not been changed. If SCALE = 0, R and L will hold the
00136 *          solutions to the homogeneous system with C = F = 0. Normally,
00137 *          SCALE = 1.
00138 *
00139 *  RDSUM   (input/output) DOUBLE PRECISION
00140 *          On entry, the sum of squares of computed contributions to
00141 *          the Dif-estimate under computation by DTGSYL, where the
00142 *          scaling factor RDSCAL (see below) has been factored out.
00143 *          On exit, the corresponding sum of squares updated with the
00144 *          contributions from the current sub-system.
00145 *          If TRANS = 'T' RDSUM is not touched.
00146 *          NOTE: RDSUM only makes sense when DTGSY2 is called by DTGSYL.
00147 *
00148 *  RDSCAL  (input/output) DOUBLE PRECISION
00149 *          On entry, scaling factor used to prevent overflow in RDSUM.
00150 *          On exit, RDSCAL is updated w.r.t. the current contributions
00151 *          in RDSUM.
00152 *          If TRANS = 'T', RDSCAL is not touched.
00153 *          NOTE: RDSCAL only makes sense when DTGSY2 is called by
00154 *                DTGSYL.
00155 *
00156 *  IWORK   (workspace) INTEGER array, dimension (M+N+2)
00157 *
00158 *  PQ      (output) INTEGER
00159 *          On exit, the number of subsystems (of size 2-by-2, 4-by-4 and
00160 *          8-by-8) solved by this routine.
00161 *
00162 *  INFO    (output) INTEGER
00163 *          On exit, if INFO is set to
00164 *            =0: Successful exit
00165 *            <0: If INFO = -i, the i-th argument had an illegal value.
00166 *            >0: The matrix pairs (A, D) and (B, E) have common or very
00167 *                close eigenvalues.
00168 *
00169 *  Further Details
00170 *  ===============
00171 *
00172 *  Based on contributions by
00173 *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
00174 *     Umea University, S-901 87 Umea, Sweden.
00175 *
00176 *  =====================================================================
00177 *  Replaced various illegal calls to DCOPY by calls to DLASET.
00178 *  Sven Hammarling, 27/5/02.
00179 *
00180 *     .. Parameters ..
00181       INTEGER            LDZ
00182       PARAMETER          ( LDZ = 8 )
00183       DOUBLE PRECISION   ZERO, ONE
00184       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00185 *     ..
00186 *     .. Local Scalars ..
00187       LOGICAL            NOTRAN
00188       INTEGER            I, IE, IERR, II, IS, ISP1, J, JE, JJ, JS, JSP1,
00189      $                   K, MB, NB, P, Q, ZDIM
00190       DOUBLE PRECISION   ALPHA, SCALOC
00191 *     ..
00192 *     .. Local Arrays ..
00193       INTEGER            IPIV( LDZ ), JPIV( LDZ )
00194       DOUBLE PRECISION   RHS( LDZ ), Z( LDZ, LDZ )
00195 *     ..
00196 *     .. External Functions ..
00197       LOGICAL            LSAME
00198       EXTERNAL           LSAME
00199 *     ..
00200 *     .. External Subroutines ..
00201       EXTERNAL           DAXPY, DCOPY, DGEMM, DGEMV, DGER, DGESC2,
00202      $                   DGETC2, DLASET, DLATDF, DSCAL, XERBLA
00203 *     ..
00204 *     .. Intrinsic Functions ..
00205       INTRINSIC          MAX
00206 *     ..
00207 *     .. Executable Statements ..
00208 *
00209 *     Decode and test input parameters
00210 *
00211       INFO = 0
00212       IERR = 0
00213       NOTRAN = LSAME( TRANS, 'N' )
00214       IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
00215          INFO = -1
00216       ELSE IF( NOTRAN ) THEN
00217          IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.2 ) ) THEN
00218             INFO = -2
00219          END IF
00220       END IF
00221       IF( INFO.EQ.0 ) THEN
00222          IF( M.LE.0 ) THEN
00223             INFO = -3
00224          ELSE IF( N.LE.0 ) THEN
00225             INFO = -4
00226          ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00227             INFO = -5
00228          ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00229             INFO = -8
00230          ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
00231             INFO = -10
00232          ELSE IF( LDD.LT.MAX( 1, M ) ) THEN
00233             INFO = -12
00234          ELSE IF( LDE.LT.MAX( 1, N ) ) THEN
00235             INFO = -14
00236          ELSE IF( LDF.LT.MAX( 1, M ) ) THEN
00237             INFO = -16
00238          END IF
00239       END IF
00240       IF( INFO.NE.0 ) THEN
00241          CALL XERBLA( 'DTGSY2', -INFO )
00242          RETURN
00243       END IF
00244 *
00245 *     Determine block structure of A
00246 *
00247       PQ = 0
00248       P = 0
00249       I = 1
00250    10 CONTINUE
00251       IF( I.GT.M )
00252      $   GO TO 20
00253       P = P + 1
00254       IWORK( P ) = I
00255       IF( I.EQ.M )
00256      $   GO TO 20
00257       IF( A( I+1, I ).NE.ZERO ) THEN
00258          I = I + 2
00259       ELSE
00260          I = I + 1
00261       END IF
00262       GO TO 10
00263    20 CONTINUE
00264       IWORK( P+1 ) = M + 1
00265 *
00266 *     Determine block structure of B
00267 *
00268       Q = P + 1
00269       J = 1
00270    30 CONTINUE
00271       IF( J.GT.N )
00272      $   GO TO 40
00273       Q = Q + 1
00274       IWORK( Q ) = J
00275       IF( J.EQ.N )
00276      $   GO TO 40
00277       IF( B( J+1, J ).NE.ZERO ) THEN
00278          J = J + 2
00279       ELSE
00280          J = J + 1
00281       END IF
00282       GO TO 30
00283    40 CONTINUE
00284       IWORK( Q+1 ) = N + 1
00285       PQ = P*( Q-P-1 )
00286 *
00287       IF( NOTRAN ) THEN
00288 *
00289 *        Solve (I, J) - subsystem
00290 *           A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
00291 *           D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
00292 *        for I = P, P - 1, ..., 1; J = 1, 2, ..., Q
00293 *
00294          SCALE = ONE
00295          SCALOC = ONE
00296          DO 120 J = P + 2, Q
00297             JS = IWORK( J )
00298             JSP1 = JS + 1
00299             JE = IWORK( J+1 ) - 1
00300             NB = JE - JS + 1
00301             DO 110 I = P, 1, -1
00302 *
00303                IS = IWORK( I )
00304                ISP1 = IS + 1
00305                IE = IWORK( I+1 ) - 1
00306                MB = IE - IS + 1
00307                ZDIM = MB*NB*2
00308 *
00309                IF( ( MB.EQ.1 ) .AND. ( NB.EQ.1 ) ) THEN
00310 *
00311 *                 Build a 2-by-2 system Z * x = RHS
00312 *
00313                   Z( 1, 1 ) = A( IS, IS )
00314                   Z( 2, 1 ) = D( IS, IS )
00315                   Z( 1, 2 ) = -B( JS, JS )
00316                   Z( 2, 2 ) = -E( JS, JS )
00317 *
00318 *                 Set up right hand side(s)
00319 *
00320                   RHS( 1 ) = C( IS, JS )
00321                   RHS( 2 ) = F( IS, JS )
00322 *
00323 *                 Solve Z * x = RHS
00324 *
00325                   CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
00326                   IF( IERR.GT.0 )
00327      $               INFO = IERR
00328 *
00329                   IF( IJOB.EQ.0 ) THEN
00330                      CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
00331      $                            SCALOC )
00332                      IF( SCALOC.NE.ONE ) THEN
00333                         DO 50 K = 1, N
00334                            CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
00335                            CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
00336    50                   CONTINUE
00337                         SCALE = SCALE*SCALOC
00338                      END IF
00339                   ELSE
00340                      CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM,
00341      $                            RDSCAL, IPIV, JPIV )
00342                   END IF
00343 *
00344 *                 Unpack solution vector(s)
00345 *
00346                   C( IS, JS ) = RHS( 1 )
00347                   F( IS, JS ) = RHS( 2 )
00348 *
00349 *                 Substitute R(I, J) and L(I, J) into remaining
00350 *                 equation.
00351 *
00352                   IF( I.GT.1 ) THEN
00353                      ALPHA = -RHS( 1 )
00354                      CALL DAXPY( IS-1, ALPHA, A( 1, IS ), 1, C( 1, JS ),
00355      $                           1 )
00356                      CALL DAXPY( IS-1, ALPHA, D( 1, IS ), 1, F( 1, JS ),
00357      $                           1 )
00358                   END IF
00359                   IF( J.LT.Q ) THEN
00360                      CALL DAXPY( N-JE, RHS( 2 ), B( JS, JE+1 ), LDB,
00361      $                           C( IS, JE+1 ), LDC )
00362                      CALL DAXPY( N-JE, RHS( 2 ), E( JS, JE+1 ), LDE,
00363      $                           F( IS, JE+1 ), LDF )
00364                   END IF
00365 *
00366                ELSE IF( ( MB.EQ.1 ) .AND. ( NB.EQ.2 ) ) THEN
00367 *
00368 *                 Build a 4-by-4 system Z * x = RHS
00369 *
00370                   Z( 1, 1 ) = A( IS, IS )
00371                   Z( 2, 1 ) = ZERO
00372                   Z( 3, 1 ) = D( IS, IS )
00373                   Z( 4, 1 ) = ZERO
00374 *
00375                   Z( 1, 2 ) = ZERO
00376                   Z( 2, 2 ) = A( IS, IS )
00377                   Z( 3, 2 ) = ZERO
00378                   Z( 4, 2 ) = D( IS, IS )
00379 *
00380                   Z( 1, 3 ) = -B( JS, JS )
00381                   Z( 2, 3 ) = -B( JS, JSP1 )
00382                   Z( 3, 3 ) = -E( JS, JS )
00383                   Z( 4, 3 ) = -E( JS, JSP1 )
00384 *
00385                   Z( 1, 4 ) = -B( JSP1, JS )
00386                   Z( 2, 4 ) = -B( JSP1, JSP1 )
00387                   Z( 3, 4 ) = ZERO
00388                   Z( 4, 4 ) = -E( JSP1, JSP1 )
00389 *
00390 *                 Set up right hand side(s)
00391 *
00392                   RHS( 1 ) = C( IS, JS )
00393                   RHS( 2 ) = C( IS, JSP1 )
00394                   RHS( 3 ) = F( IS, JS )
00395                   RHS( 4 ) = F( IS, JSP1 )
00396 *
00397 *                 Solve Z * x = RHS
00398 *
00399                   CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
00400                   IF( IERR.GT.0 )
00401      $               INFO = IERR
00402 *
00403                   IF( IJOB.EQ.0 ) THEN
00404                      CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
00405      $                            SCALOC )
00406                      IF( SCALOC.NE.ONE ) THEN
00407                         DO 60 K = 1, N
00408                            CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
00409                            CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
00410    60                   CONTINUE
00411                         SCALE = SCALE*SCALOC
00412                      END IF
00413                   ELSE
00414                      CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM,
00415      $                            RDSCAL, IPIV, JPIV )
00416                   END IF
00417 *
00418 *                 Unpack solution vector(s)
00419 *
00420                   C( IS, JS ) = RHS( 1 )
00421                   C( IS, JSP1 ) = RHS( 2 )
00422                   F( IS, JS ) = RHS( 3 )
00423                   F( IS, JSP1 ) = RHS( 4 )
00424 *
00425 *                 Substitute R(I, J) and L(I, J) into remaining
00426 *                 equation.
00427 *
00428                   IF( I.GT.1 ) THEN
00429                      CALL DGER( IS-1, NB, -ONE, A( 1, IS ), 1, RHS( 1 ),
00430      $                          1, C( 1, JS ), LDC )
00431                      CALL DGER( IS-1, NB, -ONE, D( 1, IS ), 1, RHS( 1 ),
00432      $                          1, F( 1, JS ), LDF )
00433                   END IF
00434                   IF( J.LT.Q ) THEN
00435                      CALL DAXPY( N-JE, RHS( 3 ), B( JS, JE+1 ), LDB,
00436      $                           C( IS, JE+1 ), LDC )
00437                      CALL DAXPY( N-JE, RHS( 3 ), E( JS, JE+1 ), LDE,
00438      $                           F( IS, JE+1 ), LDF )
00439                      CALL DAXPY( N-JE, RHS( 4 ), B( JSP1, JE+1 ), LDB,
00440      $                           C( IS, JE+1 ), LDC )
00441                      CALL DAXPY( N-JE, RHS( 4 ), E( JSP1, JE+1 ), LDE,
00442      $                           F( IS, JE+1 ), LDF )
00443                   END IF
00444 *
00445                ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.1 ) ) THEN
00446 *
00447 *                 Build a 4-by-4 system Z * x = RHS
00448 *
00449                   Z( 1, 1 ) = A( IS, IS )
00450                   Z( 2, 1 ) = A( ISP1, IS )
00451                   Z( 3, 1 ) = D( IS, IS )
00452                   Z( 4, 1 ) = ZERO
00453 *
00454                   Z( 1, 2 ) = A( IS, ISP1 )
00455                   Z( 2, 2 ) = A( ISP1, ISP1 )
00456                   Z( 3, 2 ) = D( IS, ISP1 )
00457                   Z( 4, 2 ) = D( ISP1, ISP1 )
00458 *
00459                   Z( 1, 3 ) = -B( JS, JS )
00460                   Z( 2, 3 ) = ZERO
00461                   Z( 3, 3 ) = -E( JS, JS )
00462                   Z( 4, 3 ) = ZERO
00463 *
00464                   Z( 1, 4 ) = ZERO
00465                   Z( 2, 4 ) = -B( JS, JS )
00466                   Z( 3, 4 ) = ZERO
00467                   Z( 4, 4 ) = -E( JS, JS )
00468 *
00469 *                 Set up right hand side(s)
00470 *
00471                   RHS( 1 ) = C( IS, JS )
00472                   RHS( 2 ) = C( ISP1, JS )
00473                   RHS( 3 ) = F( IS, JS )
00474                   RHS( 4 ) = F( ISP1, JS )
00475 *
00476 *                 Solve Z * x = RHS
00477 *
00478                   CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
00479                   IF( IERR.GT.0 )
00480      $               INFO = IERR
00481                   IF( IJOB.EQ.0 ) THEN
00482                      CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
00483      $                            SCALOC )
00484                      IF( SCALOC.NE.ONE ) THEN
00485                         DO 70 K = 1, N
00486                            CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
00487                            CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
00488    70                   CONTINUE
00489                         SCALE = SCALE*SCALOC
00490                      END IF
00491                   ELSE
00492                      CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM,
00493      $                            RDSCAL, IPIV, JPIV )
00494                   END IF
00495 *
00496 *                 Unpack solution vector(s)
00497 *
00498                   C( IS, JS ) = RHS( 1 )
00499                   C( ISP1, JS ) = RHS( 2 )
00500                   F( IS, JS ) = RHS( 3 )
00501                   F( ISP1, JS ) = RHS( 4 )
00502 *
00503 *                 Substitute R(I, J) and L(I, J) into remaining
00504 *                 equation.
00505 *
00506                   IF( I.GT.1 ) THEN
00507                      CALL DGEMV( 'N', IS-1, MB, -ONE, A( 1, IS ), LDA,
00508      $                           RHS( 1 ), 1, ONE, C( 1, JS ), 1 )
00509                      CALL DGEMV( 'N', IS-1, MB, -ONE, D( 1, IS ), LDD,
00510      $                           RHS( 1 ), 1, ONE, F( 1, JS ), 1 )
00511                   END IF
00512                   IF( J.LT.Q ) THEN
00513                      CALL DGER( MB, N-JE, ONE, RHS( 3 ), 1,
00514      $                          B( JS, JE+1 ), LDB, C( IS, JE+1 ), LDC )
00515                      CALL DGER( MB, N-JE, ONE, RHS( 3 ), 1,
00516      $                          E( JS, JE+1 ), LDE, F( IS, JE+1 ), LDF )
00517                   END IF
00518 *
00519                ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.2 ) ) THEN
00520 *
00521 *                 Build an 8-by-8 system Z * x = RHS
00522 *
00523                   CALL DLASET( 'F', LDZ, LDZ, ZERO, ZERO, Z, LDZ )
00524 *
00525                   Z( 1, 1 ) = A( IS, IS )
00526                   Z( 2, 1 ) = A( ISP1, IS )
00527                   Z( 5, 1 ) = D( IS, IS )
00528 *
00529                   Z( 1, 2 ) = A( IS, ISP1 )
00530                   Z( 2, 2 ) = A( ISP1, ISP1 )
00531                   Z( 5, 2 ) = D( IS, ISP1 )
00532                   Z( 6, 2 ) = D( ISP1, ISP1 )
00533 *
00534                   Z( 3, 3 ) = A( IS, IS )
00535                   Z( 4, 3 ) = A( ISP1, IS )
00536                   Z( 7, 3 ) = D( IS, IS )
00537 *
00538                   Z( 3, 4 ) = A( IS, ISP1 )
00539                   Z( 4, 4 ) = A( ISP1, ISP1 )
00540                   Z( 7, 4 ) = D( IS, ISP1 )
00541                   Z( 8, 4 ) = D( ISP1, ISP1 )
00542 *
00543                   Z( 1, 5 ) = -B( JS, JS )
00544                   Z( 3, 5 ) = -B( JS, JSP1 )
00545                   Z( 5, 5 ) = -E( JS, JS )
00546                   Z( 7, 5 ) = -E( JS, JSP1 )
00547 *
00548                   Z( 2, 6 ) = -B( JS, JS )
00549                   Z( 4, 6 ) = -B( JS, JSP1 )
00550                   Z( 6, 6 ) = -E( JS, JS )
00551                   Z( 8, 6 ) = -E( JS, JSP1 )
00552 *
00553                   Z( 1, 7 ) = -B( JSP1, JS )
00554                   Z( 3, 7 ) = -B( JSP1, JSP1 )
00555                   Z( 7, 7 ) = -E( JSP1, JSP1 )
00556 *
00557                   Z( 2, 8 ) = -B( JSP1, JS )
00558                   Z( 4, 8 ) = -B( JSP1, JSP1 )
00559                   Z( 8, 8 ) = -E( JSP1, JSP1 )
00560 *
00561 *                 Set up right hand side(s)
00562 *
00563                   K = 1
00564                   II = MB*NB + 1
00565                   DO 80 JJ = 0, NB - 1
00566                      CALL DCOPY( MB, C( IS, JS+JJ ), 1, RHS( K ), 1 )
00567                      CALL DCOPY( MB, F( IS, JS+JJ ), 1, RHS( II ), 1 )
00568                      K = K + MB
00569                      II = II + MB
00570    80             CONTINUE
00571 *
00572 *                 Solve Z * x = RHS
00573 *
00574                   CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
00575                   IF( IERR.GT.0 )
00576      $               INFO = IERR
00577                   IF( IJOB.EQ.0 ) THEN
00578                      CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
00579      $                            SCALOC )
00580                      IF( SCALOC.NE.ONE ) THEN
00581                         DO 90 K = 1, N
00582                            CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
00583                            CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
00584    90                   CONTINUE
00585                         SCALE = SCALE*SCALOC
00586                      END IF
00587                   ELSE
00588                      CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM,
00589      $                            RDSCAL, IPIV, JPIV )
00590                   END IF
00591 *
00592 *                 Unpack solution vector(s)
00593 *
00594                   K = 1
00595                   II = MB*NB + 1
00596                   DO 100 JJ = 0, NB - 1
00597                      CALL DCOPY( MB, RHS( K ), 1, C( IS, JS+JJ ), 1 )
00598                      CALL DCOPY( MB, RHS( II ), 1, F( IS, JS+JJ ), 1 )
00599                      K = K + MB
00600                      II = II + MB
00601   100             CONTINUE
00602 *
00603 *                 Substitute R(I, J) and L(I, J) into remaining
00604 *                 equation.
00605 *
00606                   IF( I.GT.1 ) THEN
00607                      CALL DGEMM( 'N', 'N', IS-1, NB, MB, -ONE,
00608      $                           A( 1, IS ), LDA, RHS( 1 ), MB, ONE,
00609      $                           C( 1, JS ), LDC )
00610                      CALL DGEMM( 'N', 'N', IS-1, NB, MB, -ONE,
00611      $                           D( 1, IS ), LDD, RHS( 1 ), MB, ONE,
00612      $                           F( 1, JS ), LDF )
00613                   END IF
00614                   IF( J.LT.Q ) THEN
00615                      K = MB*NB + 1
00616                      CALL DGEMM( 'N', 'N', MB, N-JE, NB, ONE, RHS( K ),
00617      $                           MB, B( JS, JE+1 ), LDB, ONE,
00618      $                           C( IS, JE+1 ), LDC )
00619                      CALL DGEMM( 'N', 'N', MB, N-JE, NB, ONE, RHS( K ),
00620      $                           MB, E( JS, JE+1 ), LDE, ONE,
00621      $                           F( IS, JE+1 ), LDF )
00622                   END IF
00623 *
00624                END IF
00625 *
00626   110       CONTINUE
00627   120    CONTINUE
00628       ELSE
00629 *
00630 *        Solve (I, J) - subsystem
00631 *             A(I, I)' * R(I, J) + D(I, I)' * L(J, J)  =  C(I, J)
00632 *             R(I, I)  * B(J, J) + L(I, J)  * E(J, J)  = -F(I, J)
00633 *        for I = 1, 2, ..., P, J = Q, Q - 1, ..., 1
00634 *
00635          SCALE = ONE
00636          SCALOC = ONE
00637          DO 200 I = 1, P
00638 *
00639             IS = IWORK( I )
00640             ISP1 = IS + 1
00641             IE = IWORK ( I+1 ) - 1
00642             MB = IE - IS + 1
00643             DO 190 J = Q, P + 2, -1
00644 *
00645                JS = IWORK( J )
00646                JSP1 = JS + 1
00647                JE = IWORK( J+1 ) - 1
00648                NB = JE - JS + 1
00649                ZDIM = MB*NB*2
00650                IF( ( MB.EQ.1 ) .AND. ( NB.EQ.1 ) ) THEN
00651 *
00652 *                 Build a 2-by-2 system Z' * x = RHS
00653 *
00654                   Z( 1, 1 ) = A( IS, IS )
00655                   Z( 2, 1 ) = -B( JS, JS )
00656                   Z( 1, 2 ) = D( IS, IS )
00657                   Z( 2, 2 ) = -E( JS, JS )
00658 *
00659 *                 Set up right hand side(s)
00660 *
00661                   RHS( 1 ) = C( IS, JS )
00662                   RHS( 2 ) = F( IS, JS )
00663 *
00664 *                 Solve Z' * x = RHS
00665 *
00666                   CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
00667                   IF( IERR.GT.0 )
00668      $               INFO = IERR
00669 *
00670                   CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
00671                   IF( SCALOC.NE.ONE ) THEN
00672                      DO 130 K = 1, N
00673                         CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
00674                         CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
00675   130                CONTINUE
00676                      SCALE = SCALE*SCALOC
00677                   END IF
00678 *
00679 *                 Unpack solution vector(s)
00680 *
00681                   C( IS, JS ) = RHS( 1 )
00682                   F( IS, JS ) = RHS( 2 )
00683 *
00684 *                 Substitute R(I, J) and L(I, J) into remaining
00685 *                 equation.
00686 *
00687                   IF( J.GT.P+2 ) THEN
00688                      ALPHA = RHS( 1 )
00689                      CALL DAXPY( JS-1, ALPHA, B( 1, JS ), 1, F( IS, 1 ),
00690      $                           LDF )
00691                      ALPHA = RHS( 2 )
00692                      CALL DAXPY( JS-1, ALPHA, E( 1, JS ), 1, F( IS, 1 ),
00693      $                           LDF )
00694                   END IF
00695                   IF( I.LT.P ) THEN
00696                      ALPHA = -RHS( 1 )
00697                      CALL DAXPY( M-IE, ALPHA, A( IS, IE+1 ), LDA,
00698      $                           C( IE+1, JS ), 1 )
00699                      ALPHA = -RHS( 2 )
00700                      CALL DAXPY( M-IE, ALPHA, D( IS, IE+1 ), LDD,
00701      $                           C( IE+1, JS ), 1 )
00702                   END IF
00703 *
00704                ELSE IF( ( MB.EQ.1 ) .AND. ( NB.EQ.2 ) ) THEN
00705 *
00706 *                 Build a 4-by-4 system Z' * x = RHS
00707 *
00708                   Z( 1, 1 ) = A( IS, IS )
00709                   Z( 2, 1 ) = ZERO
00710                   Z( 3, 1 ) = -B( JS, JS )
00711                   Z( 4, 1 ) = -B( JSP1, JS )
00712 *
00713                   Z( 1, 2 ) = ZERO
00714                   Z( 2, 2 ) = A( IS, IS )
00715                   Z( 3, 2 ) = -B( JS, JSP1 )
00716                   Z( 4, 2 ) = -B( JSP1, JSP1 )
00717 *
00718                   Z( 1, 3 ) = D( IS, IS )
00719                   Z( 2, 3 ) = ZERO
00720                   Z( 3, 3 ) = -E( JS, JS )
00721                   Z( 4, 3 ) = ZERO
00722 *
00723                   Z( 1, 4 ) = ZERO
00724                   Z( 2, 4 ) = D( IS, IS )
00725                   Z( 3, 4 ) = -E( JS, JSP1 )
00726                   Z( 4, 4 ) = -E( JSP1, JSP1 )
00727 *
00728 *                 Set up right hand side(s)
00729 *
00730                   RHS( 1 ) = C( IS, JS )
00731                   RHS( 2 ) = C( IS, JSP1 )
00732                   RHS( 3 ) = F( IS, JS )
00733                   RHS( 4 ) = F( IS, JSP1 )
00734 *
00735 *                 Solve Z' * x = RHS
00736 *
00737                   CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
00738                   IF( IERR.GT.0 )
00739      $               INFO = IERR
00740                   CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
00741                   IF( SCALOC.NE.ONE ) THEN
00742                      DO 140 K = 1, N
00743                         CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
00744                         CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
00745   140                CONTINUE
00746                      SCALE = SCALE*SCALOC
00747                   END IF
00748 *
00749 *                 Unpack solution vector(s)
00750 *
00751                   C( IS, JS ) = RHS( 1 )
00752                   C( IS, JSP1 ) = RHS( 2 )
00753                   F( IS, JS ) = RHS( 3 )
00754                   F( IS, JSP1 ) = RHS( 4 )
00755 *
00756 *                 Substitute R(I, J) and L(I, J) into remaining
00757 *                 equation.
00758 *
00759                   IF( J.GT.P+2 ) THEN
00760                      CALL DAXPY( JS-1, RHS( 1 ), B( 1, JS ), 1,
00761      $                           F( IS, 1 ), LDF )
00762                      CALL DAXPY( JS-1, RHS( 2 ), B( 1, JSP1 ), 1,
00763      $                           F( IS, 1 ), LDF )
00764                      CALL DAXPY( JS-1, RHS( 3 ), E( 1, JS ), 1,
00765      $                           F( IS, 1 ), LDF )
00766                      CALL DAXPY( JS-1, RHS( 4 ), E( 1, JSP1 ), 1,
00767      $                           F( IS, 1 ), LDF )
00768                   END IF
00769                   IF( I.LT.P ) THEN
00770                      CALL DGER( M-IE, NB, -ONE, A( IS, IE+1 ), LDA,
00771      $                          RHS( 1 ), 1, C( IE+1, JS ), LDC )
00772                      CALL DGER( M-IE, NB, -ONE, D( IS, IE+1 ), LDD,
00773      $                          RHS( 3 ), 1, C( IE+1, JS ), LDC )
00774                   END IF
00775 *
00776                ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.1 ) ) THEN
00777 *
00778 *                 Build a 4-by-4 system Z' * x = RHS
00779 *
00780                   Z( 1, 1 ) = A( IS, IS )
00781                   Z( 2, 1 ) = A( IS, ISP1 )
00782                   Z( 3, 1 ) = -B( JS, JS )
00783                   Z( 4, 1 ) = ZERO
00784 *
00785                   Z( 1, 2 ) = A( ISP1, IS )
00786                   Z( 2, 2 ) = A( ISP1, ISP1 )
00787                   Z( 3, 2 ) = ZERO
00788                   Z( 4, 2 ) = -B( JS, JS )
00789 *
00790                   Z( 1, 3 ) = D( IS, IS )
00791                   Z( 2, 3 ) = D( IS, ISP1 )
00792                   Z( 3, 3 ) = -E( JS, JS )
00793                   Z( 4, 3 ) = ZERO
00794 *
00795                   Z( 1, 4 ) = ZERO
00796                   Z( 2, 4 ) = D( ISP1, ISP1 )
00797                   Z( 3, 4 ) = ZERO
00798                   Z( 4, 4 ) = -E( JS, JS )
00799 *
00800 *                 Set up right hand side(s)
00801 *
00802                   RHS( 1 ) = C( IS, JS )
00803                   RHS( 2 ) = C( ISP1, JS )
00804                   RHS( 3 ) = F( IS, JS )
00805                   RHS( 4 ) = F( ISP1, JS )
00806 *
00807 *                 Solve Z' * x = RHS
00808 *
00809                   CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
00810                   IF( IERR.GT.0 )
00811      $               INFO = IERR
00812 *
00813                   CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
00814                   IF( SCALOC.NE.ONE ) THEN
00815                      DO 150 K = 1, N
00816                         CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
00817                         CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
00818   150                CONTINUE
00819                      SCALE = SCALE*SCALOC
00820                   END IF
00821 *
00822 *                 Unpack solution vector(s)
00823 *
00824                   C( IS, JS ) = RHS( 1 )
00825                   C( ISP1, JS ) = RHS( 2 )
00826                   F( IS, JS ) = RHS( 3 )
00827                   F( ISP1, JS ) = RHS( 4 )
00828 *
00829 *                 Substitute R(I, J) and L(I, J) into remaining
00830 *                 equation.
00831 *
00832                   IF( J.GT.P+2 ) THEN
00833                      CALL DGER( MB, JS-1, ONE, RHS( 1 ), 1, B( 1, JS ),
00834      $                          1, F( IS, 1 ), LDF )
00835                      CALL DGER( MB, JS-1, ONE, RHS( 3 ), 1, E( 1, JS ),
00836      $                          1, F( IS, 1 ), LDF )
00837                   END IF
00838                   IF( I.LT.P ) THEN
00839                      CALL DGEMV( 'T', MB, M-IE, -ONE, A( IS, IE+1 ),
00840      $                           LDA, RHS( 1 ), 1, ONE, C( IE+1, JS ),
00841      $                           1 )
00842                      CALL DGEMV( 'T', MB, M-IE, -ONE, D( IS, IE+1 ),
00843      $                           LDD, RHS( 3 ), 1, ONE, C( IE+1, JS ),
00844      $                           1 )
00845                   END IF
00846 *
00847                ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.2 ) ) THEN
00848 *
00849 *                 Build an 8-by-8 system Z' * x = RHS
00850 *
00851                   CALL DLASET( 'F', LDZ, LDZ, ZERO, ZERO, Z, LDZ )
00852 *
00853                   Z( 1, 1 ) = A( IS, IS )
00854                   Z( 2, 1 ) = A( IS, ISP1 )
00855                   Z( 5, 1 ) = -B( JS, JS )
00856                   Z( 7, 1 ) = -B( JSP1, JS )
00857 *
00858                   Z( 1, 2 ) = A( ISP1, IS )
00859                   Z( 2, 2 ) = A( ISP1, ISP1 )
00860                   Z( 6, 2 ) = -B( JS, JS )
00861                   Z( 8, 2 ) = -B( JSP1, JS )
00862 *
00863                   Z( 3, 3 ) = A( IS, IS )
00864                   Z( 4, 3 ) = A( IS, ISP1 )
00865                   Z( 5, 3 ) = -B( JS, JSP1 )
00866                   Z( 7, 3 ) = -B( JSP1, JSP1 )
00867 *
00868                   Z( 3, 4 ) = A( ISP1, IS )
00869                   Z( 4, 4 ) = A( ISP1, ISP1 )
00870                   Z( 6, 4 ) = -B( JS, JSP1 )
00871                   Z( 8, 4 ) = -B( JSP1, JSP1 )
00872 *
00873                   Z( 1, 5 ) = D( IS, IS )
00874                   Z( 2, 5 ) = D( IS, ISP1 )
00875                   Z( 5, 5 ) = -E( JS, JS )
00876 *
00877                   Z( 2, 6 ) = D( ISP1, ISP1 )
00878                   Z( 6, 6 ) = -E( JS, JS )
00879 *
00880                   Z( 3, 7 ) = D( IS, IS )
00881                   Z( 4, 7 ) = D( IS, ISP1 )
00882                   Z( 5, 7 ) = -E( JS, JSP1 )
00883                   Z( 7, 7 ) = -E( JSP1, JSP1 )
00884 *
00885                   Z( 4, 8 ) = D( ISP1, ISP1 )
00886                   Z( 6, 8 ) = -E( JS, JSP1 )
00887                   Z( 8, 8 ) = -E( JSP1, JSP1 )
00888 *
00889 *                 Set up right hand side(s)
00890 *
00891                   K = 1
00892                   II = MB*NB + 1
00893                   DO 160 JJ = 0, NB - 1
00894                      CALL DCOPY( MB, C( IS, JS+JJ ), 1, RHS( K ), 1 )
00895                      CALL DCOPY( MB, F( IS, JS+JJ ), 1, RHS( II ), 1 )
00896                      K = K + MB
00897                      II = II + MB
00898   160             CONTINUE
00899 *
00900 *
00901 *                 Solve Z' * x = RHS
00902 *
00903                   CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
00904                   IF( IERR.GT.0 )
00905      $               INFO = IERR
00906 *
00907                   CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
00908                   IF( SCALOC.NE.ONE ) THEN
00909                      DO 170 K = 1, N
00910                         CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
00911                         CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
00912   170                CONTINUE
00913                      SCALE = SCALE*SCALOC
00914                   END IF
00915 *
00916 *                 Unpack solution vector(s)
00917 *
00918                   K = 1
00919                   II = MB*NB + 1
00920                   DO 180 JJ = 0, NB - 1
00921                      CALL DCOPY( MB, RHS( K ), 1, C( IS, JS+JJ ), 1 )
00922                      CALL DCOPY( MB, RHS( II ), 1, F( IS, JS+JJ ), 1 )
00923                      K = K + MB
00924                      II = II + MB
00925   180             CONTINUE
00926 *
00927 *                 Substitute R(I, J) and L(I, J) into remaining
00928 *                 equation.
00929 *
00930                   IF( J.GT.P+2 ) THEN
00931                      CALL DGEMM( 'N', 'T', MB, JS-1, NB, ONE,
00932      $                           C( IS, JS ), LDC, B( 1, JS ), LDB, ONE,
00933      $                           F( IS, 1 ), LDF )
00934                      CALL DGEMM( 'N', 'T', MB, JS-1, NB, ONE,
00935      $                           F( IS, JS ), LDF, E( 1, JS ), LDE, ONE,
00936      $                           F( IS, 1 ), LDF )
00937                   END IF
00938                   IF( I.LT.P ) THEN
00939                      CALL DGEMM( 'T', 'N', M-IE, NB, MB, -ONE,
00940      $                           A( IS, IE+1 ), LDA, C( IS, JS ), LDC,
00941      $                           ONE, C( IE+1, JS ), LDC )
00942                      CALL DGEMM( 'T', 'N', M-IE, NB, MB, -ONE,
00943      $                           D( IS, IE+1 ), LDD, F( IS, JS ), LDF,
00944      $                           ONE, C( IE+1, JS ), LDC )
00945                   END IF
00946 *
00947                END IF
00948 *
00949   190       CONTINUE
00950   200    CONTINUE
00951 *
00952       END IF
00953       RETURN
00954 *
00955 *     End of DTGSY2
00956 *
00957       END