LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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 communicaton 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.
Date
November 2015
Contributors:
Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.

Definition at line 276 of file stgsy2.f.

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,
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 *
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:173
subroutine sger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
SGER
Definition: sger.f:132
logical function lde(RI, RJ, LR)
Definition: dblat2.f:2945
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
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:113
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:158
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:112
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:116
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:54
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53

Here is the call graph for this function:

Here is the caller graph for this function: