← Back to branch summary

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

« back to all changes in this revision

Viewing changes to lib/linalg/dlasq5.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/dlasq5.f
 
1
*> \brief \b DLASQ5 computes one dqds 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 DLASQ5 + dependencies 
 
10
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq5.f"> 
 
11
*> [TGZ]</a> 
 
12
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq5.f"> 
 
13
*> [ZIP]</a> 
 
14
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq5.f"> 
 
15
*> [TXT]</a>
 
16
*> \endhtmlonly 
 
17
*
 
18
*  Definition:
 
19
*  ===========
 
20
*
 
21
*       SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN,
 
22
*                          DNM1, DNM2, IEEE, EPS )
 
23
* 
 
24
*       .. Scalar Arguments ..
 
25
*       LOGICAL            IEEE
 
26
*       INTEGER            I0, N0, PP
 
27
*       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU, SIGMA, EPS
 
28
*       ..
 
29
*       .. Array Arguments ..
 
30
*       DOUBLE PRECISION   Z( * )
 
31
*       ..
 
32
*  
 
33
*
 
34
*> \par Purpose:
 
35
*  =============
 
36
*>
 
37
*> \verbatim
 
38
*>
 
39
*> DLASQ5 computes one dqds transform in ping-pong form, one
 
40
*> version for IEEE machines another for non IEEE machines.
 
41
*> \endverbatim
 
42
*
 
43
*  Arguments:
 
44
*  ==========
 
45
*
 
46
*> \param[in] I0
 
47
*> \verbatim
 
48
*>          I0 is INTEGER
 
49
*>        First index.
 
50
*> \endverbatim
 
51
*>
 
52
*> \param[in] N0
 
53
*> \verbatim
 
54
*>          N0 is INTEGER
 
55
*>        Last index.
 
56
*> \endverbatim
 
57
*>
 
58
*> \param[in] Z
 
59
*> \verbatim
 
60
*>          Z is DOUBLE PRECISION array, dimension ( 4*N )
 
61
*>        Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
 
62
*>        an extra argument.
 
63
*> \endverbatim
 
64
*>
 
65
*> \param[in] PP
 
66
*> \verbatim
 
67
*>          PP is INTEGER
 
68
*>        PP=0 for ping, PP=1 for pong.
 
69
*> \endverbatim
 
70
*>
 
71
*> \param[in] TAU
 
72
*> \verbatim
 
73
*>          TAU is DOUBLE PRECISION
 
74
*>        This is the shift.
 
75
*> \endverbatim
 
76
*>
 
77
*> \param[in] SIGMA
 
78
*> \verbatim
 
79
*>          SIGMA is DOUBLE PRECISION
 
80
*>        This is the accumulated shift up to this step.
 
81
*> \endverbatim
 
82
*>
 
83
*> \param[out] DMIN
 
84
*> \verbatim
 
85
*>          DMIN is DOUBLE PRECISION
 
86
*>        Minimum value of d.
 
87
*> \endverbatim
 
88
*>
 
89
*> \param[out] DMIN1
 
90
*> \verbatim
 
91
*>          DMIN1 is DOUBLE PRECISION
 
92
*>        Minimum value of d, excluding D( N0 ).
 
93
*> \endverbatim
 
94
*>
 
95
*> \param[out] DMIN2
 
96
*> \verbatim
 
97
*>          DMIN2 is DOUBLE PRECISION
 
98
*>        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
 
99
*> \endverbatim
 
100
*>
 
101
*> \param[out] DN
 
102
*> \verbatim
 
103
*>          DN is DOUBLE PRECISION
 
104
*>        d(N0), the last value of d.
 
105
*> \endverbatim
 
106
*>
 
107
*> \param[out] DNM1
 
108
*> \verbatim
 
109
*>          DNM1 is DOUBLE PRECISION
 
110
*>        d(N0-1).
 
111
*> \endverbatim
 
112
*>
 
113
*> \param[out] DNM2
 
114
*> \verbatim
 
115
*>          DNM2 is DOUBLE PRECISION
 
116
*>        d(N0-2).
 
117
*> \endverbatim
 
118
*>
 
119
*> \param[in] IEEE
 
