LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ stgsy2()

subroutine stgsy2 ( character  TRANS,
integer  IJOB,
integer  M,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( ldb, * )  B,
integer  LDB,
real, dimension( ldc, * )  C,
integer  LDC,
real, dimension( ldd, * )  D,
integer  LDD,
real, dimension( lde, * )  E,
integer  LDE,
real, dimension( ldf, * )  F,
integer  LDF,
real  SCALE,
real  RDSUM,
real  RDSCAL,
integer, dimension( * )  IWORK,
integer  PQ,
integer  INFO 
)

STGSY2 solves the generalized Sylvester equation (unblocked algorithm).

Download STGSY2 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 STGSY2 solves the generalized Sylvester equation:

             A * R - L * B = scale * C                (1)
             D * R - L * E = scale * F,

 using Level 1 and 2 BLAS. where R and L are unknown M-by-N matrices,
 (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M,
 N-by-N and M-by-N, respectively, with real entries. (A, D) and (B, E)
 must be in generalized Schur canonical form, i.e. A, B are upper
 quasi triangular and D, E are upper triangular. The solution (R, L)
 overwrites (C, F). 0 <= SCALE <= 1 is an output scaling factor
 chosen to avoid overflow.

 In matrix notation solving equation (1) corresponds to solve
 Z*x = scale*b, where Z is defined as

        Z = [ kron(In, A)  -kron(B**T, Im) ]             (2)
            [ kron(In, D)  -kron(E**T, Im) ],

 Ik is the identity matrix of size k and X**T is the transpose of X.
 kron(X, Y) is the Kronecker product between the matrices X and Y.
 In the process of solving (1), we solve a number of such systems
 where Dim(In), Dim(In) = 1 or 2.

 If TRANS = 'T', solve the transposed system Z**T*y = scale*b for y,
 which is equivalent to solve for R and L in

             A**T * R  + D**T * L   = scale * C           (3)
             R  * B**T + L  * E**T  = scale * -F

 This case is used to compute an estimate of Dif[(A, D), (B, E)] =
 sigma_min(Z) using reverse communication with SLACON.

 STGSY2 also (IJOB >= 1) contributes to the computation in STGSYL
 of an upper bound on the separation between to matrix pairs. Then
 the input (A, D), (B, E) are sub-pencils of the matrix pair in
 STGSYL. See STGSYL for details.
Parameters
[in]TRANS
          TRANS is CHARACTER*1
          = 'N': solve the generalized Sylvester equation (1).
          = 'T': solve the 'transposed' system (3).
[in]IJOB
          IJOB is INTEGER
          Specifies what kind of functionality to be performed.
          = 0: solve (1) only.
          = 1: A contribution from this subsystem to a Frobenius
               norm-based estimate of the separation between two matrix
               pairs is computed. (look ahead strategy is used).
          = 2: A contribution from this subsystem to a Frobenius
               norm-based estimate of the separation between two matrix
               pairs is computed. (SGECON on sub-systems is used.)
          Not referenced if TRANS = 'T'.
[in]M
          M is INTEGER
          On entry, M specifies the order of A and D, and the row
          dimension of C, F, R and L.
[in]N
          N is INTEGER
          On entry, N specifies the order of B and E, and the column
          dimension of C, F, R and L.
[in]A
          A is REAL array, dimension (LDA, M)
          On entry, A contains an upper quasi triangular matrix.
[in]LDA
          LDA is INTEGER
          The leading dimension of the matrix A. LDA >= max(1, M).
[in]B
          B is REAL array, dimension (LDB, N)
          On entry, B contains an upper quasi triangular matrix.
[in]LDB
          LDB is INTEGER
          The leading dimension of the matrix B. LDB >= max(1, N).
[in,out]C
          C is REAL array, dimension (LDC, N)
          On entry, C contains the right-hand-side of the first matrix
          equation in (1).
          On exit, if IJOB = 0, C has been overwritten by the
          solution R.
[in]LDC
          LDC is INTEGER
          The leading dimension of the matrix C. LDC >= max(1, M).
[in]D
          D is REAL array, dimension (LDD, M)
          On entry, D contains an upper triangular matrix.
[in]LDD
          LDD is INTEGER
          The leading dimension of the matrix D. LDD >= max(1, M).
[in]E
          E is REAL array, dimension (LDE, N)
          On entry, E contains an upper triangular matrix.
[in]LDE
          LDE is INTEGER
          The leading dimension of the matrix E. LDE >= max(1, N).
[in,out]F
          F is REAL array, dimension (LDF, N)
          On entry, F contains the right-hand-side of the second matrix
          equation in (1).
          On exit, if IJOB = 0, F has been overwritten by the
          solution L.
[in]LDF
          LDF is INTEGER
          The leading dimension of the matrix F. LDF >= max(1, M).
[out]SCALE
          SCALE is REAL
          On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions
          R and L (C and F on entry) will hold the solutions to a
          slightly perturbed system but the input matrices A, B, D and
          E have not been changed. If SCALE = 0, R and L will hold the
          solutions to the homogeneous system with C = F = 0. Normally,
          SCALE = 1.
[in,out]RDSUM
          RDSUM is REAL
          On entry, the sum of squares of computed contributions to
          the Dif-estimate under computation by STGSYL, where the
          scaling factor RDSCAL (see below) has been factored out.
          On exit, the corresponding sum of squares updated with the
          contributions from the current sub-system.
          If TRANS = 'T' RDSUM is not touched.
          NOTE: RDSUM only makes sense when STGSY2 is called by STGSYL.
[in,out]RDSCAL
          RDSCAL is REAL
          On entry, scaling factor used to prevent overflow in RDSUM.
          On exit, RDSCAL is updated w.r.t. the current contributions
          in RDSUM.
          If TRANS = 'T', RDSCAL is not touched.
          NOTE: RDSCAL only makes sense when STGSY2 is called by
                STGSYL.
[out]IWORK
          IWORK is INTEGER array, dimension (M+N+2)
[out]PQ
          PQ is INTEGER
          On exit, the number of subsystems (of size 2-by-2, 4-by-4 and
          8-by-8) solved by this routine.
[out]INFO
          INFO is INTEGER
          On exit, if INFO is set to
            =0: Successful exit
            <0: If INFO = -i, the i-th argument had an illegal value.
            >0: The matrix pairs (A, D) and (B, E) have common or very
                close eigenvalues.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.

Definition at line 271 of file stgsy2.f.

274 *
275 * -- LAPACK auxiliary routine --
276 * -- LAPACK is a software package provided by Univ. of Tennessee, --
277 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
278 *
279 * .. Scalar Arguments ..
280  CHARACTER TRANS
281  INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N,
282  $ PQ
283  REAL RDSCAL, RDSUM, SCALE
284 * ..
285 * .. Array Arguments ..
286  INTEGER IWORK( * )
287  REAL A( LDA, * ), B( LDB, * ), C( LDC, * ),
288  $ D( LDD, * ), E( LDE, * ), F( LDF, * )
289 * ..
290 *
291 * =====================================================================
292 * Replaced various illegal calls to SCOPY by calls to SLASET.
293 * Sven Hammarling, 27/5/02.
294 *
295 * .. Parameters ..
296  INTEGER LDZ
297  parameter( ldz = 8 )
298  REAL ZERO, ONE
299  parameter( zero = 0.0e+0, one = 1.0e+0 )
300 * ..
301 * .. Local Scalars ..
302  LOGICAL NOTRAN
303  INTEGER I, IE, IERR, II, IS, ISP1, J, JE, JJ, JS, JSP1,
304  $ K, MB, NB, P, Q, ZDIM
305  REAL ALPHA, SCALOC
306 * ..
307 * .. Local Arrays ..
308  INTEGER IPIV( LDZ ), JPIV( LDZ )
309  REAL RHS( LDZ ), Z( LDZ, LDZ )
310 * ..
311 * .. External Functions ..
312  LOGICAL LSAME
313  EXTERNAL lsame
314 * ..
315 * .. External Subroutines ..
316  EXTERNAL saxpy, scopy, sgemm, sgemv, sger, sgesc2,
318 * ..
319 * .. Intrinsic Functions ..
320  INTRINSIC max
321 * ..
322 * .. Executable Statements ..
323 *
324 * Decode and test input parameters
325 *
326  info = 0
327  ierr = 0
328  notran = lsame( trans, 'N' )
329  IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) ) THEN
330  info = -1
331  ELSE IF( notran ) THEN
332  IF( ( ijob.LT.0 ) .OR. ( ijob.GT.2 ) ) THEN
333  info = -2
334  END IF
335  END IF
336  IF( info.EQ.0 ) THEN
337  IF( m.LE.0 ) THEN
338  info = -3
339  ELSE IF( n.LE.0 ) THEN
340  info = -4
341  ELSE IF( lda.LT.max( 1, m ) ) THEN
342  info = -6
343  ELSE IF( ldb.LT.max( 1, n ) ) THEN
344  info = -8
345  ELSE IF( ldc.LT.max( 1, m ) ) THEN
346  info = -10
347  ELSE IF( ldd.LT.max( 1, m ) ) THEN
348  info = -12
349  ELSE IF( lde.LT.max( 1, n ) ) THEN
350  info = -14
351  ELSE IF( ldf.LT.max( 1, m ) ) THEN
352  info = -16
353  END IF
354  END IF
355  IF( info.NE.0 ) THEN
356  CALL xerbla( 'STGSY2', -info )
357  RETURN
358  END IF
359 *
360 * Determine block structure of A
361 *
362  pq = 0
363  p = 0
364  i = 1
365  10 CONTINUE
366  IF( i.GT.m )
367  $ GO TO 20
368  p = p + 1
369  iwork( p ) = i
370  IF( i.EQ.m )
371  $ GO TO 20
372  IF( a( i+1, i ).NE.zero ) THEN
373  i = i + 2
374  ELSE
375  i = i + 1
376  END IF
377  GO TO 10
378  20 CONTINUE
379  iwork( p+1 ) = m + 1
380 *
381 * Determine block structure of B
382 *
383  q = p + 1
384  j = 1
385  30 CONTINUE
386  IF( j.GT.n )
387  $ GO TO 40
388  q = q + 1
389  iwork( q ) = j
390  IF( j.EQ.n )
391  $ GO TO 40
392  IF( b( j+1, j ).NE.zero ) THEN
393  j = j + 2
394  ELSE
395  j = j + 1
396  END IF
397  GO TO 30
398  40 CONTINUE
399  iwork( q+1 ) = n + 1
400  pq = p*( q-p-1 )
401 *
402  IF( notran ) THEN
403 *
404 * Solve (I, J) - subsystem
405 * A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
406 * D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
407 * for I = P, P - 1, ..., 1; J = 1, 2, ..., Q
408 *
409  scale = one
410  scaloc = one
411  DO 120 j = p + 2, q
412  js = iwork( j )
413  jsp1 = js + 1
414  je = iwork( j+1 ) - 1
415  nb = je - js + 1
416  DO 110 i = p, 1, -1
417 *
418  is = iwork( i )
419  isp1 = is + 1
420  ie = iwork( i+1 ) - 1
421  mb = ie - is + 1
422  zdim = mb*nb*2
423 *
424  IF( ( mb.EQ.1 ) .AND. ( nb.EQ.1 ) ) THEN
425 *
426 * Build a 2-by-2 system Z * x = RHS
427 *
428  z( 1, 1 ) = a( is, is )
429  z( 2, 1 ) = d( is, is )
430  z( 1, 2 ) = -b( js, js )
431  z( 2, 2 ) = -e( js, js )
432 *
433 * Set up right hand side(s)
434 *
435  rhs( 1 ) = c( is, js )
436  rhs( 2 ) = f( is, js )
437 *
438 * Solve Z * x = RHS
439 *
440  CALL sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
441  IF( ierr.GT.0 )
442  $ info = ierr
443 *
444  IF( ijob.EQ.0 ) THEN
445  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
446  $ scaloc )
447  IF( scaloc.NE.one ) THEN
448  DO 50 k = 1, n
449  CALL sscal( m, scaloc, c( 1, k ), 1 )
450  CALL sscal( m, scaloc, f( 1, k ), 1 )
451  50 CONTINUE
452  scale = scale*scaloc
453  END IF
454  ELSE
455  CALL slatdf( ijob, zdim, z, ldz, rhs, rdsum,
456  $ rdscal, ipiv, jpiv )
457  END IF
458 *
459 * Unpack solution vector(s)
460 *
461  c( is, js ) = rhs( 1 )
462  f( is, js ) = rhs( 2 )
463 *
464 * Substitute R(I, J) and L(I, J) into remaining
465 * equation.
466 *
467  IF( i.GT.1 ) THEN
468  alpha = -rhs( 1 )
469  CALL saxpy( is-1, alpha, a( 1, is ), 1, c( 1, js ),
470  $ 1 )
471  CALL saxpy( is-1, alpha, d( 1, is ), 1, f( 1, js ),
472  $ 1 )
473  END IF
474  IF( j.LT.q ) THEN
475  CALL saxpy( n-je, rhs( 2 ), b( js, je+1 ), ldb,
476  $ c( is, je+1 ), ldc )
477  CALL saxpy( n-je, rhs( 2 ), e( js, je+1 ), lde,
478  $ f( is, je+1 ), ldf )
479  END IF
480 *
481  ELSE IF( ( mb.EQ.1 ) .AND. ( nb.EQ.2 ) ) THEN
482 *
483 * Build a 4-by-4 system Z * x = RHS
484 *
485  z( 1, 1 ) = a( is, is )
486  z( 2, 1 ) = zero
487  z( 3, 1 ) = d( is, is )
488  z( 4, 1 ) = zero
489 *
490  z( 1, 2 ) = zero
491  z( 2, 2 ) = a( is, is )
492  z( 3, 2 ) = zero
493  z( 4, 2 ) = d( is, is )
494 *
495  z( 1, 3 ) = -b( js, js )
496  z( 2, 3 ) = -b( js, jsp1 )
497  z( 3, 3 ) = -e( js, js )
498  z( 4, 3 ) = -e( js, jsp1 )
499 *
500  z( 1, 4 ) = -b( jsp1, js )
501  z( 2, 4 ) = -b( jsp1, jsp1 )
502  z( 3, 4 ) = zero
503  z( 4, 4 ) = -e( jsp1, jsp1 )
504 *
505 * Set up right hand side(s)
506 *
507  rhs( 1 ) = c( is, js )
508  rhs( 2 ) = c( is, jsp1 )
509  rhs( 3 ) = f( is, js )
510  rhs( 4 ) = f( is, jsp1 )
511 *
512 * Solve Z * x = RHS
513 *
514  CALL sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
515  IF( ierr.GT.0 )
516  $ info = ierr
517 *
518  IF( ijob.EQ.0 ) THEN
519  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
520  $ scaloc )
521  IF( scaloc.NE.one ) THEN
522  DO 60 k = 1, n
523  CALL sscal( m, scaloc, c( 1, k ), 1 )
524  CALL sscal( m, scaloc, f( 1, k ), 1 )
525  60 CONTINUE
526  scale = scale*scaloc
527  END IF
528  ELSE
529  CALL slatdf( ijob, zdim, z, ldz, rhs, rdsum,
530  $ rdscal, ipiv, jpiv )
531  END IF
532 *
533 * Unpack solution vector(s)
534 *
535  c( is, js ) = rhs( 1 )
536  c( is, jsp1 ) = rhs( 2 )
537  f( is, js ) = rhs( 3 )
538  f( is, jsp1 ) = rhs( 4 )
539 *
540 * Substitute R(I, J) and L(I, J) into remaining
541 * equation.
542 *
543  IF( i.GT.1 ) THEN
544  CALL sger( is-1, nb, -one, a( 1, is ), 1, rhs( 1 ),
545  $ 1, c( 1, js ), ldc )
546  CALL sger( is-1, nb, -one, d( 1, is ), 1, rhs( 1 ),
547  $ 1, f( 1, js ), ldf )
548  END IF
549  IF( j.LT.q ) THEN
550  CALL saxpy( n-je, rhs( 3 ), b( js, je+1 ), ldb,
551  $ c( is, je+1 ), ldc )
552  CALL saxpy( n-je, rhs( 3 ), e( js, je+1 ), lde,
553  $ f( is, je+1 ), ldf )
554  CALL saxpy( n-je, rhs( 4 ), b( jsp1, je+1 ), ldb,
555  $ c( is, je+1 ), ldc )
556  CALL saxpy( n-je, rhs( 4 ), e( jsp1, je+1 ), lde,
557  $ f( is, je+1 ), ldf )
558  END IF
559 *
560  ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.1 ) ) THEN
561 *
562 * Build a 4-by-4 system Z * x = RHS
563 *
564  z( 1, 1 ) = a( is, is )
565  z( 2, 1 ) = a( isp1, is )
566  z( 3, 1 ) = d( is, is )
567  z( 4, 1 ) = zero
568 *
569  z( 1, 2 ) = a( is, isp1 )
570  z( 2, 2 ) = a( isp1, isp1 )
571  z( 3, 2 ) = d( is, isp1 )
572  z( 4, 2 ) = d( isp1, isp1 )
573 *
574  z( 1, 3 ) = -b( js, js )
575  z( 2, 3 ) = zero
576  z( 3, 3 ) = -e( js, js )
577  z( 4, 3 ) = zero
578 *
579  z( 1, 4 ) = zero
580  z( 2, 4 ) = -b( js, js )
581  z( 3, 4 ) = zero
582  z( 4, 4 ) = -e( js, js )
583 *
584 * Set up right hand side(s)
585 *
586  rhs( 1 ) = c( is, js )
587  rhs( 2 ) = c( isp1, js )
588  rhs( 3 ) = f( is, js )
589  rhs( 4 ) = f( isp1, js )
590 *
591 * Solve Z * x = RHS
592 *
593  CALL sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
594  IF( ierr.GT.0 )
595  $ info = ierr
596  IF( ijob.EQ.0 ) THEN
597  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
598  $ scaloc )
599  IF( scaloc.NE.one ) THEN
600  DO 70 k = 1, n
601  CALL sscal( m, scaloc, c( 1, k ), 1 )
602  CALL sscal( m, scaloc, f( 1, k ), 1 )
603  70 CONTINUE
604  scale = scale*scaloc
605  END IF
606  ELSE
607  CALL slatdf( ijob, zdim, z, ldz, rhs, rdsum,
608  $ rdscal, ipiv, jpiv )
609  END IF
610 *
611 * Unpack solution vector(s)
612 *
613  c( is, js ) = rhs( 1 )
614  c( isp1, js ) = rhs( 2 )
615  f( is, js ) = rhs( 3 )
616  f( isp1, js ) = rhs( 4 )
617 *
618 * Substitute R(I, J) and L(I, J) into remaining
619 * equation.
620 *
621  IF( i.GT.1 ) THEN
622  CALL sgemv( 'N', is-1, mb, -one, a( 1, is ), lda,
623  $ rhs( 1 ), 1, one, c( 1, js ), 1 )
624  CALL sgemv( 'N', is-1, mb, -one, d( 1, is ), ldd,
625  $ rhs( 1 ), 1, one, f( 1, js ), 1 )
626  END IF
627  IF( j.LT.q ) THEN
628  CALL sger( mb, n-je, one, rhs( 3 ), 1,
629  $ b( js, je+1 ), ldb, c( is, je+1 ), ldc )
630  CALL sger( mb, n-je, one, rhs( 3 ), 1,
631  $ e( js, je+1 ), lde, f( is, je+1 ), ldf )
632  END IF
633 *
634  ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.2 ) ) THEN
635 *
636 * Build an 8-by-8 system Z * x = RHS
637 *
638  CALL slaset( 'F', ldz, ldz, zero, zero, z, ldz )
639 *
640  z( 1, 1 ) = a( is, is )
641  z( 2, 1 ) = a( isp1, is )
642  z( 5, 1 ) = d( is, is )
643 *
644  z( 1, 2 ) = a( is, isp1 )
645  z( 2, 2 ) = a( isp1, isp1 )
646  z( 5, 2 ) = d( is, isp1 )
647  z( 6, 2 ) = d( isp1, isp1 )
648 *
649  z( 3, 3 ) = a( is, is )
650  z( 4, 3 ) = a( isp1, is )
651  z( 7, 3 ) = d( is, is )
652 *
653  z( 3, 4 ) = a( is, isp1 )
654  z( 4, 4 ) = a( isp1, isp1 )
655  z( 7, 4 ) = d( is, isp1 )
656  z( 8, 4 ) = d( isp1, isp1 )
657 *
658  z( 1, 5 ) = -b( js, js )
659  z( 3, 5 ) = -b( js, jsp1 )
660  z( 5, 5 ) = -e( js, js )
661  z( 7, 5 ) = -e( js, jsp1 )
662 *
663  z( 2, 6 ) = -b( js, js )
664  z( 4, 6 ) = -b( js, jsp1 )
665  z( 6, 6 ) = -e( js, js )
666  z( 8, 6 ) = -e( js, jsp1 )
667 *
668  z( 1, 7 ) = -b( jsp1, js )
669  z( 3, 7 ) = -b( jsp1, jsp1 )
670  z( 7, 7 ) = -e( jsp1, jsp1 )
671 *
672  z( 2, 8 ) = -b( jsp1, js )
673  z( 4, 8 ) = -b( jsp1, jsp1 )
674  z( 8, 8 ) = -e( jsp1, jsp1 )
675 *
676 * Set up right hand side(s)
677 *
678  k = 1
679  ii = mb*nb + 1
680  DO 80 jj = 0, nb - 1
681  CALL scopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
682  CALL scopy( mb, f( is, js+jj ), 1, rhs( ii ), 1 )
683  k = k + mb
684  ii = ii + mb
685  80 CONTINUE
686 *
687 * Solve Z * x = RHS
688 *
689  CALL sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
690  IF( ierr.GT.0 )
691  $ info = ierr
692  IF( ijob.EQ.0 ) THEN
693  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
694  $ scaloc )
695  IF( scaloc.NE.one ) THEN
696  DO 90 k = 1, n
697  CALL sscal( m, scaloc, c( 1, k ), 1 )
698  CALL sscal( m, scaloc, f( 1, k ), 1 )
699  90 CONTINUE
700  scale = scale*scaloc
701  END IF
702  ELSE
703  CALL slatdf( ijob, zdim, z, ldz, rhs, rdsum,
704  $ rdscal, ipiv, jpiv )
705  END IF
706 *
707 * Unpack solution vector(s)
708 *
709  k = 1
710  ii = mb*nb + 1
711  DO 100 jj = 0, nb - 1
712  CALL scopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
713  CALL scopy( mb, rhs( ii ), 1, f( is, js+jj ), 1 )
714  k = k + mb
715  ii = ii + mb
716  100 CONTINUE
717 *
718 * Substitute R(I, J) and L(I, J) into remaining
719 * equation.
720 *
721  IF( i.GT.1 ) THEN
722  CALL sgemm( 'N', 'N', is-1, nb, mb, -one,
723  $ a( 1, is ), lda, rhs( 1 ), mb, one,
724  $ c( 1, js ), ldc )
725  CALL sgemm( 'N', 'N', is-1, nb, mb, -one,
726  $ d( 1, is ), ldd, rhs( 1 ), mb, one,
727  $ f( 1, js ), ldf )
728  END IF
729  IF( j.LT.q ) THEN
730  k = mb*nb + 1
731  CALL sgemm( 'N', 'N', mb, n-je, nb, one, rhs( k ),
732  $ mb, b( js, je+1 ), ldb, one,
733  $ c( is, je+1 ), ldc )
734  CALL sgemm( 'N', 'N', mb, n-je, nb, one, rhs( k ),
735  $ mb, e( js, je+1 ), lde, one,
736  $ f( is, je+1 ), ldf )
737  END IF
738 *
739  END IF
740 *
741  110 CONTINUE
742  120 CONTINUE
743  ELSE
744 *
745 * Solve (I, J) - subsystem
746 * A(I, I)**T * R(I, J) + D(I, I)**T * L(J, J) = C(I, J)
747 * R(I, I) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
748 * for I = 1, 2, ..., P, J = Q, Q - 1, ..., 1
749 *
750  scale = one
751  scaloc = one
752  DO 200 i = 1, p
753 *
754  is = iwork( i )
755  isp1 = is + 1
756  ie = iwork( i+1 ) - 1
757  mb = ie - is + 1
758  DO 190 j = q, p + 2, -1
759 *
760  js = iwork( j )
761  jsp1 = js + 1
762  je = iwork( j+1 ) - 1
763  nb = je - js + 1
764  zdim = mb*nb*2
765  IF( ( mb.EQ.1 ) .AND. ( nb.EQ.1 ) ) THEN
766 *
767 * Build a 2-by-2 system Z**T * x = RHS
768 *
769  z( 1, 1 ) = a( is, is )
770  z( 2, 1 ) = -b( js, js )
771  z( 1, 2 ) = d( is, is )
772  z( 2, 2 ) = -e( js, js )
773 *
774 * Set up right hand side(s)
775 *
776  rhs( 1 ) = c( is, js )
777  rhs( 2 ) = f( is, js )
778 *
779 * Solve Z**T * x = RHS
780 *
781  CALL sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
782  IF( ierr.GT.0 )
783  $ info = ierr
784 *
785  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
786  IF( scaloc.NE.one ) THEN
787  DO 130 k = 1, n
788  CALL sscal( m, scaloc, c( 1, k ), 1 )
789  CALL sscal( m, scaloc, f( 1, k ), 1 )
790  130 CONTINUE
791  scale = scale*scaloc
792  END IF
793 *
794 * Unpack solution vector(s)
795 *
796  c( is, js ) = rhs( 1 )
797  f( is, js ) = rhs( 2 )
798 *
799 * Substitute R(I, J) and L(I, J) into remaining
800 * equation.
801 *
802  IF( j.GT.p+2 ) THEN
803  alpha = rhs( 1 )
804  CALL saxpy( js-1, alpha, b( 1, js ), 1, f( is, 1 ),
805  $ ldf )
806  alpha = rhs( 2 )
807  CALL saxpy( js-1, alpha, e( 1, js ), 1, f( is, 1 ),
808  $ ldf )
809  END IF
810  IF( i.LT.p ) THEN
811  alpha = -rhs( 1 )
812  CALL saxpy( m-ie, alpha, a( is, ie+1 ), lda,
813  $ c( ie+1, js ), 1 )
814  alpha = -rhs( 2 )
815  CALL saxpy( m-ie, alpha, d( is, ie+1 ), ldd,
816  $ c( ie+1, js ), 1 )
817  END IF
818 *
819  ELSE IF( ( mb.EQ.1 ) .AND. ( nb.EQ.2 ) ) THEN
820 *
821 * Build a 4-by-4 system Z**T * x = RHS
822 *
823  z( 1, 1 ) = a( is, is )
824  z( 2, 1 ) = zero
825  z( 3, 1 ) = -b( js, js )
826  z( 4, 1 ) = -b( jsp1, js )
827 *
828  z( 1, 2 ) = zero
829  z( 2, 2 ) = a( is, is )
830  z( 3, 2 ) = -b( js, jsp1 )
831  z( 4, 2 ) = -b( jsp1, jsp1 )
832 *
833  z( 1, 3 ) = d( is, is )
834  z( 2, 3 ) = zero
835  z( 3, 3 ) = -e( js, js )
836  z( 4, 3 ) = zero
837 *
838  z( 1, 4 ) = zero
839  z( 2, 4 ) = d( is, is )
840  z( 3, 4 ) = -e( js, jsp1 )
841  z( 4, 4 ) = -e( jsp1, jsp1 )
842 *
843 * Set up right hand side(s)
844 *
845  rhs( 1 ) = c( is, js )
846  rhs( 2 ) = c( is, jsp1 )
847  rhs( 3 ) = f( is, js )
848  rhs( 4 ) = f( is, jsp1 )
849 *
850 * Solve Z**T * x = RHS
851 *
852  CALL sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
853  IF( ierr.GT.0 )
854  $ info = ierr
855  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
856  IF( scaloc.NE.one ) THEN
857  DO 140 k = 1, n
858  CALL sscal( m, scaloc, c( 1, k ), 1 )
859  CALL sscal( m, scaloc, f( 1, k ), 1 )
860  140 CONTINUE
861  scale = scale*scaloc
862  END IF
863 *
864 * Unpack solution vector(s)
865 *
866  c( is, js ) = rhs( 1 )
867  c( is, jsp1 ) = rhs( 2 )
868  f( is, js ) = rhs( 3 )
869  f( is, jsp1 ) = rhs( 4 )
870 *
871 * Substitute R(I, J) and L(I, J) into remaining
872 * equation.
873 *
874  IF( j.GT.p+2 ) THEN
875  CALL saxpy( js-1, rhs( 1 ), b( 1, js ), 1,
876  $ f( is, 1 ), ldf )
877  CALL saxpy( js-1, rhs( 2 ), b( 1, jsp1 ), 1,
878  $ f( is, 1 ), ldf )
879  CALL saxpy( js-1, rhs( 3 ), e( 1, js ), 1,
880  $ f( is, 1 ), ldf )
881  CALL saxpy( js-1, rhs( 4 ), e( 1, jsp1 ), 1,
882  $ f( is, 1 ), ldf )
883  END IF
884  IF( i.LT.p ) THEN
885  CALL sger( m-ie, nb, -one, a( is, ie+1 ), lda,
886  $ rhs( 1 ), 1, c( ie+1, js ), ldc )
887  CALL sger( m-ie, nb, -one, d( is, ie+1 ), ldd,
888  $ rhs( 3 ), 1, c( ie+1, js ), ldc )
889  END IF
890 *
891  ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.1 ) ) THEN
892 *
893 * Build a 4-by-4 system Z**T * x = RHS
894 *
895  z( 1, 1 ) = a( is, is )
896  z( 2, 1 ) = a( is, isp1 )
897  z( 3, 1 ) = -b( js, js )
898  z( 4, 1 ) = zero
899 *
900  z( 1, 2 ) = a( isp1, is )
901  z( 2, 2 ) = a( isp1, isp1 )
902  z( 3, 2 ) = zero
903  z( 4, 2 ) = -b( js, js )
904 *
905  z( 1, 3 ) = d( is, is )
906  z( 2, 3 ) = d( is, isp1 )
907  z( 3, 3 ) = -e( js, js )
908  z( 4, 3 ) = zero
909 *
910  z( 1, 4 ) = zero
911  z( 2, 4 ) = d( isp1, isp1 )
912  z( 3, 4 ) = zero
913  z( 4, 4 ) = -e( js, js )
914 *
915 * Set up right hand side(s)
916 *
917  rhs( 1 ) = c( is, js )
918  rhs( 2 ) = c( isp1, js )
919  rhs( 3 ) = f( is, js )
920  rhs( 4 ) = f( isp1, js )
921 *
922 * Solve Z**T * x = RHS
923 *
924  CALL sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
925  IF( ierr.GT.0 )
926  $ info = ierr
927 *
928  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
929  IF( scaloc.NE.one ) THEN
930  DO 150 k = 1, n
931  CALL sscal( m, scaloc, c( 1, k ), 1 )
932  CALL sscal( m, scaloc, f( 1, k ), 1 )
933  150 CONTINUE
934  scale = scale*scaloc
935  END IF
936 *
937 * Unpack solution vector(s)
938 *
939  c( is, js ) = rhs( 1 )
940  c( isp1, js ) = rhs( 2 )
941  f( is, js ) = rhs( 3 )
942  f( isp1, js ) = rhs( 4 )
943 *
944 * Substitute R(I, J) and L(I, J) into remaining
945 * equation.
946 *
947  IF( j.GT.p+2 ) THEN
948  CALL sger( mb, js-1, one, rhs( 1 ), 1, b( 1, js ),
949  $ 1, f( is, 1 ), ldf )
950  CALL sger( mb, js-1, one, rhs( 3 ), 1, e( 1, js ),
951  $ 1, f( is, 1 ), ldf )
952  END IF
953  IF( i.LT.p ) THEN
954  CALL sgemv( 'T', mb, m-ie, -one, a( is, ie+1 ),
955  $ lda, rhs( 1 ), 1, one, c( ie+1, js ),
956  $ 1 )
957  CALL sgemv( 'T', mb, m-ie, -one, d( is, ie+1 ),
958  $ ldd, rhs( 3 ), 1, one, c( ie+1, js ),
959  $ 1 )
960  END IF
961 *
962  ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.2 ) ) THEN
963 *
964 * Build an 8-by-8 system Z**T * x = RHS
965 *
966  CALL slaset( 'F', ldz, ldz, zero, zero, z, ldz )
967 *
968  z( 1, 1 ) = a( is, is )
969  z( 2, 1 ) = a( is, isp1 )
970  z( 5, 1 ) = -b( js, js )
971  z( 7, 1 ) = -b( jsp1, js )
972 *
973  z( 1, 2 ) = a( isp1, is )
974  z( 2, 2 ) = a( isp1, isp1 )
975  z( 6, 2 ) = -b( js, js )
976  z( 8, 2 ) = -b( jsp1, js )
977 *
978  z( 3, 3 ) = a( is, is )
979  z( 4, 3 ) = a( is, isp1 )
980  z( 5, 3 ) = -b( js, jsp1 )
981  z( 7, 3 ) = -b( jsp1, jsp1 )
982 *
983  z( 3, 4 ) = a( isp1, is )
984  z( 4, 4 ) = a( isp1, isp1 )
985  z( 6, 4 ) = -b( js, jsp1 )
986  z( 8, 4 ) = -b( jsp1, jsp1 )
987 *
988  z( 1, 5 ) = d( is, is )
989  z( 2, 5 ) = d( is, isp1 )
990  z( 5, 5 ) = -e( js, js )
991 *
992  z( 2, 6 ) = d( isp1, isp1 )
993  z( 6, 6 ) = -e( js, js )
994 *
995  z( 3, 7 ) = d( is, is )
996  z( 4, 7 ) = d( is, isp1 )
997  z( 5, 7 ) = -e( js, jsp1 )
998  z( 7, 7 ) = -e( jsp1, jsp1 )
999 *
1000  z( 4, 8 ) = d( isp1, isp1 )
1001  z( 6, 8 ) = -e( js, jsp1 )
1002  z( 8, 8 ) = -e( jsp1, jsp1 )
1003 *
1004 * Set up right hand side(s)
1005 *
1006  k = 1
1007  ii = mb*nb + 1
1008  DO 160 jj = 0, nb - 1
1009  CALL scopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
1010  CALL scopy( mb, f( is, js+jj ), 1, rhs( ii ), 1 )
1011  k = k + mb
1012  ii = ii + mb
1013  160 CONTINUE
1014 *
1015 *
1016 * Solve Z**T * x = RHS
1017 *
1018  CALL sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
1019  IF( ierr.GT.0 )
1020  $ info = ierr
1021 *
1022  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
1023  IF( scaloc.NE.one ) THEN
1024  DO 170 k = 1, n
1025  CALL sscal( m, scaloc, c( 1, k ), 1 )
1026  CALL sscal( m, scaloc, f( 1, k ), 1 )
1027  170 CONTINUE
1028  scale = scale*scaloc
1029  END IF
1030 *
1031 * Unpack solution vector(s)
1032 *
1033  k = 1
1034  ii = mb*nb + 1
1035  DO 180 jj = 0, nb - 1
1036  CALL scopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
1037  CALL scopy( mb, rhs( ii ), 1, f( is, js+jj ), 1 )
1038  k = k + mb
1039  ii = ii + mb
1040  180 CONTINUE
1041 *
1042 * Substitute R(I, J) and L(I, J) into remaining
1043 * equation.
1044 *
1045  IF( j.GT.p+2 ) THEN
1046  CALL sgemm( 'N', 'T', mb, js-1, nb, one,
1047  $ c( is, js ), ldc, b( 1, js ), ldb, one,
1048  $ f( is, 1 ), ldf )
1049  CALL sgemm( 'N', 'T', mb, js-1, nb, one,
1050  $ f( is, js ), ldf, e( 1, js ), lde, one,
1051  $ f( is, 1 ), ldf )
1052  END IF
1053  IF( i.LT.p ) THEN
1054  CALL sgemm( 'T', 'N', m-ie, nb, mb, -one,
1055  $ a( is, ie+1 ), lda, c( is, js ), ldc,
1056  $ one, c( ie+1, js ), ldc )
1057  CALL sgemm( 'T', 'N', m-ie, nb, mb, -one,
1058  $ d( is, ie+1 ), ldd, f( is, js ), ldf,
1059  $ one, c( ie+1, js ), ldc )
1060  END IF
1061 *
1062  END IF
1063 *
1064  190 CONTINUE
1065  200 CONTINUE
1066 *
1067  END IF
1068  RETURN
1069 *
1070 * End of STGSY2
1071 *
logical function lde(RI, RJ, LR)
Definition: dblat2.f:2942
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.f:110
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine sgesc2(N, A, LDA, RHS, IPIV, JPIV, SCALE)
SGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed...
Definition: sgesc2.f:114
subroutine sgetc2(N, A, LDA, IPIV, JPIV, INFO)
SGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix.
Definition: sgetc2.f:111
subroutine slatdf(IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV, JPIV)
SLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution ...
Definition: slatdf.f:171
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:89
subroutine sger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
SGER
Definition: sger.f:130
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:156
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
Here is the call graph for this function:
Here is the caller graph for this function: