← Back to branch summary

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

« back to all changes in this revision

Viewing changes to routines/control/comqr3.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/comqr3.f
 
1
      subroutine comqr3(nm,n,low,igh,hr,hi,wr,wi,zr,zi,ierr,job)
 
2
c
 
3
      integer i,j,l,n,en,ll,nm,igh,ip1,
 
4
     x        itn,its,low,lp1,enm1,iend,ierr
 
5
      double precision hr(nm,n),hi(nm,n),wr(n),wi(n),zr(nm,n),zi(nm,n)
 
6
      double precision si,sr,ti,tr,xi,xr,yi,yr,zzi,zzr,norm
 
7
      double precision pythag
 
8
c
 
9
c!originator
 
10
c     this subroutine is a translation of a unitary analogue of the
 
11
c     algol procedure  comlr2, num. math. 16, 181-204(1970) by peters
 
12
c     and wilkinson.
 
13
c     handbook for auto. comp., vol.ii-linear algebra, 372-395(1971).
 
14
c     the unitary analogue substitutes the qr algorithm of francis
 
15
c     (comp. jour. 4, 332-345(1962)) for the lr algorithm.
 
16
c
 
17
c     modified by  c. moler
 
18
c!purpose
 
19
c     this subroutine finds the eigenvalues of a complex upper 
 
20
c     hessenberg matrix by the qr method. The unitary transformation
 
21
c     can also be accumulated if  corth  has been used to reduce
 
22
c     this general matrix to hessenberg form.
 
23
c
 
24
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
25
c     MODIFICATION OF EISPACK COMQR+COMQR2 
 
26
c        1. job parameter added 
 
27
c        2. code concerning eigenvector computation deleted
 
28
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
29
c
 
30
c!calling sequence
 
31
c      subroutine comqr3(nm,n,low,igh,hr,hi,wr,wi,zr,zi,ierr
 
32
c     *                 ,job)
 
33
c
 
34
c     on input.
 
35
c
 
36
c        nm must be set to the row dimension of two-dimensional
 
37
c          array parameters as declared in the calling program
 
38
c          dimension statement.
 
39
c
 
40
c        n is the order of the matrix.
 
41
c
 
42
c        low and igh are integers determined by the balancing
 
43
c          subroutine  cbal.  if  cbal  has not been used,
 
44
c          set low=1, igh=n.
 
45
c
 
46
c        hr and hi contain the real and imaginary parts,
 
47
c          respectively, of the complex upper hessenberg matrix.
 
48
c          their lower triangles below the subdiagonal contain further
 
49
c          information about the transformations which were used in the
 
50
c          reduction by  corth, if performed.  
 
51
c
 
52
c        zr and zi contain the real and imaginary parts,respectively
 
53
c           of the unitary similarity used to put h on hessenberg form
 
54
c           or a unitary matrix ,if vectors are desired
 
55
c
 
56
c       job indicate the job to be performed: job=xy
 
57
c           if y=0 no accumulation of the unitary transformation
 
58
c           if y=1 transformation  accumulated in z
 
59
c
 
60
c     on output.
 
61
c     the upper hessenberg portions of hr and hi have been destroyed
 
62
c
 
63
c
 
64
c        wr and wi contain the real and imaginary parts,
 
65
c          respectively, of the eigenvalues.  if an error
 
66
c          exit is made, the eigenvalues should be correct
 
67
c          for indices ierr+1,...,n.
 
68
c
 
69
c        zr and zi contain the real and imaginary parts,
 
70
c          respectively, of the eigenvectors.  the eigenvectors
 
71
c          are unnormalized.  if an error exit is made, none of
 
72
c          the eigenvectors has been found.
 
73
c
 
74
c        ierr is set to
 
75
c          zero       for normal return,
 
76
c          j          if the j-th eigenvalue has not been
 
77
c                     determined after a total of 30*n iterations.
 
78
 
 
79
c
 
80
c     questions and comments should be directed to b. s. garbow,
 
81
c     applied mathematics division, argonne national laboratory
 
82
c
 
83
c!auxiliary routines
 
84
c     pythag
 
85
c!
 
86
c     ------------------------------------------------------------------
 
87
c
 
88
      ierr = 0
 
89
c*****
 
90
      jx=job/10
 
91
      jy=job-10*jx
 
92
c
 
93
c     .......... create real subdiagonal elements ..........
 
94
      iend=igh-low-1
 
95
      if(iend.lt.0) goto 180
 
96
  150 l = low + 1
 
97
c
 
98
      do 170 i = l, igh
 
99
         ll = min(i+1,igh)
 
100
         if (hi(i,i-1) .eq. 0.0d+0) go to 170
 
101
         norm = pythag(hr(i,i-1),hi(i,i-1))
 
102
         yr = hr(i,i-1)/norm
 
103
         yi = hi(i,i-1)/norm
 
104
         hr(i,i-1) = norm
 
105
         hi(i,i-1) = 0.0d+0
 
106
c
 
107
         do 155 j = i, n
 
108
            si = yr*hi(i,j) - yi*hr(i,j)
 
109
            hr(i,j) = yr*hr(i,j) + yi*hi(i,j)
 
110
            hi(i,j) = si
 
111
  155    continue
 
112
c
 
113
         do 160 j = 1, ll
 
114
            si = yr*hi(j,i) + yi*hr(j,i)
 
115
            hr(j,i) = yr*hr(j,i) - yi*hi(j,i)
 
116
            hi(j,i) = si
 
117
  160    continue
 
118
c*****
 
119
         if (jy .eq. 0) go to 170
 
120
c*****
 
121
         do 165 j = low, igh
 
122
            si = yr*zi(j,i) + yi*zr(j,i)
 
123
            zr(j,i) = yr*zr(j,i) - yi*zi(j,i)
 
124
            zi(j,i) = si
 
125
  165    continue
 
126
c
 
127
  170 continue
 
128
c     .......... store roots isolated by cbal ..........
 
129
c
 
130
  180 do 200 i = 1, n
 
131
         if (i .ge. low .and. i .le. igh) go to 200
 
132
         wr(i) = hr(i,i)
 
133
         wi(i) = hi(i,i)
 
134
  200 continue
 
135
c
 
136
  210 continue
 
137
      en = igh
 
138
      tr = 0.0d+0
 
139
      ti = 0.0d+0
 
140
      itn = 30*n
 
141
c     .......... search for next eigenvalue ..........
 
142
  220 if (en .lt. low) go to 1001
 
143
      its = 0
 
144
      enm1 = en - 1
 
145
c     .......... look for single small sub-diagonal element
 
146
c                for l=en step -1 until low do -- ..........
 
147
  240 do 260 ll = low, en
 
148
         l = en + low - ll
 
149
         if (l .eq. low) go to 300
 
150
c*****
 
151
         xr = abs(hr(l-1,l-1)) + abs(hi(l-1,l-1)
 
152
     x             + abs(hr(l,l)) +abs(hi(l,l)))
 
153
         yr = xr + abs(hr(l,l-1))
 
154
         if (xr .eq. yr) go to 300
 
155
c*****
 
156
  260 continue
 
157
c     .......... form shift ..........
 
158
  300 if (l .eq. en) go to 660
 
159
      if (itn .eq. 0) go to 1000
 
160
      if (its .eq. 10 .or. its .eq. 20) go to 320
 
161
      sr = hr(en,en)
 
162
      si = hi(en,en)
 
163
      xr = hr(enm1,en)*hr(en,enm1)
 
164
      xi = hi(enm1,en)*hr(en,enm1)
 
165
      if (xr .eq. 0.0d+0 .and. xi .eq. 0.0d+0) go to 340
 
166
      yr = (hr(enm1,enm1) - sr)/2.0d+0
 
167
      yi = (hi(enm1,enm1) - si)/2.0d+0
 
168
      call wsqrt(yr**2-yi**2+xr,2.0d+0*yr*yi+xi,zzr,zzi)
 
169
      if (yr*zzr + yi*zzi .ge. 0.0d+0) go to 310
 
170
      zzr = -zzr
 
171
      zzi = -zzi
 
172
  310 call cdiv(xr,xi,yr+zzr,yi+zzi,zzr,zzi)
 
173
      sr = sr - zzr
 
174
      si = si - zzi
 
175
      go to 340
 
176
c     .......... form exceptional shift ..........
 
177
  320 sr = abs(hr(en,enm1)) + abs(hr(enm1,en-2))
 
178
      si = 0.0d+0
 
179
c
 
180
  340 do 360 i = low, en
 
181
         hr(i,i) = hr(i,i) - sr
 
182
         hi(i,i) = hi(i,i) - si
 
183
  360 continue
 
184
c
 
185
      tr = tr + sr
 
186
      ti = ti + si
 
187
      its = its + 1
 
188
      itn = itn - 1
 
189
c     .......... reduce to triangle (rows) ..........
 
190
      lp1 = l + 1
 
191
c
 
192
      do 500 i = lp1, en
 
193
         sr = hr(i,i-1)
 
194
         hr(i,i-1) = 0.0d+0
 
195
         norm = pythag(pythag(hr(i-1,i-1),hi(i-1,i-1)),sr)
 
196
         xr = hr(i-1,i-1)/norm
 
197
         wr(i-1) = xr
 
198
         xi = hi(i-1,i-1)/norm
 
