Routine: CLAMSH()  File: SRC\clamsh.f

 
 
# lines: 260
  # code: 260
  # comment: 0
  # blank:0
# Variables:42
# Callers:0
# Callings:0
# Words:92
# Keywords:52
 

 

..
     .. Array Arguments ..
     ..
  Purpose
  =======
  CLAMSH sends multiple shifts through a small (single node) matrix to
     see how consecutive small subdiagonal elements are modified by
     subsequent shifts in an effort to maximize the number of bulges
     that can be sent through.
  CLAMSH should only be called when there are multiple shifts/bulges
     (NBULGE > 1) and the first shift is starting in the middle of an
     unreduced Hessenberg matrix because of two or more consecutive
     small subdiagonal elements.
  Arguments
  =========
  S       (local input/output) COMPLEX array, ( LDS,* )
          On entry, the matrix of shifts.  Only the 2x2 diagonal of S
          is referenced.  It is assumed that S has JBLK double shifts
             (size 2).
          On exit, the data is rearranged in the best order for
             applying.
  LDS     (local input) INTEGER
          On entry, the leading dimension of S.  Unchanged on exit.
              1 < NBULGE <= JBLK <= LDS/2
  NBULGE  (local input/output) INTEGER
          On entry, the number of bulges to send through H ( >1 ).
              NBULGE should be less than the maximum determined (JBLK).
              1 < NBULGE <= JBLK <= LDS/2
          On exit, the maximum number of bulges that can be sent
              through.
  JBLK    (local input) INTEGER
          On entry, the number of shifts determined for S.
          Unchanged on exit.
  H       (local input/output) COMPLEX array ( LDH,N )
          On entry, the local matrix to apply the shifts on.
              H should be aligned so that the starting row is 2.
          On exit, the data is destroyed.
  LDH     (local input) INTEGER
          On entry, the leading dimension of H.  Unchanged on exit.
  N       (local input) INTEGER
          On entry, the size of H.  If all the bulges are expected to
              go through, N should be at least 4*NBULGE+2.
              Otherwise, NBULGE may be reduced by this routine.
  ULP     (local input) REAL
          On entry, machine precision
          Unchanged on exit.
  Further Details
  ===============
  Implemented by:  M. Fahey, May 28, 1999
  =====================================================================
     .. Parameters ..

 
