Statistiques
| Révision :

root / src / lapack / double / dlasq6.f @ 10

Historique | Voir | Annoter | Télécharger (5,09 ko)

1 1 pfleura2
      SUBROUTINE DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN,
2 1 pfleura2
     $                   DNM1, DNM2 )
3 1 pfleura2
*
4 1 pfleura2
*  -- LAPACK routine (version 3.2)                                    --
5 1 pfleura2
*
6 1 pfleura2
*  -- Contributed by Osni Marques of the Lawrence Berkeley National   --
7 1 pfleura2
*  -- Laboratory and Beresford Parlett of the Univ. of California at  --
8 1 pfleura2
*  -- Berkeley                                                        --
9 1 pfleura2
*  -- November 2008                                                   --
10 1 pfleura2
*
11 1 pfleura2
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
12 1 pfleura2
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
13 1 pfleura2
*
14 1 pfleura2
*     .. Scalar Arguments ..
15 1 pfleura2
      INTEGER            I0, N0, PP
16 1 pfleura2
      DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2
17 1 pfleura2
*     ..
18 1 pfleura2
*     .. Array Arguments ..
19 1 pfleura2
      DOUBLE PRECISION   Z( * )
20 1 pfleura2
*     ..
21 1 pfleura2
*
22 1 pfleura2
*  Purpose
23 1 pfleura2
*  =======
24 1 pfleura2
*
25 1 pfleura2
*  DLASQ6 computes one dqd (shift equal to zero) transform in
26 1 pfleura2
*  ping-pong form, with protection against underflow and overflow.
27 1 pfleura2
*
28 1 pfleura2
*  Arguments
29 1 pfleura2
*  =========
30 1 pfleura2
*
31 1 pfleura2
*  I0    (input) INTEGER
32 1 pfleura2
*        First index.
33 1 pfleura2
*
34 1 pfleura2
*  N0    (input) INTEGER
35 1 pfleura2
*        Last index.
36 1 pfleura2
*
37 1 pfleura2
*  Z     (input) DOUBLE PRECISION array, dimension ( 4*N )
38 1 pfleura2
*        Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
39 1 pfleura2
*        an extra argument.
40 1 pfleura2
*
41 1 pfleura2
*  PP    (input) INTEGER
42 1 pfleura2
*        PP=0 for ping, PP=1 for pong.
43 1 pfleura2
*
44 1 pfleura2
*  DMIN  (output) DOUBLE PRECISION
45 1 pfleura2
*        Minimum value of d.
46 1 pfleura2
*
47 1 pfleura2
*  DMIN1 (output) DOUBLE PRECISION
48 1 pfleura2
*        Minimum value of d, excluding D( N0 ).
49 1 pfleura2
*
50 1 pfleura2
*  DMIN2 (output) DOUBLE PRECISION
51 1 pfleura2
*        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
52 1 pfleura2
*
53 1 pfleura2
*  DN    (output) DOUBLE PRECISION
54 1 pfleura2
*        d(N0), the last value of d.
55 1 pfleura2
*
56 1 pfleura2
*  DNM1  (output) DOUBLE PRECISION
57 1 pfleura2
*        d(N0-1).
58 1 pfleura2
*
59 1 pfleura2
*  DNM2  (output) DOUBLE PRECISION
60 1 pfleura2
*        d(N0-2).
61 1 pfleura2
*
62 1 pfleura2
*  =====================================================================
63 1 pfleura2
*
64 1 pfleura2
*     .. Parameter ..
65 1 pfleura2
      DOUBLE PRECISION   ZERO
66 1 pfleura2
      PARAMETER          ( ZERO = 0.0D0 )
67 1 pfleura2
*     ..
68 1 pfleura2
*     .. Local Scalars ..
69 1 pfleura2
      INTEGER            J4, J4P2
70 1 pfleura2
      DOUBLE PRECISION   D, EMIN, SAFMIN, TEMP
71 1 pfleura2
*     ..
72 1 pfleura2
*     .. External Function ..
73 1 pfleura2
      DOUBLE PRECISION   DLAMCH
74 1 pfleura2
      EXTERNAL           DLAMCH
75 1 pfleura2
*     ..
76 1 pfleura2
*     .. Intrinsic Functions ..
77 1 pfleura2
      INTRINSIC          MIN
78 1 pfleura2
*     ..
79 1 pfleura2
*     .. Executable Statements ..
80 1 pfleura2
*
81 1 pfleura2
      IF( ( N0-I0-1 ).LE.0 )
82 1 pfleura2
     $   RETURN
83 1 pfleura2
*
84 1 pfleura2
      SAFMIN = DLAMCH( 'Safe minimum' )
85 1 pfleura2
      J4 = 4*I0 + PP - 3
86 1 pfleura2
      EMIN = Z( J4+4 )
87 1 pfleura2
      D = Z( J4 )
88 1 pfleura2
      DMIN = D
89 1 pfleura2
*
90 1 pfleura2
      IF( PP.EQ.0 ) THEN
91 1 pfleura2
         DO 10 J4 = 4*I0, 4*( N0-3 ), 4
92 1 pfleura2
            Z( J4-2 ) = D + Z( J4-1 )
93 1 pfleura2
            IF( Z( J4-2 ).EQ.ZERO ) THEN
94 1 pfleura2
               Z( J4 ) = ZERO
95 1 pfleura2
               D = Z( J4+1 )
