← Back to branch summary

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

« back to all changes in this revision

Viewing changes to routines/control/qvecz.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/qvecz.f
 
1
C/MEMBR ADD NAME=QVECZ,SSI=0
 
2
      subroutine qvecz(nm,n,a,b,epsb,alfr,alfi,beta,z)
 
3
c
 
4
      integer i,j,k,m,n,en,ii,jj,na,nm,nn,isw,enm2
 
5
      double precision a(nm,n),b(nm,n),alfr(n),alfi(n),beta(n)
 
6
      double precision z(nm,n)
 
7
      double precision d,q,r,s,t,w,x,y,di,dr,ra,rr,sa,ti,tr,t1,t2,w1,
 
8
     x       x1,zz,z1,alfm,almi,almr,betm,epsb
 
9
c! purpose
 
10
c
 
11
c     this subroutine accepts a pair of real matrices, one of them in
 
12
c     quasi-triangular form (in which each 2-by-2 block corresponds to
 
13
c     a pair of complex eigenvalues) and the other in upper triangular
 
14
c     form.  it computes the eigenvectors of the triangular problem and
 
15
c     transforms the results back to the original coordinate system.
 
16
c     it is usually preceded by  qzhes,  qzit, and  qzval.
 
17
c! calling sequence
 
18
c
 
19
c      subroutine qzvec(nm,n,a,b,epsb,alfr,alfi,beta,z)
 
20
c
 
21
c
 
22
c     on input:
 
23
c
 
24
c        nm must be set to the row dimension of two-dimensional
 
25
c          array parameters as declared in the calling program
 
26
c          dimension statement;
 
27
c
 
28
c        n is the order of the matrices;
 
29
c
 
30
c        a contains a real upper quasi-triangular matrix;
 
31
c
 
32
c        b contains a real upper triangular matrix.
 
33
c          computed and saved in  qzit;
 
34
c
 
35
c        alfr, alfi, and beta  are vectors with components whose
 
36
c          ratios ((alfr+i*alfi)/beta) are the generalized
 
37
c          eigenvalues.  they are usually obtained from  qzval;
 
38
c
 
39
c        z contains the transformation matrix produced in the
 
40
c          reductions by  qzhes,  qzit, and  qzval, if performed.
 
41
c          if the eigenvectors of the triangular problem are
 
42
c          desired, z must contain the identity matrix.
 
43
c
 
44
c     on output:
 
45
c
 
46
c        a is unaltered.  its subdiagonal elements provide information
 
47
c           about the storage of the complex eigenvectors;
 
48
c
 
49
c        b has been destroyed;
 
50
c
 
51
c        alfr, alfi, and beta are unaltered;
 
52
c
 
53
c        z contains the real and imaginary parts of the eigenvectors.
 
54
c          if alfi(i) .eq. 0.0, the i-th eigenvalue is real and
 
55
c            the i-th column of z contains its eigenvector.
 
56
c          if alfi(i) .ne. 0.0, the i-th eigenvalue is complex.
 
57
c            if alfi(i) .gt. 0.0, the eigenvalue is the first of
 
58
c              a complex pair and the i-th and (i+1)-th columns
 
59
c              of z contain its eigenvector.
 
60
c            if alfi(i) .lt. 0.0, the eigenvalue is the second of
 
61
c              a complex pair and the (i-1)-th and i-th columns
 
62
c              of z contain the conjugate of its eigenvector.
 
63
c          each eigenvector is normalized so that the modulus
 
64
c          of its largest component is 1.0 .
 
65
c
 
66
c! originator
 
67
c
 
68
c     this subroutine is the optional fourth step of the qz algorithm
 
69
c     for solving generalized matrix eigenvalue problems,
 
70
c     siam j. numer. anal. 10, 241-256(1973) by moler and stewart.
 
71
c    modification de la routine qzvec de eispack concernant le
 
72
c    passage de l'argument epsb.
 
73
c!
 
74
c
 
75
c     questions and comments should be directed to b. s. garbow,
 
76
c     applied mathematics division, argonne national laboratory
 
77
c
 
78
c     ------------------------------------------------------------------
 
79
c
 
80
      isw = 1
 
81
c     :::::::::: for en=n step -1 until 1 do -- ::::::::::
 
