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