120
*> \verbatim
 
121
*>          IEEE is LOGICAL
 
122
*>        Flag for IEEE or non IEEE arithmetic.
 
123
*> \endverbatim
 
124
*
 
125
*> \param[in] EPS
 
126
*> \verbatim
 
127
*>          EPS is DOUBLE PRECISION
 
128
*>        This is the value of epsilon used.
 
129
*> \endverbatim
 
130
*>
 
131
*  Authors:
 
132
*  ========
 
133
*
 
134
*> \author Univ. of Tennessee 
 
135
*> \author Univ. of California Berkeley 
 
136
*> \author Univ. of Colorado Denver 
 
137
*> \author NAG Ltd. 
 
138
*
 
139
*> \date September 2012
 
140
*
 
141
*> \ingroup auxOTHERcomputational
 
142
*
 
143
*  =====================================================================
 
144
      SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2,
 
145
     $                   DN, DNM1, DNM2, IEEE, EPS )
 
146
*
 
147
*  -- LAPACK computational routine (version 3.4.2) --
 
148
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 
149
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 
150
*     September 2012
 
151
*
 
152
*     .. Scalar Arguments ..
 
153
      LOGICAL            IEEE
 
154
      INTEGER            I0, N0, PP
 
155
      DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU,
 
156
     $                   SIGMA, EPS
 
157
*     ..
 
158
*     .. Array Arguments ..
 
159
      DOUBLE PRECISION   Z( * )
 
160
*     ..
 
161
*
 
162
*  =====================================================================
 
163
*
 
164
*     .. Parameter ..
 
165
      DOUBLE PRECISION   ZERO, HALF
 
166
      PARAMETER          ( ZERO = 0.0D0, HALF = 0.5 )
 
167
*     ..
 
168
*     .. Local Scalars ..
 
169
      INTEGER            J4, J4P2
 
170
      DOUBLE PRECISION   D, EMIN, TEMP, DTHRESH
 
171
*     ..
 
172
*     .. Intrinsic Functions ..
 
173
      INTRINSIC          MIN
 
174
*     ..
 
175
*     .. Executable Statements ..
 
176
*
 
177
      IF( ( N0-I0-1 ).LE.0 )
 
178
     $   RETURN
 
179
*
 
180
      DTHRESH = EPS*(SIGMA+TAU)
 
181
      IF( TAU.LT.DTHRESH*HALF ) TAU = ZERO
 
182
      IF( TAU.NE.ZERO ) THEN
 
183
      J4 = 4*I0 + PP - 3
 
184
      EMIN = Z( J4+4 ) 
 
185
      D = Z( J4 ) - TAU
 
186
      DMIN = D
 
187
      DMIN1 = -Z( J4 )
 
188
*
 
189
      IF( IEEE ) THEN
 
190
*
 
191
*        Code for IEEE arithmetic.
 
192
*
 
193
         IF( PP.EQ.0 ) THEN
 
194
            DO 10 J4 = 4*I0, 4*( N0-3 ), 4
 
195
               Z( J4-2 ) = D + Z( J4-1 ) 
 
196
               TEMP = Z( J4+1 ) / Z( J4-2 )
 
197
               D = D*TEMP - TAU
 
198
               DMIN = MIN( DMIN, D )
 
199
               Z( J4 ) = Z( J4-1 )*TEMP
 
200
               EMIN = MIN( Z( J4 ), EMIN )
 
201
   10       CONTINUE
 
202
         ELSE
 
203
            DO 20 J4 = 4*I0, 4*( N0-3 ), 4
 
204
               Z( J4-3 ) = D + Z( J4 ) 
 
205
               TEMP = Z( J4+2 ) / Z( J4-3 )
 
206
               D = D*TEMP - TAU
 
207
               DMIN = MIN( DMIN, D )
 
208
               Z( J4-1 ) = Z( J4 )*TEMP
 
209
               EMIN = MIN( Z( J4-1 ), EMIN )
 
210
   20       CONTINUE
 
211
         END IF
 
212
*
 
213
*        Unroll last two steps. 
 
214
*
 
215
         DNM2 = D
 
216
         DMIN2 = DMIN
 
