Mercurial > octave-nkf

comparison libcruft/arpack/src/dstqrb.f @ 12274:9f5d2ef078e8 release-3-4-x

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
import ARPACK sources to libcruft from Debian package libarpack2 2.1+parpack96.dfsg-3+b1
author John W. Eaton <jwe@octave.org>
date Fri, 28 Jan 2011 14:04:33 -0500
parents
children
comparison
equal deleted inserted replaced
12273:83133b5bf392 12274:9f5d2ef078e8
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