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
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
|
subroutine scaleg (n,ma,a,mb,b,low,igh,cscale,cperm,wk)
c
c *****parameters:
integer igh,low,ma,mb,n
double precision a(ma,n),b(mb,n),cperm(n),cscale(n),wk(n,6)
c
c *****local variables:
integer i,ir,it,j,jc,kount,nr,nrp2
double precision alpha,basl,beta,cmax,coef,coef2,coef5,cor,
* ew,ewc,fi,fj,gamma,pgamma,sum,t,ta,tb,tc
c
c *****fortran functions:
double precision dabs, dlog10, dsign
c float
c
c *****subroutines called:
c none
c
c ---------------------------------------------------------------
c
c *****purpose:
c scales the matrices a and b in the generalized eigenvalue
c problem a*x = (lambda)*b*x such that the magnitudes of the
c elements of the submatrices of a and b (as specified by low
c and igh) are close to unity in the least squares sense.
c ref.: ward, r. c., balancing the generalized eigenvalue
c problem, siam j. sci. stat. comput., vol. 2, no. 2, june 1981,
c 141-152.
c
c *****parameter description:
c
c on input:
c
c ma,mb integer
c row dimensions of the arrays containing matrices
c a and b respectively, as declared in the main calling
c program dimension statement;
c
c n integer
c order of the matrices a and b;
c
c a real(ma,n)
c contains the a matrix of the generalized eigenproblem
c defined above;
c
c b real(mb,n)
c contains the b matrix of the generalized eigenproblem
c defined above;
c
c low integer
c specifies the beginning -1 for the rows and
c columns of a and b to be scaled;
c
c igh integer
c specifies the ending -1 for the rows and columns
c of a and b to be scaled;
c
c cperm real(n)
c work array. only locations low through igh are
c referenced and altered by this subroutine;
c
c wk real(n,6)
c work array that must contain at least 6*n locations.
c only locations low through igh, n+low through n+igh,
c ..., 5*n+low through 5*n+igh are referenced and
c altered by this subroutine.
c
c on output:
c
c a,b contain the scaled a and b matrices;
c
c cscale real(n)
c contains in its low through igh locations the integer
c exponents of 2 used for the column scaling factors.
c the other locations are not referenced;
c
c wk contains in its low through igh locations the integer
c exponents of 2 used for the row scaling factors.
c
c *****algorithm notes:
c none.
c
c *****history:
c written by r. c. ward.......
c modified 8/86 by bobby bodenheimer so that if
c sum = 0 (corresponding to the case where the matrix
c doesn't need to be scaled) the routine returns.
c
c ---------------------------------------------------------------
c
if (low .eq. igh) go to 410
do 210 i = low,igh
wk(i,1) = 0.0d0
wk(i,2) = 0.0d0
wk(i,3) = 0.0d0
wk(i,4) = 0.0d0
wk(i,5) = 0.0d0
wk(i,6) = 0.0d0
cscale(i) = 0.0d0
cperm(i) = 0.0d0
210 continue
c
c compute right side vector in resulting linear equations
c
basl = dlog10(2.0d0)
do 240 i = low,igh
do 240 j = low,igh
tb = b(i,j)
ta = a(i,j)
if (ta .eq. 0.0d0) go to 220
ta = dlog10(dabs(ta)) / basl
220 continue
if (tb .eq. 0.0d0) go to 230
tb = dlog10(dabs(tb)) / basl
230 continue
wk(i,5) = wk(i,5) - ta - tb
wk(j,6) = wk(j,6) - ta - tb
240 continue
nr = igh-low+1
coef = 1.0d0/float(2*nr)
coef2 = coef*coef
coef5 = 0.5d0*coef2
nrp2 = nr+2
beta = 0.0d0
it = 1
c
c start generalized conjugate gradient iteration
c
250 continue
ew = 0.0d0
ewc = 0.0d0
gamma = 0.0d0
do 260 i = low,igh
gamma = gamma + wk(i,5)*wk(i,5) + wk(i,6)*wk(i,6)
ew = ew + wk(i,5)
ewc = ewc + wk(i,6)
260 continue
gamma = coef*gamma - coef2*(ew**2 + ewc**2)
+ - coef5*(ew - ewc)**2
if (it .ne. 1) beta = gamma / pgamma
t = coef5*(ewc - 3.0d0*ew)
tc = coef5*(ew - 3.0d0*ewc)
do 270 i = low,igh
wk(i,2) = beta*wk(i,2) + coef*wk(i,5) + t
cperm(i) = beta*cperm(i) + coef*wk(i,6) + tc
270 continue
c
c apply matrix to vector
c
do 300 i = low,igh
kount = 0
sum = 0.0d0
do 290 j = low,igh
if (a(i,j) .eq. 0.0d0) go to 280
kount = kount+1
sum = sum + cperm(j)
280 continue
if (b(i,j) .eq. 0.0d0) go to 290
kount = kount+1
sum = sum + cperm(j)
290 continue
wk(i,3) = float(kount)*wk(i,2) + sum
300 continue
do 330 j = low,igh
kount = 0
sum = 0.0d0
do 320 i = low,igh
if (a(i,j) .eq. 0.0d0) go to 310
kount = kount+1
sum = sum + wk(i,2)
310 continue
if (b(i,j) .eq. 0.0d0) go to 320
kount = kount+1
sum = sum + wk(i,2)
320 continue
wk(j,4) = float(kount)*cperm(j) + sum
330 continue
sum = 0.0d0
do 340 i = low,igh
sum = sum + wk(i,2)*wk(i,3) + cperm(i)*wk(i,4)
340 continue
if(sum.eq.0.0d0) return
alpha = gamma / sum
c
c determine correction to current iterate
c
cmax = 0.0d0
do 350 i = low,igh
cor = alpha * wk(i,2)
if (dabs(cor) .gt. cmax) cmax = dabs(cor)
wk(i,1) = wk(i,1) + cor
cor = alpha * cperm(i)
if (dabs(cor) .gt. cmax) cmax = dabs(cor)
cscale(i) = cscale(i) + cor
350 continue
if (cmax .lt. 0.5d0) go to 370
do 360 i = low,igh
wk(i,5) = wk(i,5) - alpha*wk(i,3)
wk(i,6) = wk(i,6) - alpha*wk(i,4)
360 continue
pgamma = gamma
it = it+1
if (it .le. nrp2) go to 250
c
c end generalized conjugate gradient iteration
c
370 continue
do 380 i = low,igh
ir = wk(i,1) + dsign(0.5d0,wk(i,1))
wk(i,1) = ir
jc = cscale(i) + dsign(0.5d0,cscale(i))
cscale(i) = jc
380 continue
c
c scale a and b
c
do 400 i = 1,igh
ir = wk(i,1)
fi = 2.0d0**ir
if (i .lt. low) fi = 1.0d0
do 400 j =low,n
jc = cscale(j)
fj = 2.0d0**jc
if (j .le. igh) go to 390
if (i .lt. low) go to 400
fj = 1.0d0
390 continue
a(i,j) = a(i,j)*fi*fj
b(i,j) = b(i,j)*fi*fj
400 continue
410 continue
return
c
c last line of scaleg
c
end
|
