← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/control/scaleg.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/control/scaleg.f
 
1
      subroutine scaleg (n,ma,a,mb,b,low,igh,cscale,cperm,wk)
 
2
c
 
3
c     *****parameters:
 
4
      integer igh,low,ma,mb,n
 
5
      double precision a(ma,n),b(mb,n),cperm(n),cscale(n),wk(n,6)
 
6
c
 
7
c     *****local variables:
 
8
      integer i,ir,it,j,jc,kount,nr,nrp2
 
9
      double precision alpha,basl,beta,cmax,coef,coef2,coef5,cor,
 
10
     *                 ew,ewc,fi,fj,gamma,pgamma,sum,t,ta,tb,tc
 
11
c
 
12
c     *****fortran functions:
 
13
      double precision dabs, dlog10, dsign
 
14
c     float
 
15
c
 
16
c     *****subroutines called:
 
17
c     none
 
18
c
 
19
c     ---------------------------------------------------------------
 
20
c
 
21
c     *****purpose:
 
22
c     scales the matrices a and b in the generalized eigenvalue
 
23
c     problem a*x = (lambda)*b*x such that the magnitudes of the
 
24
c     elements of the submatrices of a and b (as specified by low
 
25
c     and igh) are close to unity in the least squares sense.
 
26
c     ref.:  ward, r. c., balancing the generalized eigenvalue
 
27
c     problem, siam j. sci. stat. comput., vol. 2, no. 2, june 1981,
 
28
c     141-152.
 
29
c
 
30
c     *****parameter description:
 
31
c
 
32
c     on input:
 
33
c
 
34
c       ma,mb   integer
 
35
c               row dimensions of the arrays containing matrices
 
36
c               a and b respectively, as declared in the main calling
 
37
c               program dimension statement;
 
38
c
 
39
c       n       integer
 
40
c               order of the matrices a and b;
 
41
c
 
42
c       a       real(ma,n)
 
43
c               contains the a matrix of the generalized eigenproblem
 
44
c               defined above;
 
45
c
 
46
c       b       real(mb,n)
 
47
c               contains the b matrix of the generalized eigenproblem
 
48
c               defined above;
 
49
c
 
50
c       low     integer
 
51
c               specifies the beginning -1 for the rows and
 
52
c               columns of a and b to be scaled;
 
53
c
 
54
c       igh     integer
 
55
c               specifies the ending -1 for the rows and columns
 
56
c               of a and b to be scaled;
 
57
c
 
58
c       cperm   real(n)
 
59
c               work array.  only locations low through igh are
 
60
c               referenced and altered by this subroutine;
 
61
c
 
62
c       wk      real(n,6)
 
63
c               work array that must contain at least 6*n locations.
 
64
c               only locations low through igh, n+low through n+igh,
 
65
c               ..., 5*n+low through 5*n+igh are referenced and
 
66
c               altered by this subroutine.
 
67
c
 
68
c     on output:
 
69
c
 
70
c       a,b     contain the scaled a and b matrices;
 
71
c
 
72
c       cscale  real(n)
 
73
c               contains in its low through igh locations the integer
 
74
c               exponents of 2 used for the column scaling factors.
 
75
c               the other locations are not referenced;
 
76
c
 
77
c       wk      contains in its low through igh locations the integer
 
78
c               exponents of 2 used for the row scaling factors.
 
79
c
 
80
c     *****algorithm notes:
 
81
c     none.
 
82
c
 
83
c     *****history:
 
84
c     written by r. c. ward.......
 
85
c     modified 8/86 by bobby bodenheimer so that if
 
