← Back to branch summary

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

« back to all changes in this revision

Viewing changes to scipy/sandbox/arpack/ARPACK/SRC/sstqrb.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:
scipy/sandbox/arpack/ARPACK/SRC/sstqrb.f
 
1
c-----------------------------------------------------------------------
 
2
c\BeginDoc
 
3
c
 
4
c\Name: sstqrb
 
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 ssteqr.
 
11
c  See Remarks.
 
12
c
 
13
c\Usage:
 
14
c  call sstqrb
 
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       Real 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       Real 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       Real 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    Real 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     saxpy   Level 1 BLAS that computes a vector triad.
 
62
c     scopy   Level 1 BLAS that copies one vector to another.
 
63
c     sswap   Level 1 BLAS that swaps the contents of two vectors.
 
64
c     lsame   LAPACK character comparison routine.
 
65
c     slae2   LAPACK routine that computes the eigenvalues of a 2-by-2 
 
66
c             symmetric matrix.
 
67
c     slaev2  LAPACK routine that eigendecomposition of a 2-by-2 symmetric 
 
68
c             matrix.
 
69
c     slamch  LAPACK routine that determines machine constants.
 
70
c     slanst  LAPACK routine that computes the norm of a matrix.
 
71
c     slapy2  LAPACK routine to compute sqrt(x**2+y**2) carefully.
 
72
c     slartg  LAPACK Givens rotation construction routine.
 
73
c     slascl  LAPACK routine for careful scaling of a matrix.
 
74
c     slaset  LAPACK matrix initialization routine.
 
75
c     slasr   LAPACK routine that applies an orthogonal transformation to 
 
76
c             a matrix.
 
77
c     slasrt  LAPACK sorting routine.
 
78
c     ssteqr  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 sstqrb ( 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
      Real
 
118
     &           d( n ), e( n-1 ), z( n ), work( 2*n-2 )
 
119
c
 
120
c     .. parameters ..
 
121
      Real               
 
122
     &                   zero, one, two, three
 
123
      parameter          ( zero = 0.0E+0, one = 1.0E+0, 
 
124
     &                     two = 2.0E+0, three = 3.0E+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
      Real               
 
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
      Real
 
139
     &                   slamch, slanst, slapy2
 
140
      external           lsame, slamch, slanst, slapy2
 
141
c     ..
 
142
c     .. external subroutines ..
 
143
      external           slae2, slaev2, slartg, slascl, slaset, slasr,
 
144
     &                   slasrt, sswap, 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 = slamch( 'e' )
 
195
      eps2 = eps**2
 
196
      safmin = slamch( '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 slaset( '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 = slanst( '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 slascl( 'g', 0, 0, anorm, ssfmax, lend-l+1, 1, d( l ), n,
 
264
     $                info )
 
265
         call slascl( '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 slascl( 'g', 0, 0, anorm, ssfmin, lend-l+1, 1, d( l ), n,
 
270
     $                info )
 
271
         call slascl( '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 slae2 or slaev2
 
308
c        to compute its eigensystem.
 
309
c
 
310
         if( m.eq.l+1 ) then
 
311
            if( icompz.gt.0 ) then
 
312
               call slaev2( d( l ), e( l ), d( l+1 ), rt1, rt2, c, s )
 
313
               work( l ) = c
 
314
               work( n-1+l ) = s
 
315
c$$$               call slasr( '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 slae2( 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 = slapy2( 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 slartg( 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 slasr( '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 slasr( '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 slae2 or slaev2
 
428
c        to compute its eigensystem.
 
429
c
 
430
         if( m.eq.l-1 ) then
 
431
            if( icompz.gt.0 ) then
 
432
               call slaev2( 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 slasr( '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 slae2( 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 = slapy2( 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 slartg( 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 slasr( '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 slasr( '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 slascl( 'g', 0, 0, ssfmax, anorm, lendsv-lsv+1, 1,
 
529
     $                d( lsv ), n, info )
 
530
         call slascl( 'g', 0, 0, ssfmax, anorm, lendsv-lsv, 1, e( lsv ),
 
531
     $                n, info )
 
532
      else if( iscale.eq.2 ) then
 
533
         call slascl( 'g', 0, 0, ssfmin, anorm, lendsv-lsv+1, 1,
 
534
     $                d( lsv ), n, info )
 
535
         call slascl( '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 slasrt( '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 sswap( 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 sstqrb |
 
592
c     %---------------%
 
593
c
 
594
      end