SRC/stgsy2.f

   1 *> \brief \b STGSY2 solves the generalized Sylvester equation (unblocked algorithm).
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download STGSY2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stgsy2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stgsy2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stgsy2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE STGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
  22 *                          LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL,
  23 *                          IWORK, PQ, INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       CHARACTER          TRANS
  27 *       INTEGER            IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N,
  28 *      $                   PQ
  29 *       REAL               RDSCAL, RDSUM, SCALE
  30 *       ..
  31 *       .. Array Arguments ..
  32 *       INTEGER            IWORK( * )
  33 *       REAL               A( LDA, * ), B( LDB, * ), C( LDC, * ),
  34 *      $                   D( LDD, * ), E( LDE, * ), F( LDF, * )
  35 *       ..
  36 *
  37 *
  38 *> \par Purpose:
  39 *  =============
  40 *>
  41 *> \verbatim
  42 *>
  43 *> STGSY2 solves the generalized Sylvester equation:
  44 *>
  45 *>             A * R - L * B = scale * C                (1)
  46 *>             D * R - L * E = scale * F,
  47 *>
  48 *> using Level 1 and 2 BLAS. where R and L are unknown M-by-N matrices,
  49 *> (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M,
  50 *> N-by-N and M-by-N, respectively, with real entries. (A, D) and (B, E)
  51 *> must be in generalized Schur canonical form, i.e. A, B are upper
  52 *> quasi triangular and D, E are upper triangular. The solution (R, L)
  53 *> overwrites (C, F). 0 <= SCALE <= 1 is an output scaling factor
  54 *> chosen to avoid overflow.
  55 *>
  56 *> In matrix notation solving equation (1) corresponds to solve
  57 *> Z*x = scale*b, where Z is defined as
  58 *>
  59 *>        Z = [ kron(In, A)  -kron(B**T, Im) ]             (2)
  60 *>            [ kron(In, D)  -kron(E**T, Im) ],
  61 *>
  62 *> Ik is the identity matrix of size k and X**T is the transpose of X.
  63 *> kron(X, Y) is the Kronecker product between the matrices X and Y.
  64 *> In the process of solving (1), we solve a number of such systems
  65 *> where Dim(In), Dim(In) = 1 or 2.
  66 *>
  67 *> If TRANS = 'T', solve the transposed system Z**T*y = scale*b for y,
  68 *> which is equivalent to solve for R and L in
  69 *>
  70 *>             A**T * R  + D**T * L   = scale * C           (3)
  71 *>             R  * B**T + L  * E**T  = scale * -F
  72 *>
  73 *> This case is used to compute an estimate of Dif[(A, D), (B, E)] =
  74 *> sigma_min(Z) using reverse communicaton with SLACON.
  75 *>
  76 *> STGSY2 also (IJOB >= 1) contributes to the computation in STGSYL
  77 *> of an upper bound on the separation between to matrix pairs. Then
  78 *> the input (A, D), (B, E) are sub-pencils of the matrix pair in
  79 *> STGSYL. See STGSYL for details.
  80 *> \endverbatim
  81 *
  82 *  Arguments:
  83 *  ==========
  84 *
  85 *> \param[in] TRANS
  86 *> \verbatim
  87 *>          TRANS is CHARACTER*1
  88 *>          = 'N', solve the generalized Sylvester equation (1).
  89 *>          = 'T': solve the 'transposed' system (3).
  90 *> \endverbatim
  91 *>
  92 *> \param[in] IJOB
  93 *> \verbatim
  94 *>          IJOB is INTEGER
  95 *>          Specifies what kind of functionality to be performed.
  96 *>          = 0: solve (1) only.
  97 *>          = 1: A contribution from this subsystem to a Frobenius
  98 *>               norm-based estimate of the separation between two matrix
  99 *>               pairs is computed. (look ahead strategy is used).
 100 *>          = 2: A contribution from this subsystem to a Frobenius
 101 *>               norm-based estimate of the separation between two matrix
 102 *>               pairs is computed. (SGECON on sub-systems is used.)
 103 *>          Not referenced if TRANS = 'T'.
 104 *> \endverbatim
 105 *>
 106 *> \param[in] M
 107 *> \verbatim
 108 *>          M is INTEGER
 109 *>          On entry, M specifies the order of A and D, and the row
 110 *>          dimension of C, F, R and L.
 111 *> \endverbatim
 112 *>
 113 *> \param[in] N
 114 *> \verbatim
 115 *>          N is INTEGER
 116 *>          On entry, N specifies the order of B and E, and the column
 117 *>          dimension of C, F, R and L.
 118 *> \endverbatim
 119 *>
 120 *> \param[in] A
 121 *> \verbatim
 122 *>          A is REAL array, dimension (LDA, M)
 123 *>          On entry, A contains an upper quasi triangular matrix.
 124 *> \endverbatim
 125 *>
 126 *> \param[in] LDA
 127 *> \verbatim
 128 *>          LDA is INTEGER
 129 *>          The leading dimension of the matrix A. LDA >= max(1, M).
 130 *> \endverbatim
 131 *>
 132 *> \param[in] B
 133 *> \verbatim
 134 *>          B is REAL array, dimension (LDB, N)
 135 *>          On entry, B contains an upper quasi triangular matrix.
 136 *> \endverbatim
 137 *>
 138 *> \param[in] LDB
 139 *> \verbatim
 140 *>          LDB is INTEGER
 141 *>          The leading dimension of the matrix B. LDB >= max(1, N).
 142 *> \endverbatim
 143 *>
 144 *> \param[in,out] C
 145 *> \verbatim
 146 *>          C is REAL array, dimension (LDC, N)
 147 *>          On entry, C contains the right-hand-side of the first matrix
 148 *>          equation in (1).
 149 *>          On exit, if IJOB = 0, C has been overwritten by the
 150 *>          solution R.
 151 *> \endverbatim
 152 *>
 153 *> \param[in] LDC
 154 *> \verbatim
 155 *>          LDC is INTEGER
 156 *>          The leading dimension of the matrix C. LDC >= max(1, M).
 157 *> \endverbatim
 158 *>
 159 *> \param[in] D
 160 *> \verbatim
 161 *>          D is REAL array, dimension (LDD, M)
 162 *>          On entry, D contains an upper triangular matrix.
 163 *> \endverbatim
 164 *>
 165 *> \param[in] LDD
 166 *> \verbatim
 167 *>          LDD is INTEGER
 168 *>          The leading dimension of the matrix D. LDD >= max(1, M).
 169 *> \endverbatim
 170 *>
 171 *> \param[in] E
 172 *> \verbatim
 173 *>          E is REAL array, dimension (LDE, N)
 174 *>          On entry, E contains an upper triangular matrix.
 175 *> \endverbatim
 176 *>
 177 *> \param[in] LDE
 178 *> \verbatim
 179 *>          LDE is INTEGER
 180 *>          The leading dimension of the matrix E. LDE >= max(1, N).
 181 *> \endverbatim
 182 *>
 183 *> \param[in,out] F
 184 *> \verbatim
 185 *>          F is REAL array, dimension (LDF, N)
 186 *>          On entry, F contains the right-hand-side of the second matrix
 187 *>          equation in (1).
 188 *>          On exit, if IJOB = 0, F has been overwritten by the
 189 *>          solution L.
 190 *> \endverbatim
 191 *>
 192 *> \param[in] LDF
 193 *> \verbatim
 194 *>          LDF is INTEGER
 195 *>          The leading dimension of the matrix F. LDF >= max(1, M).
 196 *> \endverbatim
 197 *>
 198 *> \param[out] SCALE
 199 *> \verbatim
 200 *>          SCALE is REAL
 201 *>          On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions
 202 *>          R and L (C and F on entry) will hold the solutions to a
 203 *>          slightly perturbed system but the input matrices A, B, D and
 204 *>          E have not been changed. If SCALE = 0, R and L will hold the
 205 *>          solutions to the homogeneous system with C = F = 0. Normally,
 206 *>          SCALE = 1.
 207 *> \endverbatim
 208 *>
 209 *> \param[in,out] RDSUM
 210 *> \verbatim
 211 *>          RDSUM is REAL
 212 *>          On entry, the sum of squares of computed contributions to
 213 *>          the Dif-estimate under computation by STGSYL, where the
 214 *>          scaling factor RDSCAL (see below) has been factored out.
 215 *>          On exit, the corresponding sum of squares updated with the
 216 *>          contributions from the current sub-system.
 217 *>          If TRANS = 'T' RDSUM is not touched.
 218 *>          NOTE: RDSUM only makes sense when STGSY2 is called by STGSYL.
 219 *> \endverbatim
 220 *>
 221 *> \param[in,out] RDSCAL
 222 *> \verbatim
 223 *>          RDSCAL is REAL
 224 *>          On entry, scaling factor used to prevent overflow in RDSUM.
 225 *>          On exit, RDSCAL is updated w.r.t. the current contributions
 226 *>          in RDSUM.
 227 *>          If TRANS = 'T', RDSCAL is not touched.
 228 *>          NOTE: RDSCAL only makes sense when STGSY2 is called by
 229 *>                STGSYL.
 230 *> \endverbatim
 231 *>
 232 *> \param[out] IWORK
 233 *> \verbatim
 234 *>          IWORK is INTEGER array, dimension (M+N+2)
 235 *> \endverbatim
 236 *>
 237 *> \param[out] PQ
 238 *> \verbatim
 239 *>          PQ is INTEGER
 240 *>          On exit, the number of subsystems (of size 2-by-2, 4-by-4 and
 241 *>          8-by-8) solved by this routine.
 242 *> \endverbatim
 243 *>
 244 *> \param[out] INFO
 245 *> \verbatim
 246 *>          INFO is INTEGER
 247 *>          On exit, if INFO is set to
 248 *>            =0: Successful exit
 249 *>            <0: If INFO = -i, the i-th argument had an illegal value.
 250 *>            >0: The matrix pairs (A, D) and (B, E) have common or very
 251 *>                close eigenvalues.
 252 *> \endverbatim
 253 *
 254 *  Authors:
 255 *  ========
 256 *
 257 *> \author Univ. of Tennessee
 258 *> \author Univ. of California Berkeley
 259 *> \author Univ. of Colorado Denver
 260 *> \author NAG Ltd.
 261 *
 262 *> \date November 2015
 263 *
 264 *> \ingroup realSYauxiliary
 265 *
 266 *> \par Contributors:
 267 *  ==================
 268 *>
 269 *>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
 270 *>     Umea University, S-901 87 Umea, Sweden.
 271 *
 272 *  =====================================================================
 273       SUBROUTINE STGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
 274      $                   LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL,
 275      $                   IWORK, PQ, INFO )
 276 *
 277 *  -- LAPACK auxiliary routine (version 3.6.0) --
 278 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 279 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 280 *     November 2015
 281 *
 282 *     .. Scalar Arguments ..
 283       CHARACTER          TRANS
 284       INTEGER            IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N,
 285      $                   PQ
 286       REAL               RDSCAL, RDSUM, SCALE
 287 *     ..
 288 *     .. Array Arguments ..
 289       INTEGER            IWORK( * )
 290       REAL               A( LDA, * ), B( LDB, * ), C( LDC, * ),
 291      $                   D( LDD, * ), E( LDE, * ), F( LDF, * )
 292 *     ..
 293 *
 294 *  =====================================================================
 295 *  Replaced various illegal calls to SCOPY by calls to SLASET.
 296 *  Sven Hammarling, 27/5/02.
 297 *
 298 *     .. Parameters ..
 299       INTEGER            LDZ
 300       PARAMETER          ( LDZ = 8 )
 301       REAL               ZERO, ONE
 302       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 303 *     ..
 304 *     .. Local Scalars ..
 305       LOGICAL            NOTRAN
 306       INTEGER            I, IE, IERR, II, IS, ISP1, J, JE, JJ, JS, JSP1,
 307      $                   K, MB, NB, P, Q, ZDIM
 308       REAL               ALPHA, SCALOC
 309 *     ..
 310 *     .. Local Arrays ..
 311       INTEGER            IPIV( LDZ ), JPIV( LDZ )
 312       REAL               RHS( LDZ ), Z( LDZ, LDZ )
 313 *     ..
 314 *     .. External Functions ..
 315       LOGICAL            LSAME
 316       EXTERNAL           LSAME
 317 *     ..
 318 *     .. External Subroutines ..
 319       EXTERNAL           SAXPY, SCOPY, SGEMM, SGEMV, SGER, SGESC2,
 320      $                   SGETC2, SSCAL, SLASET, SLATDF, XERBLA
 321 *     ..
 322 *     .. Intrinsic Functions ..
 323       INTRINSIC          MAX
 324 *     ..
 325 *     .. Executable Statements ..
 326 *
 327 *     Decode and test input parameters
 328 *
 329       INFO = 0
 330       IERR = 0
 331       NOTRAN = LSAME( TRANS, 'N' )
 332       IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
 333          INFO = -1
 334       ELSE IF( NOTRAN ) THEN
 335          IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.2 ) ) THEN
 336             INFO = -2
 337          END IF
 338       END IF
 339       IF( INFO.EQ.0 ) THEN
 340          IF( M.LE.0 ) THEN
 341             INFO = -3
 342          ELSE IF( N.LE.0 ) THEN
 343             INFO = -4
 344          ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 345             INFO = -6
 346          ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 347             INFO = -8
 348          ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
 349             INFO = -10
 350          ELSE IF( LDD.LT.MAX( 1, M ) ) THEN
 351             INFO = -12
 352          ELSE IF( LDE.LT.MAX( 1, N ) ) THEN
 353             INFO = -14
 354          ELSE IF( LDF.LT.MAX( 1, M ) ) THEN
 355             INFO = -16
 356          END IF
 357       END IF
 358       IF( INFO.NE.0 ) THEN
 359          CALL XERBLA( 'STGSY2', -INFO )
 360          RETURN
 361       END IF
 362 *
 363 *     Determine block structure of A
 364 *
 365       PQ = 0
 366       P = 0
 367       I = 1
 368    10 CONTINUE
 369       IF( I.GT.M )
 370      $   GO TO 20
 371       P = P + 1
 372       IWORK( P ) = I
 373       IF( I.EQ.M )
 374      $   GO TO 20
 375       IF( A( I+1, I ).NE.ZERO ) THEN
 376          I = I + 2
 377       ELSE
 378          I = I + 1
 379       END IF
 380       GO TO 10
 381    20 CONTINUE
 382       IWORK( P+1 ) = M + 1
 383 *
 384 *     Determine block structure of B
 385 *
 386       Q = P + 1
 387       J = 1
 388    30 CONTINUE
 389       IF( J.GT.N )
 390      $   GO TO 40
 391       Q = Q + 1
 392       IWORK( Q ) = J
 393       IF( J.EQ.N )
 394      $   GO TO 40
 395       IF( B( J+1, J ).NE.ZERO ) THEN
 396          J = J + 2
 397       ELSE
 398          J = J + 1
 399       END IF
 400       GO TO 30
 401    40 CONTINUE
 402       IWORK( Q+1 ) = N + 1
 403       PQ = P*( Q-P-1 )
 404 *
 405       IF( NOTRAN ) THEN
 406 *
 407 *        Solve (I, J) - subsystem
 408 *           A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
 409 *           D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
 410 *        for I = P, P - 1, ..., 1; J = 1, 2, ..., Q
 411 *
 412          SCALE = ONE
 413          SCALOC = ONE
 414          DO 120 J = P + 2, Q
 415             JS = IWORK( J )
 416             JSP1 = JS + 1
 417             JE = IWORK( J+1 ) - 1
 418             NB = JE - JS + 1
 419             DO 110 I = P, 1, -1
 420 *
 421                IS = IWORK( I )
 422                ISP1 = IS + 1
 423                IE = IWORK( I+1 ) - 1
 424                MB = IE - IS + 1
 425                ZDIM = MB*NB*2
 426 *
 427                IF( ( MB.EQ.1 ) .AND. ( NB.EQ.1 ) ) THEN
 428 *
 429 *                 Build a 2-by-2 system Z * x = RHS
 430 *
 431                   Z( 1, 1 ) = A( IS, IS )
 432                   Z( 2, 1 ) = D( IS, IS )
 433                   Z( 1, 2 ) = -B( JS, JS )
 434                   Z( 2, 2 ) = -E( JS, JS )
 435 *
 436 *                 Set up right hand side(s)
 437 *
 438                   RHS( 1 ) = C( IS, JS )
 439                   RHS( 2 ) = F( IS, JS )
 440 *
 441 *                 Solve Z * x = RHS
 442 *
 443                   CALL SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
 444                   IF( IERR.GT.0 )
 445      $               INFO = IERR
 446 *
 447                   IF( IJOB.EQ.0 ) THEN
 448                      CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
 449      $                            SCALOC )
 450                      IF( SCALOC.NE.ONE ) THEN
 451                         DO 50 K = 1, N
 452                            CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
 453                            CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
 454    50                   CONTINUE
 455                         SCALE = SCALE*SCALOC
 456                      END IF
 457                   ELSE
 458                      CALL SLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM,
 459      $                            RDSCAL, IPIV, JPIV )
 460                   END IF
 461 *
 462 *                 Unpack solution vector(s)
 463 *
 464                   C( IS, JS ) = RHS( 1 )
 465                   F( IS, JS ) = RHS( 2 )
 466 *
 467 *                 Substitute R(I, J) and L(I, J) into remaining
 468 *                 equation.
 469 *
 470                   IF( I.GT.1 ) THEN
 471                      ALPHA = -RHS( 1 )
 472                      CALL SAXPY( IS-1, ALPHA, A( 1, IS ), 1, C( 1, JS ),
 473      $                           1 )
 474                      CALL SAXPY( IS-1, ALPHA, D( 1, IS ), 1, F( 1, JS ),
 475      $                           1 )
 476                   END IF
 477                   IF( J.LT.Q ) THEN
 478                      CALL SAXPY( N-JE, RHS( 2 ), B( JS, JE+1 ), LDB,
 479      $                           C( IS, JE+1 ), LDC )
 480                      CALL SAXPY( N-JE, RHS( 2 ), E( JS, JE+1 ), LDE,
 481      $                           F( IS, JE+1 ), LDF )
 482                   END IF
 483 *
 484                ELSE IF( ( MB.EQ.1 ) .AND. ( NB.EQ.2 ) ) THEN
 485 *
 486 *                 Build a 4-by-4 system Z * x = RHS
 487 *
 488                   Z( 1, 1 ) = A( IS, IS )
 489                   Z( 2, 1 ) = ZERO
 490                   Z( 3, 1 ) = D( IS, IS )
 491                   Z( 4, 1 ) = ZERO
 492 *
 493                   Z( 1, 2 ) = ZERO
 494                   Z( 2, 2 ) = A( IS, IS )
 495                   Z( 3, 2 ) = ZERO
 496                   Z( 4, 2 ) = D( IS, IS )
 497 *
 498                   Z( 1, 3 ) = -B( JS, JS )
 499                   Z( 2, 3 ) = -B( JS, JSP1 )
 500                   Z( 3, 3 ) = -E( JS, JS )
 501                   Z( 4, 3 ) = -E( JS, JSP1 )
 502 *
 503                   Z( 1, 4 ) = -B( JSP1, JS )
 504                   Z( 2, 4 ) = -B( JSP1, JSP1 )
 505                   Z( 3, 4 ) = ZERO
 506                   Z( 4, 4 ) = -E( JSP1, JSP1 )
 507 *
 508 *                 Set up right hand side(s)
 509 *
 510                   RHS( 1 ) = C( IS, JS )
 511                   RHS( 2 ) = C( IS, JSP1 )
 512                   RHS( 3 ) = F( IS, JS )
 513                   RHS( 4 ) = F( IS, JSP1 )
 514 *
 515 *                 Solve Z * x = RHS
 516 *
 517                   CALL SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
 518                   IF( IERR.GT.0 )
 519      $               INFO = IERR
 520 *
 521                   IF( IJOB.EQ.0 ) THEN
 522                      CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
 523      $                            SCALOC )
 524                      IF( SCALOC.NE.ONE ) THEN
 525                         DO 60 K = 1, N
 526                            CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
 527                            CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
 528    60                   CONTINUE
 529                         SCALE = SCALE*SCALOC
 530                      END IF
 531                   ELSE
 532                      CALL SLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM,
 533      $                            RDSCAL, IPIV, JPIV )
 534                   END IF
 535 *
 536 *                 Unpack solution vector(s)
 537 *
 538                   C( IS, JS ) = RHS( 1 )
 539                   C( IS, JSP1 ) = RHS( 2 )
 540                   F( IS, JS ) = RHS( 3 )
 541                   F( IS, JSP1 ) = RHS( 4 )
 542 *
 543 *                 Substitute R(I, J) and L(I, J) into remaining
 544 *                 equation.
 545 *
 546                   IF( I.GT.1 ) THEN
 547                      CALL SGER( IS-1, NB, -ONE, A( 1, IS ), 1, RHS( 1 ),
 548      $                          1, C( 1, JS ), LDC )
 549                      CALL SGER( IS-1, NB, -ONE, D( 1, IS ), 1, RHS( 1 ),
 550      $                          1, F( 1, JS ), LDF )
 551                   END IF
 552                   IF( J.LT.Q ) THEN
 553                      CALL SAXPY( N-JE, RHS( 3 ), B( JS, JE+1 ), LDB,
 554      $                           C( IS, JE+1 ), LDC )
 555                      CALL SAXPY( N-JE, RHS( 3 ), E( JS, JE+1 ), LDE,
 556      $                           F( IS, JE+1 ), LDF )
 557                      CALL SAXPY( N-JE, RHS( 4 ), B( JSP1, JE+1 ), LDB,
 558      $                           C( IS, JE+1 ), LDC )
 559                      CALL SAXPY( N-JE, RHS( 4 ), E( JSP1, JE+1 ), LDE,
 560      $                           F( IS, JE+1 ), LDF )
 561                   END IF
 562 *
 563                ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.1 ) ) THEN
 564 *
 565 *                 Build a 4-by-4 system Z * x = RHS
 566 *
 567                   Z( 1, 1 ) = A( IS, IS )
 568                   Z( 2, 1 ) = A( ISP1, IS )
 569                   Z( 3, 1 ) = D( IS, IS )
 570                   Z( 4, 1 ) = ZERO
 571 *
 572                   Z( 1, 2 ) = A( IS, ISP1 )
 573                   Z( 2, 2 ) = A( ISP1, ISP1 )
 574                   Z( 3, 2 ) = D( IS, ISP1 )
 575                   Z( 4, 2 ) = D( ISP1, ISP1 )
 576 *
 577                   Z( 1, 3 ) = -B( JS, JS )
 578                   Z( 2, 3 ) = ZERO
 579                   Z( 3, 3 ) = -E( JS, JS )
 580                   Z( 4, 3 ) = ZERO
 581 *
 582                   Z( 1, 4 ) = ZERO
 583                   Z( 2, 4 ) = -B( JS, JS )
 584                   Z( 3, 4 ) = ZERO
 585                   Z( 4, 4 ) = -E( JS, JS )
 586 *
 587 *                 Set up right hand side(s)
 588 *
 589                   RHS( 1 ) = C( IS, JS )
 590                   RHS( 2 ) = C( ISP1, JS )
 591                   RHS( 3 ) = F( IS, JS )
 592                   RHS( 4 ) = F( ISP1, JS )
 593 *
 594 *                 Solve Z * x = RHS
 595 *
 596                   CALL SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
 597                   IF( IERR.GT.0 )
 598      $               INFO = IERR
 599                   IF( IJOB.EQ.0 ) THEN
 600                      CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
 601      $                            SCALOC )
 602                      IF( SCALOC.NE.ONE ) THEN
 603                         DO 70 K = 1, N
 604                            CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
 605                            CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
 606    70                   CONTINUE
 607                         SCALE = SCALE*SCALOC
 608                      END IF
 609                   ELSE
 610                      CALL SLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM,
 611      $                            RDSCAL, IPIV, JPIV )
 612                   END IF
 613 *
 614 *                 Unpack solution vector(s)
 615 *
 616                   C( IS, JS ) = RHS( 1 )
 617                   C( ISP1, JS ) = RHS( 2 )
 618                   F( IS, JS ) = RHS( 3 )
 619                   F( ISP1, JS ) = RHS( 4 )
 620 *
 621 *                 Substitute R(I, J) and L(I, J) into remaining
 622 *                 equation.
 623 *
 624                   IF( I.GT.1 ) THEN
 625                      CALL SGEMV( 'N', IS-1, MB, -ONE, A( 1, IS ), LDA,
 626      $                           RHS( 1 ), 1, ONE, C( 1, JS ), 1 )
 627                      CALL SGEMV( 'N', IS-1, MB, -ONE, D( 1, IS ), LDD,
 628      $                           RHS( 1 ), 1, ONE, F( 1, JS ), 1 )
 629                   END IF
 630                   IF( J.LT.Q ) THEN
 631                      CALL SGER( MB, N-JE, ONE, RHS( 3 ), 1,
 632      $                          B( JS, JE+1 ), LDB, C( IS, JE+1 ), LDC )
 633                      CALL SGER( MB, N-JE, ONE, RHS( 3 ), 1,
 634      $                          E( JS, JE+1 ), LDE, F( IS, JE+1 ), LDF )
 635                   END IF
 636 *
 637                ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.2 ) ) THEN
 638 *
 639 *                 Build an 8-by-8 system Z * x = RHS
 640 *
 641                   CALL SLASET( 'F', LDZ, LDZ, ZERO, ZERO, Z, LDZ )
 642 *
 643                   Z( 1, 1 ) = A( IS, IS )
 644                   Z( 2, 1 ) = A( ISP1, IS )
 645                   Z( 5, 1 ) = D( IS, IS )
 646 *
 647                   Z( 1, 2 ) = A( IS, ISP1 )
 648                   Z( 2, 2 ) = A( ISP1, ISP1 )
 649                   Z( 5, 2 ) = D( IS, ISP1 )
 650                   Z( 6, 2 ) = D( ISP1, ISP1 )
 651 *
 652                   Z( 3, 3 ) = A( IS, IS )
 653                   Z( 4, 3 ) = A( ISP1, IS )
 654                   Z( 7, 3 ) = D( IS, IS )
 655 *
 656                   Z( 3, 4 ) = A( IS, ISP1 )
 657                   Z( 4, 4 ) = A( ISP1, ISP1 )
 658                   Z( 7, 4 ) = D( IS, ISP1 )
 659                   Z( 8, 4 ) = D( ISP1, ISP1 )
 660 *
 661                   Z( 1, 5 ) = -B( JS, JS )
 662                   Z( 3, 5 ) = -B( JS, JSP1 )
 663                   Z( 5, 5 ) = -E( JS, JS )
 664                   Z( 7, 5 ) = -E( JS, JSP1 )
 665 *
 666                   Z( 2, 6 ) = -B( JS, JS )
 667                   Z( 4, 6 ) = -B( JS, JSP1 )
 668                   Z( 6, 6 ) = -E( JS, JS )
 669                   Z( 8, 6 ) = -E( JS, JSP1 )
 670 *
 671                   Z( 1, 7 ) = -B( JSP1, JS )
 672                   Z( 3, 7 ) = -B( JSP1, JSP1 )
 673                   Z( 7, 7 ) = -E( JSP1, JSP1 )
 674 *
 675                   Z( 2, 8 ) = -B( JSP1, JS )
 676                   Z( 4, 8 ) = -B( JSP1, JSP1 )
 677                   Z( 8, 8 ) = -E( JSP1, JSP1 )
 678 *
 679 *                 Set up right hand side(s)
 680 *
 681                   K = 1
 682                   II = MB*NB + 1
 683                   DO 80 JJ = 0, NB - 1
 684                      CALL SCOPY( MB, C( IS, JS+JJ ), 1, RHS( K ), 1 )
 685                      CALL SCOPY( MB, F( IS, JS+JJ ), 1, RHS( II ), 1 )
 686                      K = K + MB
 687                      II = II + MB
 688    80             CONTINUE
 689 *
 690 *                 Solve Z * x = RHS
 691 *
 692                   CALL SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
 693                   IF( IERR.GT.0 )
 694      $               INFO = IERR
 695                   IF( IJOB.EQ.0 ) THEN
 696                      CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
 697      $                            SCALOC )
 698                      IF( SCALOC.NE.ONE ) THEN
 699                         DO 90 K = 1, N
 700                            CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
 701                            CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
 702    90                   CONTINUE
 703                         SCALE = SCALE*SCALOC
 704                      END IF
 705                   ELSE
 706                      CALL SLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM,
 707      $                            RDSCAL, IPIV, JPIV )
 708                   END IF
 709 *
 710 *                 Unpack solution vector(s)
 711 *
 712                   K = 1
 713                   II = MB*NB + 1
 714                   DO 100 JJ = 0, NB - 1
 715                      CALL SCOPY( MB, RHS( K ), 1, C( IS, JS+JJ ), 1 )
 716                      CALL SCOPY( MB, RHS( II ), 1, F( IS, JS+JJ ), 1 )
 717                      K = K + MB
 718                      II = II + MB
 719   100             CONTINUE
 720 *
 721 *                 Substitute R(I, J) and L(I, J) into remaining
 722 *                 equation.
 723 *
 724                   IF( I.GT.1 ) THEN
 725                      CALL SGEMM( 'N', 'N', IS-1, NB, MB, -ONE,
 726      $                           A( 1, IS ), LDA, RHS( 1 ), MB, ONE,
 727      $                           C( 1, JS ), LDC )
 728                      CALL SGEMM( 'N', 'N', IS-1, NB, MB, -ONE,
 729      $                           D( 1, IS ), LDD, RHS( 1 ), MB, ONE,
 730      $                           F( 1, JS ), LDF )
 731                   END IF
 732                   IF( J.LT.Q ) THEN
 733                      K = MB*NB + 1
 734                      CALL SGEMM( 'N', 'N', MB, N-JE, NB, ONE, RHS( K ),
 735      $                           MB, B( JS, JE+1 ), LDB, ONE,
 736      $                           C( IS, JE+1 ), LDC )
 737                      CALL SGEMM( 'N', 'N', MB, N-JE, NB, ONE, RHS( K ),
 738      $                           MB, E( JS, JE+1 ), LDE, ONE,
 739      $                           F( IS, JE+1 ), LDF )
 740                   END IF
 741 *
 742                END IF
 743 *
 744   110       CONTINUE
 745   120    CONTINUE
 746       ELSE
 747 *
 748 *        Solve (I, J) - subsystem
 749 *             A(I, I)**T * R(I, J) + D(I, I)**T * L(J, J)  =  C(I, J)
 750 *             R(I, I)  * B(J, J) + L(I, J)  * E(J, J)  = -F(I, J)
 751 *        for I = 1, 2, ..., P, J = Q, Q - 1, ..., 1
 752 *
 753          SCALE = ONE
 754          SCALOC = ONE
 755          DO 200 I = 1, P
 756 *
 757             IS = IWORK( I )
 758             ISP1 = IS + 1
 759             IE = IWORK( I+1 ) - 1
 760             MB = IE - IS + 1
 761             DO 190 J = Q, P + 2, -1
 762 *
 763                JS = IWORK( J )
 764                JSP1 = JS + 1
 765                JE = IWORK( J+1 ) - 1
 766                NB = JE - JS + 1
 767                ZDIM = MB*NB*2
 768                IF( ( MB.EQ.1 ) .AND. ( NB.EQ.1 ) ) THEN
 769 *
 770 *                 Build a 2-by-2 system Z**T * x = RHS
 771 *
 772                   Z( 1, 1 ) = A( IS, IS )
 773                   Z( 2, 1 ) = -B( JS, JS )
 774                   Z( 1, 2 ) = D( IS, IS )
 775                   Z( 2, 2 ) = -E( JS, JS )
 776 *
 777 *                 Set up right hand side(s)
 778 *
 779                   RHS( 1 ) = C( IS, JS )
 780                   RHS( 2 ) = F( IS, JS )
 781 *
 782 *                 Solve Z**T * x = RHS
 783 *
 784                   CALL SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
 785                   IF( IERR.GT.0 )
 786      $               INFO = IERR
 787 *
 788                   CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
 789                   IF( SCALOC.NE.ONE ) THEN
 790                      DO 130 K = 1, N
 791                         CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
 792                         CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
 793   130                CONTINUE
 794                      SCALE = SCALE*SCALOC
 795                   END IF
 796 *
 797 *                 Unpack solution vector(s)
 798 *
 799                   C( IS, JS ) = RHS( 1 )
 800                   F( IS, JS ) = RHS( 2 )
 801 *
 802 *                 Substitute R(I, J) and L(I, J) into remaining
 803 *                 equation.
 804 *
 805                   IF( J.GT.P+2 ) THEN
 806                      ALPHA = RHS( 1 )
 807                      CALL SAXPY( JS-1, ALPHA, B( 1, JS ), 1, F( IS, 1 ),
 808      $                           LDF )
 809                      ALPHA = RHS( 2 )
 810                      CALL SAXPY( JS-1, ALPHA, E( 1, JS ), 1, F( IS, 1 ),
 811      $                           LDF )
 812                   END IF
 813                   IF( I.LT.P ) THEN
 814                      ALPHA = -RHS( 1 )
 815                      CALL SAXPY( M-IE, ALPHA, A( IS, IE+1 ), LDA,
 816      $                           C( IE+1, JS ), 1 )
 817                      ALPHA = -RHS( 2 )
 818                      CALL SAXPY( M-IE, ALPHA, D( IS, IE+1 ), LDD,
 819      $                           C( IE+1, JS ), 1 )
 820                   END IF
 821 *
 822                ELSE IF( ( MB.EQ.1 ) .AND. ( NB.EQ.2 ) ) THEN
 823 *
 824 *                 Build a 4-by-4 system Z**T * x = RHS
 825 *
 826                   Z( 1, 1 ) = A( IS, IS )
 827                   Z( 2, 1 ) = ZERO
 828                   Z( 3, 1 ) = -B( JS, JS )
 829                   Z( 4, 1 ) = -B( JSP1, JS )
 830 *
 831                   Z( 1, 2 ) = ZERO
 832                   Z( 2, 2 ) = A( IS, IS )
 833                   Z( 3, 2 ) = -B( JS, JSP1 )
 834                   Z( 4, 2 ) = -B( JSP1, JSP1 )
 835 *
 836                   Z( 1, 3 ) = D( IS, IS )
 837                   Z( 2, 3 ) = ZERO
 838                   Z( 3, 3 ) = -E( JS, JS )
 839                   Z( 4, 3 ) = ZERO
 840 *
 841                   Z( 1, 4 ) = ZERO
 842                   Z( 2, 4 ) = D( IS, IS )
 843                   Z( 3, 4 ) = -E( JS, JSP1 )
 844                   Z( 4, 4 ) = -E( JSP1, JSP1 )
 845 *
 846 *                 Set up right hand side(s)
 847 *
 848                   RHS( 1 ) = C( IS, JS )
 849                   RHS( 2 ) = C( IS, JSP1 )
 850                   RHS( 3 ) = F( IS, JS )
 851                   RHS( 4 ) = F( IS, JSP1 )
 852 *
 853 *                 Solve Z**T * x = RHS
 854 *
 855                   CALL SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
 856                   IF( IERR.GT.0 )
 857      $               INFO = IERR
 858                   CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
 859                   IF( SCALOC.NE.ONE ) THEN
 860                      DO 140 K = 1, N
 861                         CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
 862                         CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
 863   140                CONTINUE
 864                      SCALE = SCALE*SCALOC
 865                   END IF
 866 *
 867 *                 Unpack solution vector(s)
 868 *
 869                   C( IS, JS ) = RHS( 1 )
 870                   C( IS, JSP1 ) = RHS( 2 )
 871                   F( IS, JS ) = RHS( 3 )
 872                   F( IS, JSP1 ) = RHS( 4 )
 873 *
 874 *                 Substitute R(I, J) and L(I, J) into remaining
 875 *                 equation.
 876 *
 877                   IF( J.GT.P+2 ) THEN
 878                      CALL SAXPY( JS-1, RHS( 1 ), B( 1, JS ), 1,
 879      $                           F( IS, 1 ), LDF )
 880                      CALL SAXPY( JS-1, RHS( 2 ), B( 1, JSP1 ), 1,
 881      $                           F( IS, 1 ), LDF )
 882                      CALL SAXPY( JS-1, RHS( 3 ), E( 1, JS ), 1,
 883      $                           F( IS, 1 ), LDF )
 884                      CALL SAXPY( JS-1, RHS( 4 ), E( 1, JSP1 ), 1,
 885      $                           F( IS, 1 ), LDF )
 886                   END IF
 887                   IF( I.LT.P ) THEN
 888                      CALL SGER( M-IE, NB, -ONE, A( IS, IE+1 ), LDA,
 889      $                          RHS( 1 ), 1, C( IE+1, JS ), LDC )
 890                      CALL SGER( M-IE, NB, -ONE, D( IS, IE+1 ), LDD,
 891      $                          RHS( 3 ), 1, C( IE+1, JS ), LDC )
 892                   END IF
 893 *
 894                ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.1 ) ) THEN
 895 *
 896 *                 Build a 4-by-4 system Z**T * x = RHS
 897 *
 898                   Z( 1, 1 ) = A( IS, IS )
 899                   Z( 2, 1 ) = A( IS, ISP1 )
 900                   Z( 3, 1 ) = -B( JS, JS )
 901                   Z( 4, 1 ) = ZERO
 902 *
 903                   Z( 1, 2 ) = A( ISP1, IS )
 904                   Z( 2, 2 ) = A( ISP1, ISP1 )
 905                   Z( 3, 2 ) = ZERO
 906                   Z( 4, 2 ) = -B( JS, JS )
 907 *
 908                   Z( 1, 3 ) = D( IS, IS )
 909                   Z( 2, 3 ) = D( IS, ISP1 )
 910                   Z( 3, 3 ) = -E( JS, JS )
 911                   Z( 4, 3 ) = ZERO
 912 *
 913                   Z( 1, 4 ) = ZERO
 914                   Z( 2, 4 ) = D( ISP1, ISP1 )
 915                   Z( 3, 4 ) = ZERO
 916                   Z( 4, 4 ) = -E( JS, JS )
 917 *
 918 *                 Set up right hand side(s)
 919 *
 920                   RHS( 1 ) = C( IS, JS )
 921                   RHS( 2 ) = C( ISP1, JS )
 922                   RHS( 3 ) = F( IS, JS )
 923                   RHS( 4 ) = F( ISP1, JS )
 924 *
 925 *                 Solve Z**T * x = RHS
 926 *
 927                   CALL SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
 928                   IF( IERR.GT.0 )
 929      $               INFO = IERR
 930 *
 931                   CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
 932                   IF( SCALOC.NE.ONE ) THEN
 933                      DO 150 K = 1, N
 934                         CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
 935                         CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
 936   150                CONTINUE
 937                      SCALE = SCALE*SCALOC
 938                   END IF
 939 *
 940 *                 Unpack solution vector(s)
 941 *
 942                   C( IS, JS ) = RHS( 1 )
 943                   C( ISP1, JS ) = RHS( 2 )
 944                   F( IS, JS ) = RHS( 3 )
 945                   F( ISP1, JS ) = RHS( 4 )
 946 *
 947 *                 Substitute R(I, J) and L(I, J) into remaining
 948 *                 equation.
 949 *
 950                   IF( J.GT.P+2 ) THEN
 951                      CALL SGER( MB, JS-1, ONE, RHS( 1 ), 1, B( 1, JS ),
 952      $                          1, F( IS, 1 ), LDF )
 953                      CALL SGER( MB, JS-1, ONE, RHS( 3 ), 1, E( 1, JS ),
 954      $                          1, F( IS, 1 ), LDF )
 955                   END IF
 956                   IF( I.LT.P ) THEN
 957                      CALL SGEMV( 'T', MB, M-IE, -ONE, A( IS, IE+1 ),
 958      $                           LDA, RHS( 1 ), 1, ONE, C( IE+1, JS ),
 959      $                           1 )
 960                      CALL SGEMV( 'T', MB, M-IE, -ONE, D( IS, IE+1 ),
 961      $                           LDD, RHS( 3 ), 1, ONE, C( IE+1, JS ),
 962      $                           1 )
 963                   END IF
 964 *
 965                ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.2 ) ) THEN
 966 *
 967 *                 Build an 8-by-8 system Z**T * x = RHS
 968 *
 969                   CALL SLASET( 'F', LDZ, LDZ, ZERO, ZERO, Z, LDZ )
 970 *
 971                   Z( 1, 1 ) = A( IS, IS )
 972                   Z( 2, 1 ) = A( IS, ISP1 )
 973                   Z( 5, 1 ) = -B( JS, JS )
 974                   Z( 7, 1 ) = -B( JSP1, JS )
 975 *
 976                   Z( 1, 2 ) = A( ISP1, IS )
 977                   Z( 2, 2 ) = A( ISP1, ISP1 )
 978                   Z( 6, 2 ) = -B( JS, JS )
 979                   Z( 8, 2 ) = -B( JSP1, JS )
 980 *
 981                   Z( 3, 3 ) = A( IS, IS )
 982                   Z( 4, 3 ) = A( IS, ISP1 )
 983                   Z( 5, 3 ) = -B( JS, JSP1 )
 984                   Z( 7, 3 ) = -B( JSP1, JSP1 )
 985 *
 986                   Z( 3, 4 ) = A( ISP1, IS )
 987                   Z( 4, 4 ) = A( ISP1, ISP1 )
 988                   Z( 6, 4 ) = -B( JS, JSP1 )
 989                   Z( 8, 4 ) = -B( JSP1, JSP1 )
 990 *
 991                   Z( 1, 5 ) = D( IS, IS )
 992                   Z( 2, 5 ) = D( IS, ISP1 )
 993                   Z( 5, 5 ) = -E( JS, JS )
 994 *
 995                   Z( 2, 6 ) = D( ISP1, ISP1 )
 996                   Z( 6, 6 ) = -E( JS, JS )
 997 *
 998                   Z( 3, 7 ) = D( IS, IS )
 999                   Z( 4, 7 ) = D( IS, ISP1 )