82
      do 800 nn = 1, n
 
83
         en = n + 1 - nn
 
84
         na = en - 1
 
85
         if (isw .eq. 2) go to 795
 
86
         if (alfi(en) .ne. 0.0d+0) go to 710
 
87
c     :::::::::: real vector ::::::::::
 
88
         m = en
 
89
         b(en,en) = 1.0d+0
 
90
         if (na .eq. 0) go to 800
 
91
         alfm = alfr(m)
 
92
         betm = beta(m)
 
93
c     :::::::::: for i=en-1 step -1 until 1 do -- ::::::::::
 
94
         do 700 ii = 1, na
 
95
            i = en - ii
 
96
            w = betm * a(i,i) - alfm * b(i,i)
 
97
            r = 0.0d+0
 
98
c
 
99
            do 610 j = m, en
 
100
  610       r = r + (betm * a(i,j) - alfm * b(i,j)) * b(j,en)
 
101
c
 
102
            if (i .eq. 1 .or. isw .eq. 2) go to 630
 
103
            if (betm * a(i,i-1) .eq. 0.0d+0) go to 630
 
104
            zz = w
 
105
            s = r
 
106
            go to 690
 
107
  630       m = i
 
108
            if (isw .eq. 2) go to 640
 
109
c     :::::::::: real 1-by-1 block ::::::::::
 
110
            t = w
 
111
            if (w .eq. 0.0d+0) t = epsb
 
112
            b(i,en) = -r / t
 
113
            go to 700
 
114
c     :::::::::: real 2-by-2 block ::::::::::
 
115
  640       x = betm * a(i,i+1) - alfm * b(i,i+1)
 
116
            y = betm * a(i+1,i)
 
117
            q = w * zz - x * y
 
118
            t = (x * s - zz * r) / q
 
119
            b(i,en) = t
 
120
            if (abs(x) .le. abs(zz)) go to 650
 
121
            b(i+1,en) = (-r - w * t) / x
 
122
            go to 690
 
123
  650       b(i+1,en) = (-s - y * t) / zz
 
124
  690       isw = 3 - isw
 
125
  700    continue
 
126
c     :::::::::: end real vector ::::::::::
 
127
         go to 800
 
128
c     :::::::::: complex vector ::::::::::
 
129
  710    m = na
 
130
         almr = alfr(m)
 
131
         almi = alfi(m)
 
132
         betm = beta(m)
 
133
c     :::::::::: last vector component chosen imaginary so that
 
134
c                eigenvector matrix is triangular ::::::::::
 
135
         y = betm * a(en,na)
 
136
         b(na,na) = -almi * b(en,en) / y
 
137
         b(na,en) = (almr * b(en,en) - betm * a(en,en)) / y
 
138
         b(en,na) = 0.0d+0
 
139
         b(en,en) = 1.0d+0
 
140
         enm2 = na - 1
 
141
         if (enm2 .eq. 0) go to 795
 
142
c     :::::::::: for i=en-2 step -1 until 1 do -- ::::::::::
 
143
         do 790 ii = 1, enm2
 
144
            i = na - ii
 
145
            w = betm * a(i,i) - almr * b(i,i)
 
146
            w1 = -almi * b(i,i)
 
147
            ra = 0.0d+0
 
148
            sa = 0.0d+0
 
149
c
 
150
            do 760 j = m, en
 
151
               x = betm * a(i,j) - almr * b(i,j)
 
152
               x1 = -almi * b(i,j)
 
153
               ra = ra + x * b(j,na) - x1 * b(j,en)
 
154
               sa = sa + x * b(j,en) + x1 * b(j,na)
 
155
  760       continue
 
156
c
 
157
            if (i .eq. 1 .or. isw .eq. 2) go to 770
 
158
            if (betm * a(i,i-1) .eq. 0.0d+0) go to 770
 
159
            zz = w
 
160
            z1 = w1
 
161
            r = ra
 
162
            s = sa
 
163
            isw = 2
 
164
            go to 790
 
165
  770       m = i
 
166
            if (isw .eq. 2) go to 780
 
167
c     :::::::::: complex 1-by-1 block ::::::::::
 
168
            tr = -ra
 
169
            ti = -sa
 
170
  773       dr = w
 
171
            di = w1
 
172
c     :::::::::: complex divide (t1,t2) = (tr,ti) / (dr,di) ::::::::::
 
173
  775       if (abs(di) .gt. abs(dr)) go to 777
 
174
            rr = di / dr
 
175
            d = dr + di * rr
 
176
            t1 = (tr + ti * rr) / d
 
177
            t2 = (ti - tr * rr) / d
 
178
            go to (787,782), isw
 
179
  777       rr = dr / di
 
180
            d = dr * rr + di
 
181
            t1 = (tr * rr + ti) / d
 
182
            t2 = (ti * rr - tr) / d
 
183
            go to (787,782), isw
 
184
c     :::::::::: complex 2-by-2 block ::::::::::
 
185
  780       x = betm * a(i,i+1) - almr * b(i,i+1)
 
186
            x1 = -almi * b(i,i+1)
 
187
            y = betm * a(i+1,i)
 
188
            tr = y * ra - w * r + w1 * s
 
189
            ti = y * sa - w * s - w1 * r
 
190
            dr = w * zz - w1 * z1 - x * y
 
191
            di = w * z1 + w1 * zz - x1 * y
 
192
            if (dr .eq. 0.0d+0 .and. di .eq. 0.0d+0) dr = epsb
 
193
            go to 775
 
194
  782       b(i+1,na) = t1
 
195
            b(i+1,en) = t2
 
196
            isw = 1
 
197
            if (abs(y) .gt. abs(w) + abs(w1)) go to 785
 
198
            tr = -ra - x * b(i+1,na) + x1 * b(i+1,en)
 
199
            ti = -sa - x * b(i+1,en) - x1 * b(i+1,na)
 
200
            go to 773
 
201
  785       t1 = (-r - zz * b(i+1,na) + z1 * b(i+1,en)) / y
 
202
            t2 = (-s - zz * b(i+1,en) - z1 * b(i+1,na)) / y
 
203
  787       b(i,na) = t1
 
204
            b(i,en) = t2
 
205
  790    continue
 
206
c     :::::::::: end complex vector ::::::::::
 
207
  795    isw = 3 - isw
 
208
  800 continue
 
209
c     :::::::::: end back substitution.
 
210
c                transform to original coordinate system.
 
211
c                for j=n step -1 until 1 do -- ::::::::::
 
212
      do 880 jj = 1, n
 
213
         j = n + 1 - jj
 
214
c
 
215
         do 880 i = 1, n
 
216
            zz = 0.0d+0
 
217
c
 
218
            do 860 k = 1, j
 
219
  860       zz = zz + z(i,k) * b(k,j)
 
220
c
 
221
            z(i,j) = zz
 
222
  880 continue
 
223
c     :::::::::: normalize so that modulus of largest
 
224
c                component of each vector is 1.
 
225
c                (isw is 1 initially from before) ::::::::::
 
226
      do 950 j = 1, n
 
227
         d = 0.0d+0
 
228
         if (isw .eq. 2) go to 920
 
229
         if (alfi(j) .ne. 0.0d+0) go to 945
 
230
c
 
231
         do 890 i = 1, n
 
232
            if (abs(z(i,j)) .gt. d) d = abs(z(i,j))
 
233
  890    continue
 
234
c
 
235
         do 900 i = 1, n
 
236
  900    z(i,j) = z(i,j) / d
 
237
c
 
238
         go to 950
 
239
c
 
240
  920    do 930 i = 1, n
 
241
            r = abs(z(i,j-1)) + abs(z(i,j))
 
242
            if (r .ne. 0.0d+0) r = r * sqrt((z(i,j-1)/r)**2
 
243
     x                                     +(z(i,j)/r)**2)
 
244
            if (r .gt. d) d = r
 
245
  930    continue
 
246
c
 
247
         do 940 i = 1, n
 
248
            z(i,j-1) = z(i,j-1) / d
 
249
            z(i,j) = z(i,j) / d
 
250
  940    continue
 
251
c
 
252
  945    isw = 3 - isw
 
253
  950 continue
 
254
c
 
255
      return
 
256
c     :::::::::: last card of qzvec ::::::::::
 
257
      end