← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/sandbox/arpack/ARPACK/SRC/dstqrb.f

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/sandbox/arpack/ARPACK/SRC/dstqrb.f
1
 
c-----------------------------------------------------------------------
2
 
c\BeginDoc
3
 
c
4
 
c\Name: dstqrb
5
 
c
6
 
c\Description:
7
 
c  Computes all eigenvalues and the last component of the eigenvectors
8
 
c  of a symmetric tridiagonal matrix using the implicit QL or QR method.
9
 
c
10
 
c  This is mostly a modification of the LAPACK routine dsteqr.
11
 
c  See Remarks.
12
 
c
13
 
c\Usage:
14
 
c  call dstqrb
15
 
c     ( N, D, E, Z, WORK, INFO )
16
 
c
17
 
c\Arguments
18
 
c  N       Integer.  (INPUT)
19
 
c          The number of rows and columns in the matrix.  N >= 0.
20
 
c
21
 
c  D       Double precision array, dimension (N).  (INPUT/OUTPUT)
22
 
c          On entry, D contains the diagonal elements of the
23
 
c          tridiagonal matrix.
24
 
c          On exit, D contains the eigenvalues, in ascending order.
25
 
c          If an error exit is made, the eigenvalues are correct
26
 
c          for indices 1,2,...,INFO-1, but they are unordered and
27
 
c          may not be the smallest eigenvalues of the matrix.
28
 
c
29
 
c  E       Double precision array, dimension (N-1).  (INPUT/OUTPUT)
30
 
c          On entry, E contains the subdiagonal elements of the
31
 
c          tridiagonal matrix in positions 1 through N-1.
32
 
c          On exit, E has been destroyed.
33
 
c
34
 
c  Z       Double precision array, dimension (N).  (OUTPUT)
35
 
c          On exit, Z contains the last row of the orthonormal 
36
 
c          eigenvector matrix of the symmetric tridiagonal matrix.  
37
 
c          If an error exit is made, Z contains the last row of the
38
 
c          eigenvector matrix associated with the stored eigenvalues.
39
 
c
40
 
c  WORK    Double precision array, dimension (max(1,2*N-2)).  (WORKSPACE)
41
 
c          Workspace used in accumulating the transformation for 
42
 
c          computing the last components of the eigenvectors.
43
 
c
44
 
c  INFO    Integer.  (OUTPUT)
45
 
c          = 0:  normal return.
46
 
c          < 0:  if INFO = -i, the i-th argument had an illegal value.
47
 
c          > 0:  if INFO = +i, the i-th eigenvalue has not converged
48
 
c                              after a total of  30*N  iterations.
49
 
c
50
 
c\Remarks
51
 
c  1. None.
52
 
c
53
 
c-----------------------------------------------------------------------
54
 
c
55
 
c\BeginLib
56
 
c
57
 
c\Local variables:
58
 
c     xxxxxx  real
59
 
c
60
 
c\Routines called:
61
 
c     daxpy   Level 1 BLAS that computes a vector triad.
62
 
c     dcopy   Level 1 BLAS that copies one vector to another.
63
 
c     dswap   Level 1 BLAS that swaps the contents of two vectors.
64
 
c     lsame   LAPACK character comparison routine.
65
 
c     dlae2   LAPACK routine that computes the eigenvalues of a 2-by-2 
66
 
c             symmetric matrix.
67
 
c     dlaev2  LAPACK routine that eigendecomposition of a 2-by-2 symmetric 
68
 
c             matrix.
69
 
c     dlamch  LAPACK routine that determines machine constants.
70
 
c     dlanst  LAPACK routine that computes the norm of a matrix.
71
 
c     dlapy2  LAPACK routine to compute sqrt(x**2+y**2) carefully.
72
 
c     dlartg  LAPACK Givens rotation construction routine.
73
 
c     dlascl  LAPACK routine for careful scaling of a matrix.
74
 
c     dlaset  LAPACK matrix initialization routine.
75
 
c     dlasr   LAPACK routine that applies an orthogonal transformation to 
76
 
c             a matrix.
77
 
c     dlasrt  LAPACK sorting routine.
78
 
c     dsteqr  LAPACK routine that computes eigenvalues and eigenvectors
79
 
c             of a symmetric tridiagonal matrix.
80
 
c     xerbla  LAPACK error handler routine.
81
 
c
82
 
c\Authors
83
 
c     Danny Sorensen               Phuong Vu
84
 
c     Richard Lehoucq              CRPC / Rice University
85
 
c     Dept. of Computational &     Houston, Texas
86
 
c     Applied Mathematics
87
 
c     Rice University           
88
 
c     Houston, Texas            
89
 
c
90
 
c\SCCS Information: @(#) 
91
 
c FILE: stqrb.F   SID: 2.5   DATE OF SID: 8/27/96   RELEASE: 2
92
 
c
93
 
c\Remarks
94
 
c     1. Starting with version 2.5, this routine is a modified version
95
 
c        of LAPACK version 2.0 subroutine SSTEQR. No lines are deleted,
96
 
c        only commeted out and new lines inserted.
97
 
c        All lines commented out have "c$$$" at the beginning.
98
 
c        Note that the LAPACK version 1.0 subroutine SSTEQR contained
99
 
c        bugs. 
100
 
c
101
 
c\EndLib
102
 
c
103
 
c-----------------------------------------------------------------------
104
 
c
105
 
      subroutine dstqrb ( n, d, e, z, work, info )
106
 
c
107
 
c     %------------------%
108
 
c     | Scalar Arguments |
109
 
c     %------------------%
110
 
c
111
 
      integer    info, n
112
 
c
113
 
c     %-----------------%
114
 
c     | Array Arguments |
115
 
c     %-----------------%
116
 
c
117
 
      Double precision
118
 
     &           d( n ), e( n-1 ), z( n ), work( 2*n-2 )
119
 
c
120
 
c     .. parameters ..
121
 
      Double precision               
122
 
     &                   zero, one, two, three
123
 
      parameter          ( zero = 0.0D+0, one = 1.0D+0, 
124
 
     &                     two = 2.0D+0, three = 3.0D+0 )
125
 
      integer            maxit
126
 
      parameter          ( maxit = 30 )
127
 
c     ..
128
 
c     .. local scalars ..
129
 
      integer            i, icompz, ii, iscale, j, jtot, k, l, l1, lend,
130
 
     &                   lendm1, lendp1, lendsv, lm1, lsv, m, mm, mm1,
131
 
     &                   nm1, nmaxit
132
 
      Double precision               
133
 
     &                   anorm, b, c, eps, eps2, f, g, p, r, rt1, rt2,
134
 
     &                   s, safmax, safmin, ssfmax, ssfmin, tst
135
 
c     ..
136
 
c     .. external functions ..
137
 
      logical            lsame
138
 
      Double precision
139
 
     &                   dlamch, dlanst, dlapy2
140
 
      external           lsame, dlamch, dlanst, dlapy2
141
 
c     ..
142
 
c     .. external subroutines ..
143
 
      external           dlae2, dlaev2, dlartg, dlascl, dlaset, dlasr,
144
 
     &                   dlasrt, dswap, xerbla
145
 
c     ..
146
 
c     .. intrinsic functions ..
147
 
      intrinsic          abs, max, sign, sqrt
148
 
c     ..
149
 
c     .. executable statements ..
150
 
c
151
 
c     test the input parameters.
152
 
c
153
 
      info = 0
154
 
c
155
 
c$$$      IF( LSAME( COMPZ, 'N' ) ) THEN
156
 
c$$$         ICOMPZ = 0
157
 
c$$$      ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
158
 
c$$$         ICOMPZ = 1
159
 
c$$$      ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
160
 
c$$$         ICOMPZ = 2
161
 
c$$$      ELSE
162
 
c$$$         ICOMPZ = -1
163
 
c$$$      END IF
164
 
c$$$      IF( ICOMPZ.LT.0 ) THEN
165
 
c$$$         INFO = -1
166
 
c$$$      ELSE IF( N.LT.0 ) THEN
167
 
c$$$         INFO = -2
168
 
c$$$      ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1,
169
 
c$$$     $         N ) ) ) THEN
170
 
c$$$         INFO = -6
171
 
c$$$      END IF
172
 
c$$$      IF( INFO.NE.0 ) THEN
173
 
c$$$         CALL XERBLA( 'SSTEQR', -INFO )
174
 
c$$$         RETURN
175
 
c$$$      END IF
176
 
c
177
 
c    *** New starting with version 2.5 ***
178
 
c
179
 
      icompz = 2
180
 
c    *************************************
181
 
c
182
 
c     quick return if possible
183
 
c
184
 
      if( n.eq.0 )
185
 
     $   return
186
 
c
187
 
      if( n.eq.1 ) then
188
 
         if( icompz.eq.2 )  z( 1 ) = one
189
 
         return
190
 
      end if
191
 
c
192
 
c     determine the unit roundoff and over/underflow thresholds.
193
 
c
194
 
      eps = dlamch( 'e' )
195
 
      eps2 = eps**2
196
 
      safmin = dlamch( 's' )
197
 
      safmax = one / safmin
198
 
      ssfmax = sqrt( safmax ) / three
199
 
      ssfmin = sqrt( safmin ) / eps2
200
 
c
201
 
c     compute the eigenvalues and eigenvectors of the tridiagonal
202
 
c     matrix.
203
 
c
204
 
c$$      if( icompz.eq.2 )
205
 
c$$$     $   call dlaset( 'full', n, n, zero, one, z, ldz )
206
 
c
207
 
c     *** New starting with version 2.5 ***
208
 
c
209
 
      if ( icompz .eq. 2 ) then
210
 
         do 5 j = 1, n-1
211
 
            z(j) = zero
212
 
  5      continue
213
 
         z( n ) = one
214
 
      end if
215
 
c     *************************************
216
 
c
217
 
      nmaxit = n*maxit
218
 
      jtot = 0
219
 
c
220
 
c     determine where the matrix splits and choose ql or qr iteration
221
 
c     for each block, according to whether top or bottom diagonal
222
 
c     element is smaller.
223
 
c
224
 
      l1 = 1
225
 
      nm1 = n - 1
226
 
c
227
 
   10 continue
228
 
      if( l1.gt.n )
229
 
     $   go to 160
230
 
      if( l1.gt.1 )
231
 
     $   e( l1-1 ) = zero
232
 
      if( l1.le.nm1 ) then
233
 
         do 20 m = l1, nm1
234
 
            tst = abs( e( m ) )
235
 
            if( tst.eq.zero )
236
 
     $         go to 30
237
 
            if( tst.le.( sqrt( abs( d( m ) ) )*sqrt( abs( d( m+
238
 
     $          1 ) ) ) )*eps ) then
239
 
               e( m ) = zero
240
 
               go to 30
241
 
            end if
242
 
   20    continue
243
 
      end if
244
 
      m = n
245
 
c
246
 
   30 continue
247
 
      l = l1
248
 
      lsv = l
249
 
      lend = m
250
 
      lendsv = lend
251
 
      l1 = m + 1
252
 
      if( lend.eq.l )
253
 
     $   go to 10
254
 
c
255
 
c     scale submatrix in rows and columns l to lend
256
 
c
257
 
      anorm = dlanst( 'i', lend-l+1, d( l ), e( l ) )
258
 
      iscale = 0
259
 
      if( anorm.eq.zero )
260
 
     $   go to 10
261
 
      if( anorm.gt.ssfmax ) then
262
 
         iscale = 1
263
 
         call dlascl( 'g', 0, 0, anorm, ssfmax, lend-l+1, 1, d( l ), n,
264
 
     $                info )
265
 
         call dlascl( 'g', 0, 0, anorm, ssfmax, lend-l, 1, e( l ), n,
266
 
     $                info )
267
 
      else if( anorm.lt.ssfmin ) then
268
 
         iscale = 2
269
 
         call dlascl( 'g', 0, 0, anorm, ssfmin, lend-l+1, 1, d( l ), n,
270
 
     $                info )
271
 
         call dlascl( 'g', 0, 0, anorm, ssfmin, lend-l, 1, e( l ), n,
272
 
     $                info )
273
 
      end if
274
 
c
275
 
c     choose between ql and qr iteration
276
 
c
277
 
      if( abs( d( lend ) ).lt.abs( d( l ) ) ) then
278
 
         lend = lsv
279
 
         l = lendsv
280
 
      end if
281
 
c
282
 
      if( lend.gt.l ) then
283
 
c
284
 
c        ql iteration
285
 
c
286
 
c        look for small subdiagonal element.
287
 
c
288
 
   40    continue
289
 
         if( l.ne.lend ) then
290
 
            lendm1 = lend - 1
291
 
            do 50 m = l, lendm1
292
 
               tst = abs( e( m ) )**2
293
 
               if( tst.le.( eps2*abs( d( m ) ) )*abs( d( m+1 ) )+
294
 
     $             safmin )go to 60
295
 
   50       continue
296
 
         end if
297
 
c
298
 
         m = lend
299
 
c
300
 
   60    continue
301
 
         if( m.lt.lend )
302
 
     $      e( m ) = zero
303
 
         p = d( l )
304
 
         if( m.eq.l )
305
 
     $      go to 80
306
 
c
307
 
c        if remaining matrix is 2-by-2, use dlae2 or dlaev2
308
 
c        to compute its eigensystem.
309
 
c
310
 
         if( m.eq.l+1 ) then
311
 
            if( icompz.gt.0 ) then
312
 
               call dlaev2( d( l ), e( l ), d( l+1 ), rt1, rt2, c, s )
313
 
               work( l ) = c
314
 
               work( n-1+l ) = s
315
 
c$$$               call dlasr( 'r', 'v', 'b', n, 2, work( l ),
316
 
c$$$     $                     work( n-1+l ), z( 1, l ), ldz )
317
 
c
318
 
c              *** New starting with version 2.5 ***
319
 
c
320
 
               tst      = z(l+1)
321
 
               z(l+1) = c*tst - s*z(l)
322
 
               z(l)   = s*tst + c*z(l)
323
 
c              *************************************
324
 
            else
325
 
               call dlae2( d( l ), e( l ), d( l+1 ), rt1, rt2 )
326
 
            end if
327
 
            d( l ) = rt1
328
 
            d( l+1 ) = rt2
329
 
            e( l ) = zero
330
 
            l = l + 2
331
 
            if( l.le.lend )
332
 
     $         go to 40
333
 
            go to 140
334
 
         end if
335
 
c
336
 
         if( jtot.eq.nmaxit )
337
 
     $      go to 140
338
 
         jtot = jtot + 1
339
 
c
340
 
c        form shift.
341
 
c
342
 
         g = ( d( l+1 )-p ) / ( two*e( l ) )
343
 
         r = dlapy2( g, one )
344
 
         g = d( m ) - p + ( e( l ) / ( g+sign( r, g ) ) )
345
 
c
346
 
         s = one
347
 
         c = one
348
 
         p = zero
349
 
c
350
 
c        inner loop
351
 
c
352
 
         mm1 = m - 1
353
 
         do 70 i = mm1, l, -1
354
 
            f = s*e( i )
355
 
            b = c*e( i )
356
 
            call dlartg( g, f, c, s, r )
357
 
            if( i.ne.m-1 )
358
 
     $         e( i+1 ) = r
359
 
            g = d( i+1 ) - p
360
 
            r = ( d( i )-g )*s + two*c*b
361
 
            p = s*r
362
 
            d( i+1 ) = g + p
363
 
            g = c*r - b
364
 
c
365
 
c           if eigenvectors are desired, then save rotations.
366
 
c
367
 
            if( icompz.gt.0 ) then
368
 
               work( i ) = c
369
 
               work( n-1+i ) = -s
370
 
            end if
371
 
c
372
 
   70    continue
373
 
c
374
 
c        if eigenvectors are desired, then apply saved rotations.
375
 
c
376
 
         if( icompz.gt.0 ) then
377
 
            mm = m - l + 1
378
 
