← Back to branch summary

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

« back to all changes in this revision

Viewing changes to routines/control/corth.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/corth.f
 
1
      subroutine corth(nm,n,low,igh,ar,ai,ortr,orti)
 
2
c
 
3
      integer i,j,m,n,ii,jj,la,mp,nm,igh,kp1,low
 
4
      double precision ar(nm,n),ai(nm,n),ortr(igh),orti(igh)
 
5
      double precision f,g,h,fi,fr,scale
 
6
c
 
7
c!purpose
 
8
c
 
9
c     given a complex general matrix, this subroutine
 
10
c     reduces a submatrix situated in rows and columns
 
11
c     low through igh to upper hessenberg form by
 
12
c     unitary similarity transformations.
 
13
c
 
14
c!calling sequence
 
15
c     subroutine corth(nm,n,low,igh,ar,ai,ortr,orti)
 
16
c
 
17
c     integer i,j,m,n,ii,jj,la,mp,nm,igh,kp1,low
 
18
c     double precision ar(nm,n),ai(nm,n),ortr(igh),orti(igh)
 
19
c     double precision f,g,h,fi,fr,scale
 
20
c
 
21
c     on input:
 
22
c
 
23
c        nm must be set to the row dimension of two-dimensional
 
24
c          array parameters as declared in the calling program
 
25
c          dimension statement;
 
26
c
 
27
c        n is the order of the matrix;
 
28
c
 
29
c        low and igh are integers determined by the balancing
 
30
c          subroutine  cbal.  if  cbal  has not been used,
 
31
c          set low=1, igh=n;
 
32
c
 
33
c        ar and ai contain the real and imaginary parts,
 
34
c          respectively, of the complex input matrix.
 
35
c
 
36
c     on output:
 
37
c
 
38
c        ar and ai contain the real and imaginary parts,
 
39
c          respectively, of the hessenberg matrix.  information
 
40
c          about the unitary transformations used in the reduction
 
41
c          is stored in the remaining triangles under the
 
42
c          hessenberg matrix;
 
43
c
 
44
c        ortr and orti contain further information about the
 
45
c          transformations.  only elements low through igh are used.
 
46
c
 
47
c!originator
 
48
c
 
49
c     this subroutine is a translation of a complex analogue of
 
50
c     the algol procedure orthes, num. math. 12, 349-368(1968)
 
51
c     by martin and wilkinson.
 
52
c     handbook for auto. comp., vol.ii-linear algebra, 339-358(1971).
 
53
c     questions and comments should be directed to b. s. garbow,
 
54
c     applied mathematics division, argonne national laboratory
 
55
c
 
56
c!
 
57
c     ------------------------------------------------------------------
 
58
c
 
59
      la = igh - 1
 
60
      kp1 = low + 1
 
61
      if (la .lt. kp1) go to 200
 
62
c
 
63
      do 180 m = kp1, la
 
64
         h = 0.0d+0
 
65
         ortr(m) = 0.0d+0
 
66
         orti(m) = 0.0d+0
 
67
         scale = 0.0d+0
 
68
c     :::::::::: scale column (algol tol then not needed) ::::::::::
 
69
         do 90 i = m, igh
 
70
   90    scale = scale + abs(ar(i,m-1)) + abs(ai(i,m-1))
 
71
c
 
72
         if (scale .eq. 0.0d+0) go to 180
 
73
         mp = m + igh
 
74
c     :::::::::: for i=igh step -1 until m do -- ::::::::::
 
75
         do 100 ii = m, igh
 
76
            i = mp - ii
 
77
            ortr(i) = ar(i,m-1) / scale
 
78
            orti(i) = ai(i,m-1) / scale
 
79
            h = h + ortr(i) * ortr(i) + orti(i) * orti(i)
 
80
  100    continue
 
81
c
 
82
         g = sqrt(h)
 
83
        f = sqrt(ortr(m)*ortr(m)+orti(m)*orti(m))
 
84
         if (f .eq. 0.0d+0) go to 103
 
85
         h = h + f * g
 
86
         g = g / f
 
87
         ortr(m) = (1.0d+0 + g) * ortr(m)
 
88
         orti(m) = (1.0d+0 + g) * orti(m)
 
89
         go to 105
 
90
c
 
91
  103    ortr(m) = g
 
92
         ar(m,m-1) = scale
 
93
c     :::::::::: form (i-(u*ut)/h) * a ::::::::::
 
94
  105    do 130 j = m, n
 
95
            fr = 0.0d+0
 
96
            fi = 0.0d+0
 
97
c     :::::::::: for i=igh step -1 until m do -- ::::::::::
 
98
            do 110 ii = m, igh
 
99
               i = mp - ii
 
100
               fr = fr + ortr(i) * ar(i,j) + orti(i) * ai(i,j)
 
101
               fi = fi + ortr(i) * ai(i,j) - orti(i) * ar(i,j)
 
102
  110       continue
 
103
c
 
104
            fr = fr / h
 
105
            fi = fi / h
 
106
c
 
107
            do 120 i = m, igh
 
108
               ar(i,j) = ar(i,j) - fr * ortr(i) + fi * orti(i)
 
109
               ai(i,j) = ai(i,j) - fr * orti(i) - fi * ortr(i)
 
110
  120       continue
 
111
c
 
112
  130    continue
 
113
c     :::::::::: form (i-(u*ut)/h)*a*(i-(u*ut)/h) ::::::::::
 
114
         do 160 i = 1, igh
 
115
            fr = 0.0d+0
 
116
            fi = 0.0d+0
 
117
c     :::::::::: for j=igh step -1 until m do -- ::::::::::
 
118
            do 140 jj = m, igh
 
119
               j = mp - jj
 
120
               fr = fr + ortr(j) * ar(i,j) - orti(j) * ai(i,j)
 
121
               fi = fi + ortr(j) * ai(i,j) + orti(j) * ar(i,j)
 
122
  140       continue
 
123
c
 
124
            fr = fr / h
 
125
            fi = fi / h
 
126
c
 
127
            do 150 j = m, igh
 
128
               ar(i,j) = ar(i,j) - fr * ortr(j) - fi * orti(j)
 
129
               ai(i,j) = ai(i,j) + fr * orti(j) - fi * ortr(j)
 
130
  150       continue
 
131
c
 
132
  160    continue
 
133
c
 
134
         ortr(m) = scale * ortr(m)
 
135
         orti(m) = scale * orti(m)
 
136
         ar(m,m-1) = -g * ar(m,m-1)
 
137
         ai(m,m-1) = -g * ai(m,m-1)
 
138
  180 continue
 
139
c
 
140
  200 return
 
141
c     :::::::::: last card of corth ::::::::::
 
142
      end