← Back to branch summary

~ubuntu-branches/debian/sid/lammps/sid

« back to all changes in this revision

Viewing changes to lib/linalg/dlasq6.f

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2015-04-29 23:44:49 UTC
  • mfrom: (5.1.3 experimental)
  • Revision ID: package-import@ubuntu.com-20150429234449-mbhy9utku6hp6oq8
Tags: 0~20150313.gitfa668e1-1
Upload into unstable.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/linalg/dlasq6.f
 
1
*> \brief \b DLASQ6 computes one dqd transform in ping-pong form. Used by sbdsqr and sstegr.
 
2
*
 
3
*  =========== DOCUMENTATION ===========
 
4
*
 
5
* Online html documentation available at 
 
6
*            http://www.netlib.org/lapack/explore-html/ 
 
7
*
 
8
*> \htmlonly
 
9
*> Download DLASQ6 + dependencies 
 
10
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq6.f"> 
 
11
*> [TGZ]</a> 
 
12
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq6.f"> 
 
13
*> [ZIP]</a> 
 
14
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq6.f"> 
 
15
*> [TXT]</a>
 
16
*> \endhtmlonly 
 
17
*
 
18
*  Definition:
 
19
*  ===========
 
20
*
 
21
*       SUBROUTINE DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN,
 
22
*                          DNM1, DNM2 )
 
23
* 
 
24
*       .. Scalar Arguments ..
 
25
*       INTEGER            I0, N0, PP
 
26
*       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2
 
27
*       ..
 
28
*       .. Array Arguments ..
 
29
*       DOUBLE PRECISION   Z( * )
 
30
*       ..
 
31
*  
 
32
*
 
33
*> \par Purpose:
 
34
*  =============
 
35
*>
 
36
*> \verbatim
 
37
*>
 
38
*> DLASQ6 computes one dqd (shift equal to zero) transform in
 
39
*> ping-pong form, with protection against underflow and overflow.
 
40
*> \endverbatim
 
41
*
 
42
*  Arguments:
 
43
*  ==========
 
44
*
 
45
*> \param[in] I0
 
46
*> \verbatim
 
47
*>          I0 is INTEGER
 
48
*>        First index.
 
49
*> \endverbatim
 
50
*>
 
51
*> \param[in] N0
 
52
*> \verbatim
 
53
*>          N0 is INTEGER
 
54
*>        Last index.
 
55
*> \endverbatim
 
56
*>
 
57
*> \param[in] Z
 
58
*> \verbatim
 
59
*>          Z is DOUBLE PRECISION array, dimension ( 4*N )
 
60
*>        Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
 
61
*>        an extra argument.
 
62
*> \endverbatim
 
63
*>
 
64
*> \param[in] PP
 
65
*> \verbatim
 
66
*>          PP is INTEGER
 
67
*>        PP=0 for ping, PP=1 for pong.
 
68
*> \endverbatim
 
69
*>
 
70
*> \param[out] DMIN
 
71
*> \verbatim
 
72
*>          DMIN is DOUBLE PRECISION
 
73
*>        Minimum value of d.
 
74
*> \endverbatim
 
75
*>
 
76
*> \param[out] DMIN1
 
77
*> \verbatim
 
78
*>          DMIN1 is DOUBLE PRECISION
 
79
*>        Minimum value of d, excluding D( N0 ).
 
80
*> \endverbatim
 
81
*>
 
82
*> \param[out] DMIN2
 
83
*> \verbatim
 
84
*>          DMIN2 is DOUBLE PRECISION
 
85
*>        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
 
86
*> \endverbatim
 
87
*>
 
88
*> \param[out] DN
 
89
*> \verbatim
 
90
*>          DN is DOUBLE PRECISION
 
91
*>        d(N0), the last value of d.
 
92
*> \endverbatim
 
93
*>
 
94
*> \param[out] DNM1
 
95
*> \verbatim
 
96
*>          DNM1 is DOUBLE PRECISION
 
97
*>        d(N0-1).
 
98
*> \endverbatim
 
99
*>
 
100
*> \param[out] DNM2
 
101
*> \verbatim
 
102
*>          DNM2 is DOUBLE PRECISION
 
103
*>        d(N0-2).
 
104
*> \endverbatim
 
105
*
 
106
*  Authors:
 
107
*  ========
 
108
*
 
109
*> \author Univ. of Tennessee 
 
110
*> \author Univ. of California Berkeley 
 
111
*> \author Univ. of Colorado Denver 
 
112
*> \author NAG Ltd. 
 
113
*
 
114
*> \date September 2012
 
115
*
 
116
*> \ingroup auxOTHERcomputational
 
117
*
 
118
*  =====================================================================
 
119
      SUBROUTINE DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN,
 
120
     $                   DNM1, DNM2 )
 
121
*
 
122
*  -- LAPACK computational routine (version 3.4.2) --
 
123
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 
124
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 
125
*     September 2012
 
126
*
 
127
*     .. Scalar Arguments ..
 
128
      INTEGER            I0, N0, PP
 
129
      DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2
 
130
*     ..
 
131
*     .. Array Arguments ..
 
132
      DOUBLE PRECISION   Z( * )
 
133
*     ..
 
134
*
 
135
*  =====================================================================
 
136
*
 
137
*     .. Parameter ..
 
138
      DOUBLE PRECISION   ZERO
 
139
      PARAMETER          ( ZERO = 0.0D0 )
 
140
*     ..
 