c$$$            call dlasr( 'r', 'v', 'b', n, mm, work( l ), work( n-1+l ),
379
 
c$$$     $                  z( 1, l ), ldz )
380
 
c
381
 
c             *** New starting with version 2.5 ***
382
 
c
383
 
              call dlasr( 'r', 'v', 'b', 1, mm, work( l ), 
384
 
     &                    work( n-1+l ), z( l ), 1 )
385
 
c             *************************************                             
386
 
         end if
387
 
c
388
 
         d( l ) = d( l ) - p
389
 
         e( l ) = g
390
 
         go to 40
391
 
c
392
 
c        eigenvalue found.
393
 
c
394
 
   80    continue
395
 
         d( l ) = p
396
 
c
397
 
         l = l + 1
398
 
         if( l.le.lend )
399
 
     $      go to 40
400
 
         go to 140
401
 
c
402
 
      else
403
 
c
404
 
c        qr iteration
405
 
c
406
 
c        look for small superdiagonal element.
407
 
c
408
 
   90    continue
409
 
         if( l.ne.lend ) then
410
 
            lendp1 = lend + 1
411
 
            do 100 m = l, lendp1, -1
412
 
               tst = abs( e( m-1 ) )**2
413
 
               if( tst.le.( eps2*abs( d( m ) ) )*abs( d( m-1 ) )+
414
 
     $             safmin )go to 110
415
 
  100       continue
416
 
         end if
417
 
c
418
 
         m = lend
419
 
c
420
 
  110    continue
421
 
         if( m.gt.lend )
422
 
     $      e( m-1 ) = zero
423
 
         p = d( l )
424
 
         if( m.eq.l )
425
 
     $      go to 130
426
 
c
427
 
c        if remaining matrix is 2-by-2, use dlae2 or dlaev2
428
 
c        to compute its eigensystem.
429
 
c
430
 
         if( m.eq.l-1 ) then
431
 
            if( icompz.gt.0 ) then
432
 
               call dlaev2( d( l-1 ), e( l-1 ), d( l ), rt1, rt2, c, s )
433
 
c$$$               work( m ) = c
434
 
c$$$               work( n-1+m ) = s
435
 
c$$$               call dlasr( 'r', 'v', 'f', n, 2, work( m ),
436
 
c$$$     $                     work( n-1+m ), z( 1, l-1 ), ldz )
437
 
c
438
 
c               *** New starting with version 2.5 ***
439
 
c
440
 
                tst      = z(l)
441
 
                z(l)   = c*tst - s*z(l-1)
442
 
                z(l-1) = s*tst + c*z(l-1)
443
 
c               ************************************* 
444
 
            else
445
 
               call dlae2( d( l-1 ), e( l-1 ), d( l ), rt1, rt2 )
446
 
            end if
447
 
            d( l-1 ) = rt1
448
 
            d( l ) = rt2
449
 
            e( l-1 ) = zero
450
 
            l = l - 2
451
 
            if( l.ge.lend )
452
 
     $         go to 90
453
 
            go to 140
454
 
         end if
455
 
c
456
 
         if( jtot.eq.nmaxit )
457
 
     $      go to 140
458
 
         jtot = jtot + 1
459
 
c
460
 
c        form shift.
461
 
c
462
 
         g = ( d( l-1 )-p ) / ( two*e( l-1 ) )
463
 
         r = dlapy2( g, one )
464
 
         g = d( m ) - p + ( e( l-1 ) / ( g+sign( r, g ) ) )
465
 
c
466
 
         s = one
467
 
         c = one
468
 
         p = zero
469
 
c
470
 
c        inner loop
471
 
c
472
 
         lm1 = l - 1
473
 
         do 120 i = m, lm1
474
 
            f = s*e( i )
475
 
            b = c*e( i )
476
 
            call dlartg( g, f, c, s, r )
477
 
            if( i.ne.m )
478
 
     $         e( i-1 ) = r
479
 
            g = d( i ) - p
480
 
            r = ( d( i+1 )-g )*s + two*c*b
481
 
            p = s*r
482
 
            d( i ) = g + p