1000                   Z( 5, 7 ) = -E( JS, JSP1 )
1001                   Z( 7, 7 ) = -E( JSP1, JSP1 )
1002 *
1003                   Z( 4, 8 ) = D( ISP1, ISP1 )
1004                   Z( 6, 8 ) = -E( JS, JSP1 )
1005                   Z( 8, 8 ) = -E( JSP1, JSP1 )
1006 *
1007 *                 Set up right hand side(s)
1008 *
1009                   K = 1
1010                   II = MB*NB + 1
1011                   DO 160 JJ = 0, NB - 1
1012                      CALL SCOPY( MB, C( IS, JS+JJ ), 1, RHS( K ), 1 )
1013                      CALL SCOPY( MB, F( IS, JS+JJ ), 1, RHS( II ), 1 )
1014                      K = K + MB
1015                      II = II + MB
1016   160             CONTINUE
1017 *
1018 *
1019 *                 Solve Z**T * x = RHS
1020 *
1021                   CALL SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
1022                   IF( IERR.GT.0 )
1023      $               INFO = IERR
1024 *
1025                   CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
1026                   IF( SCALOC.NE.ONE ) THEN
1027                      DO 170 K = 1, N
1028                         CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
1029                         CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
1030   170                CONTINUE
1031                      SCALE = SCALE*SCALOC
1032                   END IF
1033 *
1034 *                 Unpack solution vector(s)
1035 *
1036                   K = 1
1037                   II = MB*NB + 1
1038                   DO 180 JJ = 0, NB - 1
1039                      CALL SCOPY( MB, RHS( K ), 1, C( IS, JS+JJ ), 1 )
1040                      CALL SCOPY( MB, RHS( II ), 1, F( IS, JS+JJ ), 1 )
1041                      K = K + MB
1042                      II = II + MB
1043   180             CONTINUE
1044 *
1045 *                 Substitute R(I, J) and L(I, J) into remaining
1046 *                 equation.
1047 *
1048                   IF( J.GT.P+2 ) THEN
1049                      CALL SGEMM( 'N', 'T', MB, JS-1, NB, ONE,
1050      $                           C( IS, JS ), LDC, B( 1, JS ), LDB, ONE,
1051      $                           F( IS, 1 ), LDF )
1052                      CALL SGEMM( 'N', 'T', MB, JS-1, NB, ONE,
1053      $                           F( IS, JS ), LDF, E( 1, JS ), LDE, ONE,
1054      $                           F( IS, 1 ), LDF )
1055                   END IF
1056                   IF( I.LT.P ) THEN
1057                      CALL SGEMM( 'T', 'N', M-IE, NB, MB, -ONE,
1058      $                           A( IS, IE+1 ), LDA, C( IS, JS ), LDC,
1059      $                           ONE, C( IE+1, JS ), LDC )
1060                      CALL SGEMM( 'T', 'N', M-IE, NB, MB, -ONE,
1061      $                           D( IS, IE+1 ), LDD, F( IS, JS ), LDF,
1062      $                           ONE, C( IE+1, JS ), LDC )
1063                   END IF
1064 *
1065                END IF
1066 *
1067   190       CONTINUE
1068   200    CONTINUE
1069 *
1070       END IF
1071       RETURN
1072 *
1073 *     End of STGSY2
1074 *
1075       END