Routine: DLAMSH()  File: SRC\dlamsh.f

 
 
# lines: 236
  # code: 236
  # comment: 0
  # blank:0
# Variables:39
# Callers:0
# Callings:0
# Words:115
# Keywords:48
 

 

..
     .. Array Arguments ..
     ..
  Purpose
  =======
  DLAMSH 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.
  DLAMSH 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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.
  LDS     (local input) INTEGER
          On entry, the leading dimension of S.  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) DOUBLE PRECISION
          On entry, machine precision
          Unchanged on exit.
  Implemented by:  G. Henry, May 1, 1997
  =====================================================================
     .. Parameters ..

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

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