483
 
            g = c*r - b
484
 
c
485
 
c           if eigenvectors are desired, then save rotations.
486
 
c
487
 
            if( icompz.gt.0 ) then
488
 
               work( i ) = c
489
 
               work( n-1+i ) = s
490
 
            end if
491
 
c
492
 
  120    continue
493
 
c
494
 
c        if eigenvectors are desired, then apply saved rotations.
495
 
c
496
 
         if( icompz.gt.0 ) then
497
 
            mm = l - m + 1
498
 
c$$$            call dlasr( 'r', 'v', 'f', n, mm, work( m ), work( n-1+m ),
499
 
c$$$     $                  z( 1, m ), ldz )
500
 
c
501
 
c           *** New starting with version 2.5 ***
502
 
c
503
 
            call dlasr( 'r', 'v', 'f', 1, mm, work( m ), work( n-1+m ),
504
 
     &                  z( m ), 1 )
505
 
c           *************************************                             
506
 
         end if
507
 
c
508
 
         d( l ) = d( l ) - p
509
 
         e( lm1 ) = g
510
 
         go to 90
511
 
c
512
 
c        eigenvalue found.
513
 
c
514
 
  130    continue
515
 
         d( l ) = p
516
 
c
517
 
         l = l - 1
518
 
         if( l.ge.lend )
519
 
     $      go to 90
520
 
         go to 140
521
 
c
522
 
      end if
523
 
c
524
 
c     undo scaling if necessary
525
 
c
526
 
  140 continue
527
 
      if( iscale.eq.1 ) then
528
 
         call dlascl( 'g', 0, 0, ssfmax, anorm, lendsv-lsv+1, 1,
529
 
     $                d( lsv ), n, info )
530
 
         call dlascl( 'g', 0, 0, ssfmax, anorm, lendsv-lsv, 1, e( lsv ),
531
 
     $                n, info )
532
 
      else if( iscale.eq.2 ) then
533
 
         call dlascl( 'g', 0, 0, ssfmin, anorm, lendsv-lsv+1, 1,
534
 
     $                d( lsv ), n, info )
535
 
         call dlascl( 'g', 0, 0, ssfmin, anorm, lendsv-lsv, 1, e( lsv ),
536
 
     $                n, info )
537
 
      end if
538
 
c
539
 
c     check for no convergence to an eigenvalue after a total
540
 
c     of n*maxit iterations.
541
 
c
542
 
      if( jtot.lt.nmaxit )
543
 
     $   go to 10
544
 
      do 150 i = 1, n - 1
545
 
         if( e( i ).ne.zero )
546
 
     $      info = info + 1
547
 
  150 continue
548
 
      go to 190
549
 
c
550
 
c     order eigenvalues and eigenvectors.
551
 
c
552
 
  160 continue
553
 
      if( icompz.eq.0 ) then
554
 
c
555
 
c        use quick sort
556
 
c
557
 
         call dlasrt( 'i', n, d, info )
558
 
c
559
 
      else
560
 
c
561
 
c        use selection sort to minimize swaps of eigenvectors
562
 
c
563
 
         do 180 ii = 2, n
564
 
            i = ii - 1
565
 
            k = i
566
 
            p = d( i )
567
 
            do 170 j = ii, n
568
 
               if( d( j ).lt.p ) then
569
 
                  k = j
570
 
                  p = d( j )
571
 
               end if
572
 
  170       continue
573
 
            if( k.ne.i ) then
574
 
               d( k ) = d( i )
575
 
               d( i ) = p
576
 
c$$$               call dswap( n, z( 1, i ), 1, z( 1, k ), 1 )
577
 
c           *** New starting with version 2.5 ***
578
 
c
579
 
               p    = z(k)
580
 
               z(k) = z(i)
581
 
               z(i) = p
582
 
c           *************************************
583
 
            end if
584
 
  180    continue
585
 
      end if
586
 
c
587
 
  190 continue
588
 
      return
589
 
c
590
 
c     %---------------%
591
 
c     | End of dstqrb |
592
 
c     %---------------%
593
 
c
594
 
      end