Display dynamic version Find AutoScroll Reload FontSize: - + Hide Comments Hide Blanks Frame FullScreen MailPrint

 
001        SUBROUTINE CLAMSH( S , LDS , NBULGE , JBLK , H , LDH , N , ULP )
002  
003  *     -- ScaLAPACK routine(version 1.7) --
004  *     Univ. of Tennessee , Univ. of California Berkeley , NAG Ltd. ,
005  *     Courant Institute , Argonne National Lab , and Rice University
006  *     May 28 , 1999
007  
008  *     .. Scalar Arguments ..
009        INTEGER JBLK , LDH , LDS , N , NBULGE
010        REAL ULP
011        REAL RONE , TEN
012        PARAMETER( RONE = 1.0E + 0 , TEN = 10.0E + 0 )
013        COMPLEX ZERO
014        PARAMETER( ZERO =( 0.0E + 0 , 0.0E + 0 ) )
015  *     ..
016  *     .. Local Scalars ..
017        INTEGER I , IBULGE , IVAL , J , K , M , NR
018        REAL DVAL , S1 , TST1
019        COMPLEX CDUM , H00 , H10 , H11 , H12 , H21 , H22 , H33 , H33S ,
020       $H43H34 , H44 , H44S , SUM , T1 , T2 , T3 , V1 , V2 , V3
021  *     ..
022  *     .. Local Arrays ..
023        COMPLEX V( 3 )
024  *     ..
025  *     .. External Subroutines ..
026        EXTERNAL CCOPY , CLARFG
027  *     ..
028  *     .. Intrinsic Functions ..
029        INTRINSIC ABS , REAL , CONJG , AIMAG , MAX , MIN
030  *     ..
031  *     .. Statement Functions ..
032        REAL CABS1
033  *     ..
034  *     .. Statement Function definitions ..
035        CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
036  *     ..
037  *     .. Executable Statements ..
038  
039        M = 2
040        DO 50 IBULGE = 1 , NBULGE
041            H44 = S( 2*JBLK - 2*IBULGE + 2 , 2*JBLK - 2*IBULGE + 2 )
042            H33 = S( 2*JBLK - 2*IBULGE + 1 , 2*JBLK - 2*IBULGE + 1 )
043            H43H34 = S( 2*JBLK - 2*IBULGE + 1 , 2*JBLK - 2*IBULGE + 2 )*
044       $    S( 2*JBLK - 2*IBULGE + 2 , 2*JBLK - 2*IBULGE + 1 )
045            H11 = H( M , M )
046            H22 = H( M + 1 , M + 1 )
047            H21 = H( M + 1 , M )
048            H12 = H( M , M + 1 )
049            H44S = H44 - H11
050            H33S = H33 - H11
051            V1 =( H33S*H44S - H43H34 ) / H21 + H12
052            V2 = H22 - H11 - H33S - H44S
053            V3 = H( M + 2 , M + 1 )
054            S1 = CABS1( V1 ) + CABS1( V2 ) + CABS1( V3 )
055            V1 = V1 / S1
056            V2 = V2 / S1
057            V3 = V3 / S1
058            V( 1 ) = V1
059            V( 2 ) = V2
060            V( 3 ) = V3
061            H00 = H( M - 1 , M - 1 )
062            H10 = H( M , M - 1 )
063            TST1 = CABS1( V1 )*( CABS1( H00 ) + CABS1( H11 ) + CABS1( H22 ) )
064            IF( CABS1( H10 )*( CABS1( V2 ) + CABS1( V3 ) ).GT.ULP*TST1 ) THEN
065  *             Find minimum
066                DVAL =( CABS1( H10 )*( CABS1( V2 ) + CABS1( V3 ) ) ) /
067       $( ULP*TST1 )
068                IVAL = IBULGE
069                DO 10 I = IBULGE + 1 , NBULGE
070                    H44 = S( 2*JBLK - 2*I + 2 , 2*JBLK - 2*I + 2 )
071                    H33 = S( 2*JBLK - 2*I + 1 , 2*JBLK - 2*I + 1 )
072                    H43H34 = S( 2*JBLK - 2*I + 1 , 2*JBLK - 2*I + 2 )*
073       $            S( 2*JBLK - 2*I + 2 , 2*JBLK - 2*I + 1 )
074                    H11 = H( M , M )
075                    H22 = H( M + 1 , M + 1 )
076                    H21 = H( M + 1 , M )
077                    H12 = H( M , M + 1 )
078                    H44S = H44 - H11
079                    H33S = H33 - H11
080                    V1 =( H33S*H44S - H43H34 ) / H21 + H12
081                    V2 = H22 - H11 - H33S - H44S
082                    V3 = H( M + 2 , M + 1 )
083                    S1 = CABS1( V1 ) + CABS1( V2 ) + CABS1( V3 )
084                    V1 = V1 / S1
085                    V2 = V2 / S1
086                    V3 = V3 / S1
087                    V( 1 ) = V1
088                    V( 2 ) = V2
089                    V( 3 ) = V3
090                    H00 = H( M - 1 , M - 1 )
091                    H10 = H( M , M - 1 )
092                    TST1 = CABS1( V1 )*( CABS1( H00 ) + CABS1( H11 ) +
093       $            CABS1( H22 ) )
094                    IF(( DVAL.GT.( CABS1( H10 )*( CABS1( V2 ) +
095       $            CABS1( V3 ) ) ) / ( ULP*TST1 ) ) .AND.
096       $( DVAL.GT.RONE ) ) THEN
097                    DVAL =( CABS1( H10 )*( CABS1( V2 ) + CABS1( V3 ) ) ) /
098       $( ULP*TST1 )
099                    IVAL = I
100                END IF
101     10         CONTINUE
102                IF(( DVAL.LT.TEN ) .AND.( IVAL.NE.IBULGE ) ) THEN
103                    H44 = S( 2*JBLK - 2*IVAL + 2 , 2*JBLK - 2*IVAL + 2 )
104                    H33 = S( 2*JBLK - 2*IVAL + 1 , 2*JBLK - 2*IVAL + 1 )
105                    H43H34 = S( 2*JBLK - 2*IVAL + 1 , 2*JBLK - 2*IVAL + 2 )
106                    H10 = S( 2*JBLK - 2*IVAL + 2 , 2*JBLK - 2*IVAL + 1 )
107                    S( 2*JBLK - 2*IVAL + 2 , 2*JBLK - 2*IVAL + 2 ) = S( 2*JBLK - 2*
108       $            IBULGE + 2 , 2*JBLK - 2*IBULGE + 2 )
109                    S( 2*JBLK - 2*IVAL + 1 , 2*JBLK - 2*IVAL + 1 ) = S( 2*JBLK - 2*
110       $            IBULGE + 1 , 2*JBLK - 2*IBULGE + 1 )
111                    S( 2*JBLK - 2*IVAL + 1 , 2*JBLK - 2*IVAL + 2 ) = S( 2*JBLK - 2*
112       $            IBULGE + 1 , 2*JBLK - 2*IBULGE + 2 )
113                    S( 2*JBLK - 2*IVAL + 2 , 2*JBLK - 2*IVAL + 1 ) = S( 2*JBLK - 2*
114       $            IBULGE + 2 , 2*JBLK - 2*IBULGE + 1 )
115                    S( 2*JBLK - 2*IBULGE + 2 , 2*JBLK - 2*IBULGE + 2 ) = H44
116                    S( 2*JBLK - 2*IBULGE + 1 , 2*JBLK - 2*IBULGE + 1 ) = H33
117                    S( 2*JBLK - 2*IBULGE + 1 , 2*JBLK - 2*IBULGE + 2 ) = H43H34
118                    S( 2*JBLK - 2*IBULGE + 2 , 2*JBLK - 2*IBULGE + 1 ) = H10
119                END IF
120                H44 = S( 2*JBLK - 2*IBULGE + 2 , 2*JBLK - 2*IBULGE + 2 )
121                H33 = S( 2*JBLK - 2*IBULGE + 1 , 2*JBLK - 2*IBULGE + 1 )
122                H43H34 = S( 2*JBLK - 2*IBULGE + 1 , 2*JBLK - 2*IBULGE + 2 )*
123       $        S( 2*JBLK - 2*IBULGE + 2 , 2*JBLK - 2*IBULGE + 1 )
124                H11 = H( M , M )
125                H22 = H( M + 1 , M + 1 )
126                H21 = H( M + 1 , M )
127                H12 = H( M , M + 1 )
128                H44S = H44 - H11
129                H33S = H33 - H11
130                V1 =( H33S*H44S - H43H34 ) / H21 + H12
131                V2 = H22 - H11 - H33S - H44S
132                V3 = H( M + 2 , M + 1 )
133                S1 = CABS1( V1 ) + CABS1( V2 ) + CABS1( V3 )
134                V1 = V1 / S1
135                V2 = V2 / S1
136                V3 = V3 / S1
137                V( 1 ) = V1
138                V( 2 ) = V2
139                V( 3 ) = V3
140                H00 = H( M - 1 , M - 1 )
141                H10 = H( M , M - 1 )
142                TST1 = CABS1( V1 )*( CABS1( H00 ) + CABS1( H11 ) +
143       $        CABS1( H22 ) )
144            END IF
145            IF( CABS1( H10 )*( CABS1( V2 ) + CABS1( V3 ) ).GT.TEN*ULP*TST1 )
146       $        THEN
147  *             IBULGE better not be 1 here or we have a bug !
148                NBULGE = MAX( IBULGE - 1 , 1 )
149                RETURN
150            END IF
151            DO 40 K = M , N - 1
152                NR = MIN( 3 , N - K + 1 )
153                IF( K.GT.M )
154       $            CALL CCOPY( NR , H( K , K - 1 ) , 1 , V , 1 )
155                    CALL CLARFG( NR , V( 1 ) , V( 2 ) , 1 , T1 )
156                    IF( K.GT.M ) THEN
157                        H( K , K - 1 ) = V( 1 )
158                        H( K + 1 , K - 1 ) = ZERO
159                        IF( K.LT.N - 1 )
160       $                    H( K + 2 , K - 1 ) = ZERO
161                        ELSE
162  *                         H(m , m - 1) must be updated ,
163  
164                            H( K , K - 1 ) = H( K , K - 1 ) - CONJG( T1 )*H( K , K - 1 )
165                        END IF
166                        V2 = V( 2 )
167                        T2 = T1*V2
168                        IF( NR.EQ.3 ) THEN
169                            V3 = V( 3 )
170                            T3 = T1*V3
171                            DO 20 J = K , N
172                                SUM = CONJG( T1 )*H( K , J ) +
173       $                        CONJG( T2 )*H( K + 1 , J ) +
174       $                        CONJG( T3 )*H( K + 2 , J )
175                                H( K , J ) = H( K , J ) - SUM
176                                H( K + 1 , J ) = H( K + 1 , J ) - SUM*V2
177                                H( K + 2 , J ) = H( K + 2 , J ) - SUM*V3
178     20                     CONTINUE
179                            DO 30 J = 1 , MIN( K + 3 , N )
180                                SUM = T1*H( J , K ) + T2*H( J , K + 1 ) + T3*H( J , K + 2 )
181                                H( J , K ) = H( J , K ) - SUM
182                                H( J , K + 1 ) = H( J , K + 1 ) - SUM*CONJG( V2 )
183                                H( J , K + 2 ) = H( J , K + 2 ) - SUM*CONJG( V3 )
184     30                     CONTINUE
185                        END IF
186     40     CONTINUE
187     50 CONTINUE
188  
189        RETURN
190  
191  *     End of CLAMSH
192  
193        END