« back to all changes in this revision
Viewing changes to SRC/slaqrb.f
- browse files at revision 27
- compare with another revision
- download diff
- download tarball
- view history from revision 27
- Committer: Package Import Robot
- Author(s): Sylvestre Ledru
- Date: 2015-10-12 09:33:05 UTC
- mfrom: (1.3.5)
- Revision ID: package-import@ubuntu.com-20151012093305-w41sv3g9nzbaodfz
Tags: 3.3.0-1
* New upstream release
* Update of the homepage (moved to github)
* libparpack links against libarpack (instead of static build)
* Update the symbol files. No need to change of SONAME:
- About the parpack change: https://github.com/opencollab/arpack-ng/pull/21
- Transfer the call to the library:
https://github.com/opencollab/arpack-ng/pull/19
* Standards-Version updated to version 3.9.6
* Update of the homepage (moved to github)
* libparpack links against libarpack (instead of static build)
* Update the symbol files. No need to change of SONAME:
- About the parpack change: https://github.com/opencollab/arpack-ng/pull/21
- Transfer the call to the library:
https://github.com/opencollab/arpack-ng/pull/19
* Standards-Version updated to version 3.9.6
- files added:
- .gitignore
- files removed:
- EXAMPLES/BAND/.deps
- files modified:
- CHANGES
added
removed
SRC/slaqrb.f
1
c-----------------------------------------------------------------------
2
c\BeginDoc
3
c
4
c\Name: slaqrb
5
c
6
c\Description:
7
c Compute the eigenvalues and the Schur decomposition of an upper
8
c Hessenberg submatrix in rows and columns ILO to IHI. Only the
9
c last component of the Schur vectors are computed.
10
c
11
c This is mostly a modification of the LAPACK routine slahqr.
12
c
13
c\Usage:
14
c call slaqrb
15
c ( WANTT, N, ILO, IHI, H, LDH, WR, WI, Z, INFO )
16
c
17
c\Arguments
18
c WANTT Logical variable. (INPUT)
19
c = .TRUE. : the full Schur form T is required;
20
c = .FALSE.: only eigenvalues are required.
21
c
22
c N Integer. (INPUT)
23
c The order of the matrix H. N >= 0.
24
c
25
c ILO Integer. (INPUT)
26
c IHI Integer. (INPUT)
27
c It is assumed that H is already upper quasi-triangular in
28
c rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless
29
c ILO = 1). SLAQRB works primarily with the Hessenberg
30
c submatrix in rows and columns ILO to IHI, but applies
31
c transformations to all of H if WANTT is .TRUE..
32
c 1 <= ILO <= max(1,IHI); IHI <= N.
33
c
34
c H Real array, dimension (LDH,N). (INPUT/OUTPUT)
35
c On entry, the upper Hessenberg matrix H.
36
c On exit, if WANTT is .TRUE., H is upper quasi-triangular in
37
c rows and columns ILO:IHI, with any 2-by-2 diagonal blocks in
38
c standard form. If WANTT is .FALSE., the contents of H are
39
c unspecified on exit.
40
c
41
c LDH Integer. (INPUT)
42
c The leading dimension of the array H. LDH >= max(1,N).
43
c
44
c WR Real array, dimension (N). (OUTPUT)
45
c WI Real array, dimension (N). (OUTPUT)
46
c The real and imaginary parts, respectively, of the computed
47
c eigenvalues ILO to IHI are stored in the corresponding
48
c elements of WR and WI. If two eigenvalues are computed as a
49
c complex conjugate pair, they are stored in consecutive
50
c elements of WR and WI, say the i-th and (i+1)th, with
51
c WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the
52
c eigenvalues are stored in the same order as on the diagonal
53
c of the Schur form returned in H, with WR(i) = H(i,i), and, if
54
c H(i:i+1,i:i+1) is a 2-by-2 diagonal block,
55
c WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i).
56
c
57
c Z Real array, dimension (N). (OUTPUT)
58
c On exit Z contains the last components of the Schur vectors.
59
c
60
c INFO Integer. (OUPUT)
61
c = 0: successful exit
62
c > 0: SLAQRB failed to compute all the eigenvalues ILO to IHI
63
c in a total of 30*(IHI-ILO+1) iterations; if INFO = i,
64
c elements i+1:ihi of WR and WI contain those eigenvalues
65
c which have been successfully computed.
66
c
67
c\Remarks
68
c 1. None.
69
c
70
c-----------------------------------------------------------------------
71
c
72
c\BeginLib
73
c
74
c\Local variables:
75
c xxxxxx real
76
c
77
c\Routines called:
78
c slabad LAPACK routine that computes machine constants.
79
c slamch LAPACK routine that determines machine constants.
80
c slanhs LAPACK routine that computes various norms of a matrix.
81
c slanv2 LAPACK routine that computes the Schur factorization of
82
c 2 by 2 nonsymmetric matrix in standard form.
83
c slarfg LAPACK Householder reflection construction routine.
84
c scopy Level 1 BLAS that copies one vector to another.
85
c srot Level 1 BLAS that applies a rotation to a 2 by 2 matrix.
86
87
c
88
c\Author
89
c Danny Sorensen Phuong Vu
90
c Richard Lehoucq CRPC / Rice University
91
c Dept. of Computational & Houston, Texas
92
c Applied Mathematics
93
c Rice University
94
c Houston, Texas
95
c
96
c\Revision history:
97
c xx/xx/92: Version ' 2.4'
98
c Modified from the LAPACK routine slahqr so that only the
99
c last component of the Schur vectors are computed.
100
c
101
c\SCCS Information: @(#)
102
c FILE: laqrb.F SID: 2.2 DATE OF SID: 8/27/96 RELEASE: 2
103
c
104
c\Remarks
105
c 1. None
106
c
107
c\EndLib
108
c
109
c-----------------------------------------------------------------------
110
c
111
subroutine slaqrb ( wantt, n, ilo, ihi, h, ldh, wr, wi,
112
& z, info )
113
c
114
c %------------------%
115
c | Scalar Arguments |
116
c %------------------%
117
c
118
logical wantt
119
integer ihi, ilo, info, ldh, n
120
c
121
c %-----------------%
122
c | Array Arguments |
123
c %-----------------%
124
c
125
Real
126
& h( ldh, * ), wi( * ), wr( * ), z( * )
127
c
128
c %------------%
129
c | Parameters |
130
c %------------%
131
c
132
Real
133
& zero, one, dat1, dat2
134
parameter (zero = 0.0E+0, one = 1.0E+0, dat1 = 7.5E-1,
135
& dat2 = -4.375E-1)
136
c
137
c %------------------------%
138
c | Local Scalars & Arrays |
139
c %------------------------%
140
c
141
integer i, i1, i2, itn, its, j, k, l, m, nh, nr
142
Real
143
& cs, h00, h10, h11, h12, h21, h22, h33, h33s,
144
& h43h34, h44, h44s, ovfl, s, smlnum, sn, sum,
145
& t1, t2, t3, tst1, ulp, unfl, v1, v2, v3
146
Real
147
& v( 3 ), work( 1 )
148
c
149
c %--------------------%
150
c | External Functions |
151
c %--------------------%
152
c
153
Real
154
& slamch, slanhs
155
external slamch, slanhs
156
c
157
c %----------------------%
158
c | External Subroutines |
159
c %----------------------%
160
c
161
external scopy, slabad, slanv2, slarfg, srot
162
c
163
c %-----------------------%
164
c | Executable Statements |
165
c %-----------------------%
166
c
167
info = 0
168
c
169
c %--------------------------%
170
c | Quick return if possible |
171
c %--------------------------%
172
c
173
if( n.eq.0 )
174
& return
175
if( ilo.eq.ihi ) then
176
wr( ilo ) = h( ilo, ilo )
177
wi( ilo ) = zero
178
return
179
end if
180
c
181
c %---------------------------------------------%
182
c | Initialize the vector of last components of |
183
c | the Schur vectors for accumulation. |
184
c %---------------------------------------------%
185
c
186
do 5 j = 1, n-1
187
z(j) = zero
188
5 continue
189
z(n) = one
190
c
191
nh = ihi - ilo + 1
192
c
193
c %-------------------------------------------------------------%
194
c | Set machine-dependent constants for the stopping criterion. |
195
c | If norm(H) <= sqrt(OVFL), overflow should not occur. |
196
c %-------------------------------------------------------------%
197
c
198
unfl = slamch( 'safe minimum' )
199
ovfl = one / unfl
200
call slabad( unfl, ovfl )
201
ulp = slamch( 'precision' )
202
smlnum = unfl*( nh / ulp )
203
c
204
c %---------------------------------------------------------------%
205
c | I1 and I2 are the indices of the first row and last column |
206
c | of H to which transformations must be applied. If eigenvalues |
207
c | only are computed, I1 and I2 are set inside the main loop. |
208
c | Zero out H(J+2,J) = ZERO for J=1:N if WANTT = .TRUE. |
209
c | else H(J+2,J) for J=ILO:IHI-ILO-1 if WANTT = .FALSE. |
210
c %---------------------------------------------------------------%
211
c
212
if( wantt ) then
213
i1 = 1
214
i2 = n
215
do 8 i=1,i2-2
216
h(i1+i+1,i) = zero
217
8 continue
218
else
219
do 9 i=1, ihi-ilo-1
220
h(ilo+i+1,ilo+i-1) = zero
221
9 continue
222
end if
223
c
224
c %---------------------------------------------------%
225
c | ITN is the total number of QR iterations allowed. |
226
c %---------------------------------------------------%
227
c
228
itn = 30*nh
229
c
230
c ------------------------------------------------------------------
231
c The main loop begins here. I is the loop index and decreases from
232
c IHI to ILO in steps of 1 or 2. Each iteration of the loop works
233
c with the active submatrix in rows and columns L to I.
234
c Eigenvalues I+1 to IHI have already converged. Either L = ILO or
235
c H(L,L-1) is negligible so that the matrix splits.
236
c ------------------------------------------------------------------
237
c
238
i = ihi
239
10 continue
240
l = ilo
241
if( i.lt.ilo )
242
& go to 150
243
244
c %--------------------------------------------------------------%
245
c | Perform QR iterations on rows and columns ILO to I until a |
246
c | submatrix of order 1 or 2 splits off at the bottom because a |
247
c | subdiagonal element has become negligible. |
248
c %--------------------------------------------------------------%
249
250
do 130 its = 0, itn
251
c
252
c %----------------------------------------------%
253
c | Look for a single small subdiagonal element. |
254
c %----------------------------------------------%
255
c
256
do 20 k = i, l + 1, -1
257
tst1 = abs( h( k-1, k-1 ) ) + abs( h( k, k ) )
258
if( tst1.eq.zero )
259
& tst1 = slanhs( '1', i-l+1, h( l, l ), ldh, work )
260
if( abs( h( k, k-1 ) ).le.max( ulp*tst1, smlnum ) )
261
& go to 30
262
20 continue
263
30 continue
264
l = k
265
if( l.gt.ilo ) then
266
c
267
c %------------------------%
268
c | H(L,L-1) is negligible |
269
c %------------------------%
270
c
271
h( l, l-1 ) = zero
272
end if
273
c
274
c %-------------------------------------------------------------%
275
c | Exit from loop if a submatrix of order 1 or 2 has split off |
276
c %-------------------------------------------------------------%
277
c
278
if( l.ge.i-1 )
279
& go to 140
280
c
281
c %---------------------------------------------------------%
282
c | Now the active submatrix is in rows and columns L to I. |
283
c | If eigenvalues only are being computed, only the active |
284
c | submatrix need be transformed. |
285
c %---------------------------------------------------------%
286
c
287
if( .not.wantt ) then
288
i1 = l
289
i2 = i
290
end if
291
c
292
if( its.eq.10 .or. its.eq.20 ) then
293
c
294
c %-------------------%
295
c | Exceptional shift |
296
c %-------------------%
297
c
298
s = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) )
299
h44 = dat1*s
300
h33 = h44
301
h43h34 = dat2*s*s
302
c
303
else
304
c
305
c %-----------------------------------------%
306
c | Prepare to use Wilkinson's double shift |
307
c %-----------------------------------------%
308
c
309
h44 = h( i, i )
310
h33 = h( i-1, i-1 )
311
h43h34 = h( i, i-1 )*h( i-1, i )
312
end if
313
c
314
c %-----------------------------------------------------%
315
c | Look for two consecutive small subdiagonal elements |
316
c %-----------------------------------------------------%
317
c
318
do 40 m = i - 2, l, -1
319
c
320
c %---------------------------------------------------------%
321
c | Determine the effect of starting the double-shift QR |
322
c | iteration at row M, and see if this would make H(M,M-1) |
323
c | negligible. |
324
c %---------------------------------------------------------%
325
c
326
h11 = h( m, m )
327
h22 = h( m+1, m+1 )
328
h21 = h( m+1, m )
329
h12 = h( m, m+1 )
330
h44s = h44 - h11
331
h33s = h33 - h11
332
v1 = ( h33s*h44s-h43h34 ) / h21 + h12
333
v2 = h22 - h11 - h33s - h44s
334
v3 = h( m+2, m+1 )
335
s = abs( v1 ) + abs( v2 ) + abs( v3 )
336
v1 = v1 / s
337
v2 = v2 / s
338
v3 = v3 / s
339
v( 1 ) = v1
340
v( 2 ) = v2
341
v( 3 ) = v3
342
if( m.eq.l )
343
& go to 50
344
h00 = h( m-1, m-1 )
345
h10 = h( m, m-1 )
346
tst1 = abs( v1 )*( abs( h00 )+abs( h11 )+abs( h22 ) )
347
if( abs( h10 )*( abs( v2 )+abs( v3 ) ).le.ulp*tst1 )
348
& go to 50
349
40 continue
350
50 continue
351
c
352
c %----------------------%
353
c | Double-shift QR step |
354
c %----------------------%
355
c
356
do 120 k = m, i - 1
357
c
358
c ------------------------------------------------------------
359
c The first iteration of this loop determines a reflection G
360
c from the vector V and applies it from left and right to H,
361
c thus creating a nonzero bulge below the subdiagonal.
362
c
363
c Each subsequent iteration determines a reflection G to
364
c restore the Hessenberg form in the (K-1)th column, and thus
365
c chases the bulge one step toward the bottom of the active
366
c submatrix. NR is the order of G.
367
c ------------------------------------------------------------
368
c
369
nr = min( 3, i-k+1 )
370
if( k.gt.m )
371
& call scopy( nr, h( k, k-1 ), 1, v, 1 )
372
call slarfg( nr, v( 1 ), v( 2 ), 1, t1 )
373
if( k.gt.m ) then
374
h( k, k-1 ) = v( 1 )
375
h( k+1, k-1 ) = zero
376
if( k.lt.i-1 )
377
& h( k+2, k-1 ) = zero
378
else if( m.gt.l ) then
379
h( k, k-1 ) = -h( k, k-1 )
380
end if
381
v2 = v( 2 )
382
t2 = t1*v2
383
if( nr.eq.3 ) then
384
v3 = v( 3 )
385
t3 = t1*v3
386
c
387
c %------------------------------------------------%
388
c | Apply G from the left to transform the rows of |
389
c | the matrix in columns K to I2. |
390
c %------------------------------------------------%
391
c
392
do 60 j = k, i2
393
sum = h( k, j ) + v2*h( k+1, j ) + v3*h( k+2, j )
394
h( k, j ) = h( k, j ) - sum*t1
395
h( k+1, j ) = h( k+1, j ) - sum*t2
396
h( k+2, j ) = h( k+2, j ) - sum*t3
397
60 continue
398
c
399
c %----------------------------------------------------%
400
c | Apply G from the right to transform the columns of |
401
c | the matrix in rows I1 to min(K+3,I). |
402
c %----------------------------------------------------%
403
c
404
do 70 j = i1, min( k+3, i )
405
sum = h( j, k ) + v2*h( j, k+1 ) + v3*h( j, k+2 )
406
h( j, k ) = h( j, k ) - sum*t1
407
h( j, k+1 ) = h( j, k+1 ) - sum*t2
408
h( j, k+2 ) = h( j, k+2 ) - sum*t3
409
70 continue
410
c
411
c %----------------------------------%
412
c | Accumulate transformations for Z |
413
c %----------------------------------%
414
c
415
sum = z( k ) + v2*z( k+1 ) + v3*z( k+2 )
416
z( k ) = z( k ) - sum*t1
417
z( k+1 ) = z( k+1 ) - sum*t2
418
z( k+2 ) = z( k+2 ) - sum*t3
419
420
else if( nr.eq.2 ) then
421
c
422
c %------------------------------------------------%
423
c | Apply G from the left to transform the rows of |
424
c | the matrix in columns K to I2. |
425
c %------------------------------------------------%
426
c
427
do 90 j = k, i2
428
sum = h( k, j ) + v2*h( k+1, j )
429
h( k, j ) = h( k, j ) - sum*t1
430
h( k+1, j ) = h( k+1, j ) - sum*t2
431
90 continue
432
c
433
c %----------------------------------------------------%
434
c | Apply G from the right to transform the columns of |
435
c | the matrix in rows I1 to min(K+3,I). |
436
c %----------------------------------------------------%
437
c
438
do 100 j = i1, i
439
sum = h( j, k ) + v2*h( j, k+1 )
440
h( j, k ) = h( j, k ) - sum*t1
441
h( j, k+1 ) = h( j, k+1 ) - sum*t2
442
100 continue
443
c
444
c %----------------------------------%
445
c | Accumulate transformations for Z |
446
c %----------------------------------%
447
c
448
sum = z( k ) + v2*z( k+1 )
449
z( k ) = z( k ) - sum*t1
450
z( k+1 ) = z( k+1 ) - sum*t2
451
end if
452
120 continue
453
454
130 continue
455
c
456
c %-------------------------------------------------------%
457
c | Failure to converge in remaining number of iterations |
458
c %-------------------------------------------------------%
459
c
460
info = i
461
return
462
463
140 continue
464
465
if( l.eq.i ) then
466
c
467
c %------------------------------------------------------%
468
c | H(I,I-1) is negligible: one eigenvalue has converged |
469
c %------------------------------------------------------%
470
c
471
wr( i ) = h( i, i )
472
wi( i ) = zero
473
474
else if( l.eq.i-1 ) then
475
c
476
c %--------------------------------------------------------%
477
c | H(I-1,I-2) is negligible; |
478
c | a pair of eigenvalues have converged. |
479
c | |
480
c | Transform the 2-by-2 submatrix to standard Schur form, |
481
c | and compute and store the eigenvalues. |
482
c %--------------------------------------------------------%
483
c
484
call slanv2( h( i-1, i-1 ), h( i-1, i ), h( i, i-1 ),
485
& h( i, i ), wr( i-1 ), wi( i-1 ), wr( i ), wi( i ),
486
& cs, sn )
487
488
if( wantt ) then
489
c
490
c %-----------------------------------------------------%
491
c | Apply the transformation to the rest of H and to Z, |
492
c | as required. |
493
c %-----------------------------------------------------%
494
c
495
if( i2.gt.i )
496
& call srot( i2-i, h( i-1, i+1 ), ldh, h( i, i+1 ), ldh,
497
& cs, sn )
498
call srot( i-i1-1, h( i1, i-1 ), 1, h( i1, i ), 1, cs, sn )
499
sum = cs*z( i-1 ) + sn*z( i )
500
z( i ) = cs*z( i ) - sn*z( i-1 )
501
z( i-1 ) = sum
502
end if
503
end if
504
c
505
c %---------------------------------------------------------%
506
c | Decrement number of remaining iterations, and return to |
507
c | start of the main loop with new value of I. |
508
c %---------------------------------------------------------%
509
c
510
itn = itn - its
511
i = l - 1
512
go to 10
513
514
150 continue
515
return
516
c
517
c %---------------%
518
c | End of slaqrb |
519
c %---------------%
520
c
521
end
Loggerhead is a web-based interface for Breezy
Version: 2.0.1