199
         wi(i-1) = xi
 
200
         hr(i-1,i-1) = norm
 
201
         hi(i-1,i-1) = 0.0d+0
 
202
         hi(i,i-1) = sr/norm
 
203
c
 
204
         do 490 j = i, n
 
205
            yr = hr(i-1,j)
 
206
            yi = hi(i-1,j)
 
207
            zzr = hr(i,j)
 
208
            zzi = hi(i,j)
 
209
            hr(i-1,j) = xr*yr + xi*yi + hi(i,i-1)*zzr
 
210
            hi(i-1,j) = xr*yi - xi*yr + hi(i,i-1)*zzi
 
211
            hr(i,j) = xr*zzr - xi*zzi - hi(i,i-1)*yr
 
212
            hi(i,j) = xr*zzi + xi*zzr - hi(i,i-1)*yi
 
213
  490    continue
 
214
c
 
215
  500 continue
 
216
c
 
217
      si = hi(en,en)
 
218
      if (si .eq. 0.0d+0) go to 540
 
219
      norm = pythag(hr(en,en),si)
 
220
      sr = hr(en,en)/norm
 
221
      si = si/norm
 
222
      hr(en,en) = norm
 
223
      hi(en,en) = 0.0d+0
 
224
      if (en .eq. n) go to 540
 
225
      ip1 = en + 1
 
226
c
 
227
      do 520 j = ip1, n
 
228
         yr = hr(en,j)
 
229
         yi = hi(en,j)
 
230
         hr(en,j) = sr*yr + si*yi
 
231
         hi(en,j) = sr*yi - si*yr
 
232
  520 continue
 
233
c     .......... inverse operation (columns) ..........
 
234
  540 do 600 j = lp1, en
 
235
         xr = wr(j-1)
 
236
         xi = wi(j-1)
 
237
c
 
238
         do 580 i = 1, j
 
239
            yr = hr(i,j-1)
 
240
            yi = 0.0d+0
 
241
            zzr = hr(i,j)
 
242
            zzi = hi(i,j)
 
243
            if (i .eq. j) go to 560
 
244
            yi = hi(i,j-1)
 
245
            hi(i,j-1) = xr*yi + xi*yr + hi(j,j-1)*zzi
 
246
  560       hr(i,j-1) = xr*yr - xi*yi + hi(j,j-1)*zzr
 
247
            hr(i,j) = xr*zzr + xi*zzi - hi(j,j-1)*yr
 
248
            hi(i,j) = xr*zzi - xi*zzr - hi(j,j-1)*yi
 
249
  580    continue
 
250
c*****
 
251
         if (jy .eq. 0) go to 600
 
252
c*****
 
253
         do 590 i = low, igh
 
254
            yr = zr(i,j-1)
 
255
            yi = zi(i,j-1)
 
256
            zzr = zr(i,j)
 
257
            zzi = zi(i,j)
 
258
            zr(i,j-1) = xr*yr - xi*yi + hi(j,j-1)*zzr
 
259
            zi(i,j-1) = xr*yi + xi*yr + hi(j,j-1)*zzi
 
260
            zr(i,j) = xr*zzr + xi*zzi - hi(j,j-1)*yr
 
261
            zi(i,j) = xr*zzi - xi*zzr - hi(j,j-1)*yi
 
262
  590    continue
 
263
c
 
264
  600 continue
 
265
c
 
266
      if (si .eq. 0.0d+0) go to 240
 
267
c
 
268
      do 630 i = 1, en
 
269
         yr = hr(i,en)
 
270
         yi = hi(i,en)
 
271
         hr(i,en) = sr*yr - si*yi
 
272
         hi(i,en) = sr*yi + si*yr
 
273
  630 continue
 
274
c*****
 
275
      if (jy .eq. 0) go to 240
 
276
c*****
 
277
      do 640 i = low, igh
 
278
         yr = zr(i,en)
 
279
         yi = zi(i,en)
 
280
         zr(i,en) = sr*yr - si*yi
 
281
         zi(i,en) = sr*yi + si*yr
 
282
  640 continue
 
283
c
 
284
      go to 240
 
285
c     .......... a root found ..........
 
286
  660 hr(en,en) = hr(en,en) + tr
 
287
      wr(en) = hr(en,en)
 
288
      hi(en,en) = hi(en,en) + ti
 
289
      wi(en) = hi(en,en)
 
290
      en = enm1
 
291
      go to 220
 
292
c     .......... all roots found.  ..........
 
293
 
 
294
c      go to 1001
 
295
c
 
296
c     .......... set error -- no convergence to an
 
297
c                eigenvalue after 30 iterations ..........
 
298
 1000 ierr = en
 
299
 1001 return
 
300
      end