86
c       sum = 0 (corresponding to the case where the matrix
 
87
c       doesn't need to be scaled) the routine returns.
 
88
c
 
89
c     ---------------------------------------------------------------
 
90
c
 
91
      if (low .eq. igh) go to 410
 
92
      do 210 i = low,igh
 
93
         wk(i,1) = 0.0d0
 
94
         wk(i,2) = 0.0d0
 
95
         wk(i,3) = 0.0d0
 
96
         wk(i,4) = 0.0d0
 
97
         wk(i,5) = 0.0d0
 
98
         wk(i,6) = 0.0d0
 
99
         cscale(i) = 0.0d0
 
100
         cperm(i) = 0.0d0
 
101
  210 continue
 
102
c
 
103
c     compute right side vector in resulting linear equations
 
104
c
 
105
      basl = dlog10(2.0d0)
 
106
      do 240 i = low,igh
 
107
         do 240 j = low,igh
 
108
            tb = b(i,j)
 
109
            ta = a(i,j)
 
110
            if (ta .eq. 0.0d0) go to 220
 
111
            ta = dlog10(dabs(ta)) / basl
 
112
  220       continue
 
113
            if (tb .eq. 0.0d0) go to 230
 
114
            tb = dlog10(dabs(tb)) / basl
 
115
  230       continue
 
116
            wk(i,5) = wk(i,5) - ta - tb
 
117
            wk(j,6) = wk(j,6) - ta - tb
 
118
  240 continue
 
119
      nr = igh-low+1
 
120
      coef = 1.0d0/float(2*nr)
 
121
      coef2 = coef*coef
 
122
      coef5 = 0.5d0*coef2
 
123
      nrp2 = nr+2
 
124
      beta = 0.0d0
 
125
      it = 1
 
126
c
 
127
c     start generalized conjugate gradient iteration
 
128
c
 
129
  250 continue
 
130
      ew = 0.0d0
 
131
      ewc = 0.0d0
 
132
      gamma = 0.0d0
 
133
      do 260 i = low,igh
 
134
         gamma = gamma + wk(i,5)*wk(i,5) + wk(i,6)*wk(i,6)
 
135
         ew = ew + wk(i,5)
 
136
         ewc = ewc + wk(i,6)
 
137
  260 continue
 
138
      gamma = coef*gamma - coef2*(ew**2 + ewc**2)
 
139
     +        - coef5*(ew - ewc)**2
 
140
      if (it .ne. 1) beta = gamma / pgamma
 
141
      t = coef5*(ewc - 3.0d0*ew)
 
142
      tc = coef5*(ew - 3.0d0*ewc)
 
143
      do 270 i = low,igh
 
144
         wk(i,2) = beta*wk(i,2) + coef*wk(i,5) + t
 
145
         cperm(i) = beta*cperm(i) + coef*wk(i,6) + tc
 
146
  270 continue
 
147
c
 
148
c     apply matrix to vector
 
149
c
 
150
      do 300 i = low,igh
 
151
         kount = 0
 
152
         sum = 0.0d0
 
153
         do 290 j = low,igh
 
154
            if (a(i,j) .eq. 0.0d0) go to 280
 
155
            kount = kount+1
 
156
            sum = sum + cperm(j)
 
157
  280       continue
 
158
            if (b(i,j) .eq. 0.0d0) go to 290
 
159
            kount = kount+1
 
160
            sum = sum + cperm(j)
 
161
  290    continue
 
162
         wk(i,3) = float(kount)*wk(i,2) + sum
 
163
  300 continue
 
164
      do 330 j = low,igh
 
165
         kount = 0
 
166
         sum = 0.0d0
 
167
         do 320 i = low,igh
 
168
            if (a(i,j) .eq. 0.0d0) go to 310
 
169
            kount = kount+1
 
170
            sum = sum + wk(i,2)
 
171
  310       continue
 
172
            if (b(i,j) .eq. 0.0d0) go to 320
 
173
            kount = kount+1
 
174
            sum = sum + wk(i,2)
 
175
  320    continue
 
176
         wk(j,4) = float(kount)*cperm(j) + sum
 
177
  330 continue
 
178
      sum = 0.0d0
 
179
      do 340 i = low,igh
 
180
         sum = sum + wk(i,2)*wk(i,3) + cperm(i)*wk(i,4)
 
181
  340 continue
 
182
      if(sum.eq.0.0d0) return
 
183
      alpha = gamma / sum
 
184
c
 
185
c     determine correction to current iterate
 
186
c
 
187
      cmax = 0.0d0
 
188
      do 350 i = low,igh
 
189
         cor = alpha * wk(i,2)
 
190
         if (dabs(cor) .gt. cmax) cmax = dabs(cor)
 
191
         wk(i,1) = wk(i,1) + cor
 
192
         cor = alpha * cperm(i)
 
193
         if (dabs(cor) .gt. cmax) cmax = dabs(cor)
 
194
         cscale(i) = cscale(i) + cor
 
195
  350 continue
 
196
      if (cmax .lt. 0.5d0) go to 370
 
197
      do 360 i = low,igh
 
198
         wk(i,5) = wk(i,5) - alpha*wk(i,3)
 
199
         wk(i,6) = wk(i,6) - alpha*wk(i,4)
 
200
  360 continue
 
201
      pgamma = gamma
 
202
      it = it+1
 
203
      if (it .le. nrp2) go to 250
 
204
c
 
205
c     end generalized conjugate gradient iteration
 
206
c
 
207
  370 continue
 
208
      do 380 i = low,igh
 
209
         ir = wk(i,1) + dsign(0.5d0,wk(i,1))
 
210
         wk(i,1) = ir
 
211
         jc = cscale(i) + dsign(0.5d0,cscale(i))
 
212
         cscale(i) = jc
 
213
  380 continue
 
214
c
 
215
c     scale a and b
 
216
c
 
217
      do 400 i = 1,igh
 
218
         ir = wk(i,1)
 
219
         fi = 2.0d0**ir
 
220
         if (i .lt. low) fi = 1.0d0
 
221
         do 400 j =low,n
 
222
            jc = cscale(j)
 
223
            fj = 2.0d0**jc
 
224
            if (j .le. igh) go to 390
 
225
            if (i .lt. low) go to 400
 
226
            fj = 1.0d0
 
227
  390       continue
 
228
            a(i,j) = a(i,j)*fi*fj
 
229
            b(i,j) = b(i,j)*fi*fj
 
230
  400 continue
 
231
  410 continue
 
232
      return
 
233
c
 
234
c     last line of scaleg
 
235
c
 
236
      end