1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
|
SUBROUTINE SLARRF( N, D, L, LD, LLD, IFIRST, ILAST, W, DPLUS,
$ LPLUS, WORK, IWORK, INFO )
*
* -- LAPACK auxiliary routine (instru to count ops, version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* June 30, 1999
*
* .. Scalar Arguments ..
INTEGER IFIRST, ILAST, INFO, N
* ..
* .. Array Arguments ..
INTEGER IWORK( * )
REAL D( * ), DPLUS( * ), L( * ), LD( * ), LLD( * ),
$ LPLUS( * ), W( * ), WORK( * )
* ..
* Common block to return operation count
* .. Common blocks ..
COMMON / LATIME / OPS, ITCNT
* ..
* .. Scalars in Common ..
REAL ITCNT, OPS
* ..
*
* Purpose
* =======
*
* Given the initial representation L D L^T and its cluster of close
* eigenvalues (in a relative measure), W( IFIRST ), W( IFIRST+1 ), ...
* W( ILAST ), SLARRF finds a new relatively robust representation
* L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the
* eigenvalues of L(+) D(+) L(+)^T is relatively isolated.
*
* Arguments
* =========
*
* N (input) INTEGER
* The order of the matrix.
*
* D (input) REAL array, dimension (N)
* The n diagonal elements of the diagonal matrix D.
*
* L (input) REAL array, dimension (N-1)
* The (n-1) subdiagonal elements of the unit bidiagonal
* matrix L.
*
* LD (input) REAL array, dimension (N-1)
* The n-1 elements L(i)*D(i).
*
* LLD (input) REAL array, dimension (N-1)
* The n-1 elements L(i)*L(i)*D(i).
*
* IFIRST (input) INTEGER
* The index of the first eigenvalue in the cluster.
*
* ILAST (input) INTEGER
* The index of the last eigenvalue in the cluster.
*
* W (input/output) REAL array, dimension (N)
* On input, the eigenvalues of L D L^T in ascending order.
* W( IFIRST ) through W( ILAST ) form the cluster of relatively
* close eigenalues.
* On output, W( IFIRST ) thru' W( ILAST ) are estimates of the
* corresponding eigenvalues of L(+) D(+) L(+)^T.
*
* SIGMA (input) REAL
* The shift used to form L(+) D(+) L(+)^T.
*
* DPLUS (output) REAL array, dimension (N)
* The n diagonal elements of the diagonal matrix D(+).
*
* LPLUS (output) REAL array, dimension (N)
* The first (n-1) elements of LPLUS contain the subdiagonal
* elements of the unit bidiagonal matrix L(+). LPLUS( N ) is
* set to SIGMA.
*
* WORK (input) REAL array, dimension (???)
* Workspace.
*
* Further Details
* ===============
*
* Based on contributions by
* Inderjit Dhillon, IBM Almaden, USA
* Osni Marques, LBNL/NERSC, USA
*
* =====================================================================
*
* .. Parameters ..
REAL ZERO, TWO
PARAMETER ( ZERO = 0.0E0, TWO = 2.0E0 )
* ..
* .. Local Scalars ..
INTEGER I
REAL DELTA, EPS, S, SIGMA
* ..
* .. External Functions ..
REAL SLAMCH
EXTERNAL SLAMCH
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, REAL
* ..
* .. Executable Statements ..
*
INFO = 0
EPS = SLAMCH( 'Precision' )
IF( IFIRST.EQ.1 ) THEN
SIGMA = W( IFIRST )
ELSE IF( ILAST.EQ.N ) THEN
SIGMA = W( ILAST )
ELSE
INFO = 1
RETURN
END IF
*
* Compute the new relatively robust representation (RRR)
*
OPS = OPS + REAL( 3 )
DELTA = TWO*EPS
10 CONTINUE
IF( IFIRST.EQ.1 ) THEN
SIGMA = SIGMA - ABS( SIGMA )*DELTA
ELSE
SIGMA = SIGMA + ABS( SIGMA )*DELTA
END IF
S = -SIGMA
OPS = OPS + REAL( 5*(N-1)+1 )
DO 20 I = 1, N - 1
DPLUS( I ) = D( I ) + S
LPLUS( I ) = LD( I ) / DPLUS( I )
S = S*LPLUS( I )*L( I ) - SIGMA
20 CONTINUE
DPLUS( N ) = D( N ) + S
IF( IFIRST.EQ.1 ) THEN
DO 30 I = 1, N
IF( DPLUS( I ).LT.ZERO ) THEN
OPS = OPS + REAL( 1 )
DELTA = TWO*DELTA
GO TO 10
END IF
30 CONTINUE
ELSE
DO 40 I = 1, N
IF( DPLUS( I ).GT.ZERO ) THEN
OPS = OPS + REAL( 1 )
DELTA = TWO*DELTA
GO TO 10
END IF
40 CONTINUE
END IF
DO 50 I = IFIRST, ILAST
OPS = OPS + REAL( 1 )
W( I ) = W( I ) - SIGMA
50 CONTINUE
LPLUS( N ) = SIGMA
*
RETURN
*
* End of SLARRF
*
END
|
