dd/d6a/bslaexc_8f_source.html

      SUBROUTINE bslaexc( N, T, LDT, J1, N1, N2, ITRAF, DTRAF, WORK,

     $                    INFO )

      IMPLICIT NONE

*

*  -- ScaLAPACK auxiliary routine (version 2.0.2) --

*     Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver

*     May 1 2012

*

*     .. Scalar Arguments ..

      INTEGER            INFO, J1, LDT, N, N1, N2

*     ..

*     .. Array Arguments ..

      INTEGER            ITRAF( * )

      REAL               DTRAF( * ), T( LDT, * ), WORK( * )

*     ..

*

*  Purpose

*  =======

*

*  BSLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in

*  an upper quasi-triangular matrix T by an orthogonal similarity

*  transformation.

*

*  In contrast to the LAPACK routine DLAEXC, the orthogonal

*  transformation matrix Q is not explicitly constructed but

*  represented by paramaters contained in the arrays ITRAF and DTRAF,

*  see the description of BSTREXC for more details.

*

*  T must be in Schur canonical form, that is, block upper triangular

*  with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block

*  has its diagonal elemnts equal and its off-diagonal elements of

*  opposite sign.

*

*  Arguments

*  =========

*

*  N       (input) INTEGER

*          The order of the matrix T. N >= 0.

*

*  T       (input/output) REAL             array, dimension (LDT,N)

*          On entry, the upper quasi-triangular matrix T, in Schur

*          canonical form.

*          On exit, the updated matrix T, again in Schur canonical form.

*

*  LDT     (input)  INTEGER

*          The leading dimension of the array T. LDT >= max(1,N).

*

*  J1      (input) INTEGER

*          The index of the first row of the first block T11.

*

*  N1      (input) INTEGER

*          The order of the first block T11. N1 = 0, 1 or 2.

*

*  N2      (input) INTEGER

*          The order of the second block T22. N2 = 0, 1 or 2.

*

*  ITRAF   (output) INTEGER array, length k, where

*             k = 1, if N1+N2 = 2;

*             k = 2, if N1+N2 = 3;

*             k = 4, if N1+N2 = 4.

*          List of parameters for representing the transformation

*          matrix Q, see BSTREXC.

*

*  DTRAF   (output) REAL             array, length k, where

*             k =  2, if N1+N2 = 2;

*             k =  5, if N1+N2 = 3;

*             k = 10, if N1+N2 = 4.

*          List of parameters for representing the transformation

*          matrix Q, see BSTREXC.

*

*  WORK    (workspace) REAL             array, dimension (N)

*

*  INFO    (output) INTEGER

*          = 0: successful exit

*          = 1: the transformed matrix T would be too far from Schur

*               form; the blocks are not swapped and T and Q are

*               unchanged.

*

*  =====================================================================

*

*     .. Parameters ..

      REAL               ZERO, ONE

      parameter( zero = 0.0e+0, one = 1.0e+0 )

      REAL               TEN

      parameter( ten = 10.0 )

      INTEGER            LDD, LDX

      parameter( ldd = 4, ldx = 2 )

*     ..

*     .. Local Scalars ..

      INTEGER            IERR, J2, J3, J4, K, LD, LI, ND

      REAL               CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22,

     $                   t33, tau, tau1, tau2, temp, thresh, wi1, wi2,

     $                   wr1, wr2, xnorm

*     ..

*     .. Local Arrays ..

      REAL               D( LDD, 4 ), X( LDX, 2 )

*     ..

*     .. External Functions ..

      REAL               SLAMCH, SLANGE

      EXTERNAL           slamch, slange

*     ..

*     .. External Subroutines ..

      EXTERNAL           slamov, slanv2, slarfg, slarfx, slartg, slasy2,

     $                   srot

*     ..

*     .. Intrinsic Functions ..

      INTRINSIC          abs, max

*     ..

*     .. Executable Statements ..

*

      info = 0

*

*     Quick return if possible

*

      IF( n.EQ.0 .OR. n1.EQ.0 .OR. n2.EQ.0 )

     $   RETURN

      IF( j1+n1.GT.n )

     $   RETURN

*

      j2 = j1 + 1

      j3 = j1 + 2

      j4 = j1 + 3

*

      IF( n1.EQ.1 .AND. n2.EQ.1 ) THEN

*

*        Swap two 1-by-1 blocks.

*

         t11 = t( j1, j1 )

         t22 = t( j2, j2 )

*

*        Determine the transformation to perform the interchange.

*

         CALL slartg( t( j1, j2 ), t22-t11, cs, sn, temp )

*

*        Apply transformation to the matrix T.

*

         IF( j3.LE.n )

     $      CALL srot( n-j1-1, t( j1, j3 ), ldt, t( j2, j3 ), ldt, cs,

     $                 sn )

         CALL srot( j1-1, t( 1, j1 ), 1, t( 1, j2 ), 1, cs, sn )

*

         t( j1, j1 ) = t22

         t( j2, j2 ) = t11

*

         itraf( 1 ) = j1

         dtraf( 1 ) = cs

         dtraf( 2 ) = sn

*

      ELSE

*

*        Swapping involves at least one 2-by-2 block.

*

*        Copy the diagonal block of order N1+N2 to the local array D

*        and compute its norm.

*

         nd = n1 + n2

         CALL slamov( 'Full', nd, nd, t( j1, j1 ), ldt, d, ldd )

         dnorm = slange( 'Max', nd, nd, d, ldd, work )

*

*        Compute machine-dependent threshold for test for accepting

*        swap.

*

         eps = slamch( 'P' )

         smlnum = slamch( 'S' ) / eps

         thresh = max( ten*eps*dnorm, smlnum )

*

*        Solve T11*X - X*T22 = scale*T12 for X.

*

         CALL slasy2( .false., .false., -1, n1, n2, d, ldd,

     $                d( n1+1, n1+1 ), ldd, d( 1, n1+1 ), ldd, scale, x,

     $                ldx, xnorm, ierr )

*

*        Swap the adjacent diagonal blocks.

*

         k = n1 + n1 + n2 - 3

         GO TO ( 10, 20, 30 )k

*

   10    CONTINUE

*

*        N1 = 1, N2 = 2: generate elementary reflector H so that:

*

*        ( scale, X11, X12 ) H = ( 0, 0, * )

*

         dtraf( 1 ) = scale

         dtraf( 2 ) = x( 1, 1 )

         dtraf( 3 ) = x( 1, 2 )

         CALL slarfg( 3, dtraf( 3 ), dtraf, 1, tau )

         dtraf( 3 ) = one

         t11 = t( j1, j1 )

*

*        Perform swap provisionally on diagonal block in D.

*

         CALL slarfx( 'Left', 3, 3, dtraf, tau, d, ldd, work )

         CALL slarfx( 'Right', 3, 3, dtraf, tau, d, ldd, work )

*

*        Test whether to reject swap.

*

         IF( max( abs( d( 3, 1 ) ), abs( d( 3, 2 ) ), abs( d( 3,

     $       3 )-t11 ) ).GT.thresh )GO TO 50

*

*        Accept swap: apply transformation to the entire matrix T.

*

         CALL slarfx( 'Left', 3, n-j1+1, dtraf, tau, t( j1, j1 ), ldt,

     $                work )

         CALL slarfx( 'Right', j2, 3, dtraf, tau, t( 1, j1 ), ldt,

     $                work )

*

         t( j3, j1 ) = zero

         t( j3, j2 ) = zero

         t( j3, j3 ) = t11

*

         itraf( 1 ) = 2*n + j1

         li = 2

         dtraf( 3 ) = tau

         ld = 4

         GO TO 40

*

   20    CONTINUE

*

*        N1 = 2, N2 = 1: generate elementary reflector H so that:

*

*        H (  -X11 ) = ( * )

*          (  -X21 ) = ( 0 )

*          ( scale ) = ( 0 )

*

         dtraf( 1 ) = -x( 1, 1 )

         dtraf( 2 ) = -x( 2, 1 )

         dtraf( 3 ) = scale

         CALL slarfg( 3, dtraf( 1 ), dtraf( 2 ), 1, tau )

         dtraf( 1 ) = one

         t33 = t( j3, j3 )

*

*        Perform swap provisionally on diagonal block in D.

*

         CALL slarfx( 'Left', 3, 3, dtraf, tau, d, ldd, work )

         CALL slarfx( 'Right', 3, 3, dtraf, tau, d, ldd, work )

*

*        Test whether to reject swap.

*

         IF( max( abs( d( 2, 1 ) ), abs( d( 3, 1 ) ), abs( d( 1,

     $       1 )-t33 ) ).GT.thresh )GO TO 50

*

*        Accept swap: apply transformation to the entire matrix T.

*

         CALL slarfx( 'Right', j3, 3, dtraf, tau, t( 1, j1 ), ldt,

     $                work )

         CALL slarfx( 'Left', 3, n-j1, dtraf, tau, t( j1, j2 ), ldt,

     $                work )

*

         t( j1, j1 ) = t33

         t( j2, j1 ) = zero

         t( j3, j1 ) = zero

*

         itraf( 1 ) = n + j1

         li = 2

         dtraf( 1 ) = tau

         ld = 4

         GO TO 40

*

   30    CONTINUE

*

*        N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so

*        that:

*

*        H(2) H(1) (  -X11  -X12 ) = (  *  * )

*                  (  -X21  -X22 )   (  0  * )

*                  ( scale    0  )   (  0  0 )

*                  (    0  scale )   (  0  0 )

*

         dtraf( 1 ) = -x( 1, 1 )

         dtraf( 2 ) = -x( 2, 1 )

         dtraf( 3 ) = scale

         CALL slarfg( 3, dtraf( 1 ), dtraf( 2 ), 1, tau1 )

         dtraf( 1 ) = one

*

         temp = -tau1*( x( 1, 2 )+dtraf( 2 )*x( 2, 2 ) )

         dtraf( 4 ) = -temp*dtraf( 2 ) - x( 2, 2 )

         dtraf( 5 ) = -temp*dtraf( 3 )

         dtraf( 6 ) = scale

         CALL slarfg( 3, dtraf( 4 ), dtraf( 5 ), 1, tau2 )

         dtraf( 4 ) = one

*

*        Perform swap provisionally on diagonal block in D.

*

         CALL slarfx( 'Left', 3, 4, dtraf, tau1, d, ldd, work )

         CALL slarfx( 'Right', 4, 3, dtraf, tau1, d, ldd, work )

         CALL slarfx( 'Left', 3, 4, dtraf( 4 ), tau2, d( 2, 1 ), ldd,

     $                work )

         CALL slarfx( 'Right', 4, 3, dtraf( 4 ), tau2, d( 1, 2 ), ldd,

     $                work )

*

*        Test whether to reject swap.

*

         IF( max( abs( d( 3, 1 ) ), abs( d( 3, 2 ) ), abs( d( 4, 1 ) ),

     $       abs( d( 4, 2 ) ) ).GT.thresh )GO TO 50

*

*        Accept swap: apply transformation to the entire matrix T.

*

         CALL slarfx( 'Left', 3, n-j1+1, dtraf, tau1, t( j1, j1 ), ldt,

     $                work )

         CALL slarfx( 'Right', j4, 3, dtraf, tau1, t( 1, j1 ), ldt,

     $                work )

         CALL slarfx( 'Left', 3, n-j1+1, dtraf( 4 ), tau2, t( j2, j1 ),

     $                ldt, work )

         CALL slarfx( 'Right', j4, 3, dtraf( 4 ), tau2, t( 1, j2 ), ldt,

     $                work )

*

         t( j3, j1 ) = zero

         t( j3, j2 ) = zero

         t( j4, j1 ) = zero

         t( j4, j2 ) = zero

*

         itraf( 1 ) = n + j1

         itraf( 2 ) = n + j2

         li = 3

         dtraf( 1 ) = tau1

         dtraf( 4 ) = tau2

         ld = 7

         GO TO 40

*

   40    CONTINUE

*

         IF( n2.EQ.2 ) THEN

*

*           Standardize new 2-by-2 block T11

*

            CALL slanv2( t( j1, j1 ), t( j1, j2 ), t( j2, j1 ),

     $                   t( j2, j2 ), wr1, wi1, wr2, wi2, cs, sn )

            CALL srot( n-j1-1, t( j1, j1+2 ), ldt, t( j2, j1+2 ), ldt,

     $                 cs, sn )

            CALL srot( j1-1, t( 1, j1 ), 1, t( 1, j2 ), 1, cs, sn )

            itraf( li ) = j1

            li = li + 1

            dtraf( ld ) = cs

            dtraf( ld+1 ) = sn

            ld = ld + 2

         END IF

*

         IF( n1.EQ.2 ) THEN

*

*           Standardize new 2-by-2 block T22

*

            j3 = j1 + n2

            j4 = j3 + 1

            CALL slanv2( t( j3, j3 ), t( j3, j4 ), t( j4, j3 ),

     $                   t( j4, j4 ), wr1, wi1, wr2, wi2, cs, sn )

            IF( j3+2.LE.n )

     $         CALL srot( n-j3-1, t( j3, j3+2 ), ldt, t( j4, j3+2 ),

     $                    ldt, cs, sn )

            CALL srot( j3-1, t( 1, j3 ), 1, t( 1, j4 ), 1, cs, sn )

            itraf( li ) = j3

            dtraf( ld ) = cs

            dtraf( ld+1 ) = sn

         END IF

*

      END IF

      RETURN

*

*     Exit with INFO = 1 if swap was rejected.

*

   50 CONTINUE

      info = 1

      RETURN

*

*     End of BSLAEXC

*


      END

bslaexc
subroutine bslaexc(n, t, ldt, j1, n1, n2, itraf, dtraf, work, info)
Definition bslaexc.f:3

max
#define max(A, B)
Definition pcgemr.c:180