96 1 pfleura2
               DMIN = D
97 1 pfleura2
               EMIN = ZERO
98 1 pfleura2
            ELSE IF( SAFMIN*Z( J4+1 ).LT.Z( J4-2 ) .AND.
99 1 pfleura2
     $               SAFMIN*Z( J4-2 ).LT.Z( J4+1 ) ) THEN
100 1 pfleura2
               TEMP = Z( J4+1 ) / Z( J4-2 )
101 1 pfleura2
               Z( J4 ) = Z( J4-1 )*TEMP
102 1 pfleura2
               D = D*TEMP
103 1 pfleura2
            ELSE
104 1 pfleura2
               Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
105 1 pfleura2
               D = Z( J4+1 )*( D / Z( J4-2 ) )
106 1 pfleura2
            END IF
107 1 pfleura2
            DMIN = MIN( DMIN, D )
108 1 pfleura2
            EMIN = MIN( EMIN, Z( J4 ) )
109 1 pfleura2
   10    CONTINUE
110 1 pfleura2
      ELSE
111 1 pfleura2
         DO 20 J4 = 4*I0, 4*( N0-3 ), 4
112 1 pfleura2
            Z( J4-3 ) = D + Z( J4 )
113 1 pfleura2
            IF( Z( J4-3 ).EQ.ZERO ) THEN
114 1 pfleura2
               Z( J4-1 ) = ZERO
115 1 pfleura2
               D = Z( J4+2 )
116 1 pfleura2
               DMIN = D
117 1 pfleura2
               EMIN = ZERO
118 1 pfleura2
            ELSE IF( SAFMIN*Z( J4+2 ).LT.Z( J4-3 ) .AND.
119 1 pfleura2
     $               SAFMIN*Z( J4-3 ).LT.Z( J4+2 ) ) THEN
120 1 pfleura2
               TEMP = Z( J4+2 ) / Z( J4-3 )
121 1 pfleura2
               Z( J4-1 ) = Z( J4 )*TEMP
122 1 pfleura2
               D = D*TEMP
123 1 pfleura2
            ELSE
124 1 pfleura2
               Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
125 1 pfleura2
               D = Z( J4+2 )*( D / Z( J4-3 ) )
126 1 pfleura2
            END IF
127 1 pfleura2
            DMIN = MIN( DMIN, D )
128 1 pfleura2
            EMIN = MIN( EMIN, Z( J4-1 ) )
129 1 pfleura2
   20    CONTINUE
130 1 pfleura2
      END IF
131 1 pfleura2
*
132 1 pfleura2
*     Unroll last two steps.
133 1 pfleura2
*
134 1 pfleura2
      DNM2 = D
135 1 pfleura2
      DMIN2 = DMIN
136 1 pfleura2
      J4 = 4*( N0-2 ) - PP
137 1 pfleura2
      J4P2 = J4 + 2*PP - 1
138 1 pfleura2
      Z( J4-2 ) = DNM2 + Z( J4P2 )
139 1 pfleura2
      IF( Z( J4-2 ).EQ.ZERO ) THEN
140 1 pfleura2
         Z( J4 ) = ZERO
141 1 pfleura2
         DNM1 = Z( J4P2+2 )
142 1 pfleura2
         DMIN = DNM1
143 1 pfleura2
         EMIN = ZERO
144 1 pfleura2
      ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND.
145 1 pfleura2
     $         SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN
146 1 pfleura2
         TEMP = Z( J4P2+2 ) / Z( J4-2 )
147 1 pfleura2
         Z( J4 ) = Z( J4P2 )*TEMP
148 1 pfleura2
         DNM1 = DNM2*TEMP
149 1 pfleura2
      ELSE
150 1 pfleura2
         Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
151 1 pfleura2
         DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) )
152 1 pfleura2
      END IF
153 1 pfleura2
      DMIN = MIN( DMIN, DNM1 )
154 1 pfleura2
*
155 1 pfleura2
      DMIN1 = DMIN
156 1 pfleura2
      J4 = J4 + 4
157 1 pfleura2
      J4P2 = J4 + 2*PP - 1
158 1 pfleura2
      Z( J4-2 ) = DNM1 + Z( J4P2 )
159 1 pfleura2
      IF( Z( J4-2 ).EQ.ZERO ) THEN
160 1 pfleura2
         Z( J4 ) = ZERO
161 1 pfleura2
         DN = Z( J4P2+2 )
162 1 pfleura2
         DMIN = DN
163 1 pfleura2
         EMIN = ZERO
164 1 pfleura2
      ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND.
165 1 pfleura2
     $         SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN
166 1 pfleura2
         TEMP = Z( J4P2+2 ) / Z( J4-2 )
167 1 pfleura2
         Z( J4 ) = Z( J4P2 )*TEMP
168 1 pfleura2
         DN = DNM1*TEMP
169 1 pfleura2
      ELSE
170 1 pfleura2
         Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
171 1 pfleura2
         DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) )
172 1 pfleura2
      END IF
173 1 pfleura2
      DMIN = MIN( DMIN, DN )
174 1 pfleura2
*
175 1 pfleura2
      Z( J4+2 ) = DN
176 1 pfleura2
      Z( 4*N0-PP ) = EMIN
177 1 pfleura2
      RETURN
178 1 pfleura2
*
179 1 pfleura2
*     End of DLASQ6
180 1 pfleura2
*
181 1 pfleura2
      END