217
         J4 = 4*( N0-2 ) - PP
 
218
         J4P2 = J4 + 2*PP - 1
 
219
         Z( J4-2 ) = DNM2 + Z( J4P2 )
 
220
         Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 
221
         DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
 
222
         DMIN = MIN( DMIN, DNM1 )
 
223
*
 
224
         DMIN1 = DMIN
 
225
         J4 = J4 + 4
 
226
         J4P2 = J4 + 2*PP - 1
 
227
         Z( J4-2 ) = DNM1 + Z( J4P2 )
 
228
         Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 
229
         DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
 
230
         DMIN = MIN( DMIN, DN )
 
231
*
 
232
      ELSE
 
233
*
 
234
*        Code for non IEEE arithmetic.
 
235
*
 
236
         IF( PP.EQ.0 ) THEN
 
237
            DO 30 J4 = 4*I0, 4*( N0-3 ), 4
 
238
               Z( J4-2 ) = D + Z( J4-1 ) 
 
239
               IF( D.LT.ZERO ) THEN
 
240
                  RETURN
 
241
               ELSE 
 
242
                  Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
 
243
                  D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
 
244
               END IF
 
245
               DMIN = MIN( DMIN, D )
 
246
               EMIN = MIN( EMIN, Z( J4 ) )
 
247
   30       CONTINUE
 
248
         ELSE
 
249
            DO 40 J4 = 4*I0, 4*( N0-3 ), 4
 
250
               Z( J4-3 ) = D + Z( J4 ) 
 
251
               IF( D.LT.ZERO ) THEN
 
252
                  RETURN
 
253
               ELSE 
 
254
                  Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
 
255
                  D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
 
256
               END IF
 
257
               DMIN = MIN( DMIN, D )
 
258
               EMIN = MIN( EMIN, Z( J4-1 ) )
 
259
   40       CONTINUE
 
260
         END IF
 
261
*
 
262
*        Unroll last two steps. 
 
263
*
 
264
         DNM2 = D
 
265
         DMIN2 = DMIN
 
266
         J4 = 4*( N0-2 ) - PP
 
267
         J4P2 = J4 + 2*PP - 1
 
268
         Z( J4-2 ) = DNM2 + Z( J4P2 )
 
269
         IF( DNM2.LT.ZERO ) THEN
 
270
            RETURN
 
271
         ELSE
 
272
            Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 
273
            DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
 
274
         END IF
 
275
         DMIN = MIN( DMIN, DNM1 )
 
276
*
 
277
         DMIN1 = DMIN
 
278
         J4 = J4 + 4
 
279
         J4P2 = J4 + 2*PP - 1
 
280
         Z( J4-2 ) = DNM1 + Z( J4P2 )
 
281
         IF( DNM1.LT.ZERO ) THEN
 
282
            RETURN
 
283
         ELSE
 
284
            Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 
285
            DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
 
286
         END IF
 
287
         DMIN = MIN( DMIN, DN )
 
288
*
 
289
      END IF
 
290
      ELSE
 
291
*     This is the version that sets d's to zero if they are small enough
 
292
         J4 = 4*I0 + PP - 3
 
293
         EMIN = Z( J4+4 ) 
 
294
         D = Z( J4 ) - TAU
 
295
         DMIN = D
 
296
         DMIN1 = -Z( J4 )
 
297
         IF( IEEE ) THEN
 
298
*     
 
299
*     Code for IEEE arithmetic.
 
300
*     
 
301
            IF( PP.EQ.0 ) THEN
 
302
               DO 50 J4 = 4*I0, 4*( N0-3 ), 4
 
303
                  Z( J4-2 ) = D + Z( J4-1 ) 
 
304
                  TEMP = Z( J4+1 ) / Z( J4-2 )
 
305
                  D = D*TEMP - TAU
 
306
                  IF( D.LT.DTHRESH ) D = ZERO
 
307
                  DMIN = MIN( DMIN, D )
 
308
                  Z( J4 ) = Z( J4-1 )*TEMP
 
309
                  EMIN = MIN( Z( J4 ), EMIN )
 
310
 50            CONTINUE
 
311
            ELSE
 