141
*     .. Local Scalars ..
 
142
      INTEGER            J4, J4P2
 
143
      DOUBLE PRECISION   D, EMIN, SAFMIN, TEMP
 
144
*     ..
 
145
*     .. External Function ..
 
146
      DOUBLE PRECISION   DLAMCH
 
147
      EXTERNAL           DLAMCH
 
148
*     ..
 
149
*     .. Intrinsic Functions ..
 
150
      INTRINSIC          MIN
 
151
*     ..
 
152
*     .. Executable Statements ..
 
153
*
 
154
      IF( ( N0-I0-1 ).LE.0 )
 
155
     $   RETURN
 
156
*
 
157
      SAFMIN = DLAMCH( 'Safe minimum' )
 
158
      J4 = 4*I0 + PP - 3
 
159
      EMIN = Z( J4+4 ) 
 
160
      D = Z( J4 )
 
161
      DMIN = D
 
162
*
 
163
      IF( PP.EQ.0 ) THEN
 
164
         DO 10 J4 = 4*I0, 4*( N0-3 ), 4
 
165
            Z( J4-2 ) = D + Z( J4-1 ) 
 
166
            IF( Z( J4-2 ).EQ.ZERO ) THEN
 
167
               Z( J4 ) = ZERO
 
168
               D = Z( J4+1 )
 
169
               DMIN = D
 
170
               EMIN = ZERO
 
171
            ELSE IF( SAFMIN*Z( J4+1 ).LT.Z( J4-2 ) .AND.
 
172
     $               SAFMIN*Z( J4-2 ).LT.Z( J4+1 ) ) THEN
 
173
               TEMP = Z( J4+1 ) / Z( J4-2 )
 
174
               Z( J4 ) = Z( J4-1 )*TEMP
 
175
               D = D*TEMP
 
176
            ELSE 
 
177
               Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
 
178
               D = Z( J4+1 )*( D / Z( J4-2 ) )
 
179
            END IF
 
180
            DMIN = MIN( DMIN, D )
 
181
            EMIN = MIN( EMIN, Z( J4 ) )
 
182
   10    CONTINUE
 
183
      ELSE
 
184
         DO 20 J4 = 4*I0, 4*( N0-3 ), 4
 
185
            Z( J4-3 ) = D + Z( J4 ) 
 
186
            IF( Z( J4-3 ).EQ.ZERO ) THEN
 
187
               Z( J4-1 ) = ZERO
 
188
               D = Z( J4+2 )
 
189
               DMIN = D
 
190
               EMIN = ZERO
 
191
            ELSE IF( SAFMIN*Z( J4+2 ).LT.Z( J4-3 ) .AND.
 
192
     $               SAFMIN*Z( J4-3 ).LT.Z( J4+2 ) ) THEN
 
193
               TEMP = Z( J4+2 ) / Z( J4-3 )
 
194
               Z( J4-1 ) = Z( J4 )*TEMP
 
195
               D = D*TEMP
 
196
            ELSE 
 
197
               Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
 
198
               D = Z( J4+2 )*( D / Z( J4-3 ) )
 
199
            END IF
 
200
            DMIN = MIN( DMIN, D )
 
201
            EMIN = MIN( EMIN, Z( J4-1 ) )
 
202
   20    CONTINUE
 
203
      END IF
 
204
*
 
205
*     Unroll last two steps. 
 
206
*
 
207
      DNM2 = D
 
208
      DMIN2 = DMIN
 
209
      J4 = 4*( N0-2 ) - PP
 
210
      J4P2 = J4 + 2*PP - 1
 
211
      Z( J4-2 ) = DNM2 + Z( J4P2 )
 
212
      IF( Z( J4-2 ).EQ.ZERO ) THEN
 
213
         Z( J4 ) = ZERO
 
214
         DNM1 = Z( J4P2+2 )
 
215
         DMIN = DNM1
 
216
         EMIN = ZERO
 
217
      ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND.
 
218
     $         SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN
 
219
         TEMP = Z( J4P2+2 ) / Z( J4-2 )
 
220
         Z( J4 ) = Z( J4P2 )*TEMP
 
221
         DNM1 = DNM2*TEMP
 
222
      ELSE
 
223
         Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 
224
         DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) )
 
225
      END IF
 
226
      DMIN = MIN( DMIN, DNM1 )
 
227
*
 
228
      DMIN1 = DMIN
 
229
      J4 = J4 + 4
 
230
      J4P2 = J4 + 2*PP - 1
 
231
      Z( J4-2 ) = DNM1 + Z( J4P2 )
 
232
      IF( Z( J4-2 ).EQ.ZERO ) THEN
 
233
         Z( J4 ) = ZERO
 
234
         DN = Z( J4P2+2 )
 
235
         DMIN = DN
 
236
         EMIN = ZERO
 
237
      ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND.
 
238
     $         SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN
 
239
         TEMP = Z( J4P2+2 ) / Z( J4-2 )
 
240
         Z( J4 ) = Z( J4P2 )*TEMP
 
241
         DN = DNM1*TEMP
 
242
      ELSE
 
243
         Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 
244
         DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) )
 
245
      END IF
 
246
      DMIN = MIN( DMIN, DN )
 
247
*
 
248
      Z( J4+2 ) = DN
 
249
      Z( 4*N0-PP ) = EMIN
 
250
      RETURN
 
251
*
 
252
*     End of DLASQ6
 
253
*
 
254
      END