312
               DO 60 J4 = 4*I0, 4*( N0-3 ), 4
 
313
                  Z( J4-3 ) = D + Z( J4 ) 
 
314
                  TEMP = Z( J4+2 ) / Z( J4-3 )
 
315
                  D = D*TEMP - TAU
 
316
                  IF( D.LT.DTHRESH ) D = ZERO
 
317
                  DMIN = MIN( DMIN, D )
 
318
                  Z( J4-1 ) = Z( J4 )*TEMP
 
319
                  EMIN = MIN( Z( J4-1 ), EMIN )
 
320
 60            CONTINUE
 
321
            END IF
 
322
*     
 
323
*     Unroll last two steps. 
 
324
*     
 
325
            DNM2 = D
 
326
            DMIN2 = DMIN
 
327
            J4 = 4*( N0-2 ) - PP
 
328
            J4P2 = J4 + 2*PP - 1
 
329
            Z( J4-2 ) = DNM2 + Z( J4P2 )
 
330
            Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 
331
            DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
 
332
            DMIN = MIN( DMIN, DNM1 )
 
333
*     
 
334
            DMIN1 = DMIN
 
335
            J4 = J4 + 4
 
336
            J4P2 = J4 + 2*PP - 1
 
337
            Z( J4-2 ) = DNM1 + Z( J4P2 )
 
338
            Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 
339
            DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
 
340
            DMIN = MIN( DMIN, DN )
 
341
*     
 
342
         ELSE
 
343
*     
 
344
*     Code for non IEEE arithmetic.
 
345
*     
 
346
            IF( PP.EQ.0 ) THEN
 
347
               DO 70 J4 = 4*I0, 4*( N0-3 ), 4
 
348
                  Z( J4-2 ) = D + Z( J4-1 ) 
 
349
                  IF( D.LT.ZERO ) THEN
 
350
                     RETURN
 
351
                  ELSE 
 
352
                     Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
 
353
                     D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
 
354
                  END IF
 
355
                  IF( D.LT.DTHRESH) D = ZERO
 
356
                  DMIN = MIN( DMIN, D )
 
357
                  EMIN = MIN( EMIN, Z( J4 ) )
 
358
 70            CONTINUE
 
359
            ELSE
 
360
               DO 80 J4 = 4*I0, 4*( N0-3 ), 4
 
361
                  Z( J4-3 ) = D + Z( J4 ) 
 
362
                  IF( D.LT.ZERO ) THEN
 
363
                     RETURN
 
364
                  ELSE 
 
365
                     Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
 
366
                     D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
 
367
                  END IF
 
368
                  IF( D.LT.DTHRESH) D = ZERO
 
369
                  DMIN = MIN( DMIN, D )
 
370
                  EMIN = MIN( EMIN, Z( J4-1 ) )
 
371
 80            CONTINUE
 
372
            END IF
 
373
*     
 
374
*     Unroll last two steps. 
 
375
*     
 
376
            DNM2 = D
 
377
            DMIN2 = DMIN
 
378
            J4 = 4*( N0-2 ) - PP
 
379
            J4P2 = J4 + 2*PP - 1
 
380
            Z( J4-2 ) = DNM2 + Z( J4P2 )
 
381
            IF( DNM2.LT.ZERO ) THEN
 
382
               RETURN
 
383
            ELSE
 
384
               Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 
385
               DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
 
386
            END IF
 
387
            DMIN = MIN( DMIN, DNM1 )
 
388
*     
 
389
            DMIN1 = DMIN
 
390
            J4 = J4 + 4
 
391
            J4P2 = J4 + 2*PP - 1
 
392
            Z( J4-2 ) = DNM1 + Z( J4P2 )
 
393
            IF( DNM1.LT.ZERO ) THEN
 
394
               RETURN
 
395
            ELSE
 
396
               Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
 
397
               DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
 
398
            END IF
 
399
            DMIN = MIN( DMIN, DN )
 
400
*     
 
401
         END IF
 
402
      END IF
 
403
*     
 
404
      Z( J4+2 ) = DN
 
405
      Z( 4*N0-PP ) = EMIN
 
406
      RETURN
 
407
*
 
408
*     End of DLASQ5
 
409
*
 
410
      END