Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pycrth.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pycrth.f
1 
2 C*********************************************************************
3 
4 C...PYCRTH
5 C...Auxiliary to PYEICG.
6 C
7 C THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF
8 C THE ALGOL PROCEDURE ORTHES, NUM. MATH. 12, 349-368(1968)
9 C BY MARTIN AND WILKINSON.
10 C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
11 C
12 C GIVEN A COMPLEX GENERAL MATRIX, THIS SUBROUTINE
13 C REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS
14 C LOW THROUGH IGH TO UPPER HESSENBERG FORM BY
15 C UNITARY SIMILARITY TRANSFORMATIONS.
16 C
17 C ON INPUT
18 C
19 C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
20 C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
21 C DIMENSION STATEMENT.
22 C
23 C N IS THE ORDER OF THE MATRIX.
24 C
25 C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
26 C SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED,
27 C SET LOW=1, IGH=N.
28 C
29 C AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
30 C RESPECTIVELY, OF THE COMPLEX INPUT MATRIX.
31 C
32 C ON OUTPUT
33 C
34 C AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
35 C RESPECTIVELY, OF THE HESSENBERG MATRIX. INFORMATION
36 C ABOUT THE UNITARY TRANSFORMATIONS USED IN THE REDUCTION
37 C IS STORED IN THE REMAINING TRIANGLES UNDER THE
38 C HESSENBERG MATRIX.
39 C
40 C ORTR AND ORTI CONTAIN FURTHER INFORMATION ABOUT THE
41 C TRANSFORMATIONS. ONLY ELEMENTS LOW THROUGH IGH ARE USED.
42 C
43 C CALLS PYTHAG FOR DSQRT(A*A + B*B) .
44 C
45 C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
46 C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
47 C
48 C THIS VERSION DATED AUGUST 1983.
49 C
50 
51  SUBROUTINE pycrth(NM,N,LOW,IGH,AR,AI,ORTR,ORTI)
52 
53  INTEGER i,j,m,n,ii,jj,la,mp,nm,igh,kp1,low
54  DOUBLE PRECISION ar(4,4),ai(4,4),ortr(4),orti(4)
55  DOUBLE PRECISION f,g,h,fi,fr,scale,pythag
56 
57  la = igh - 1
58  kp1 = low + 1
59  IF (la .LT. kp1) goto 210
60 C
61  DO 200 m = kp1, la
62  h = 0.0d0
63  ortr(m) = 0.0d0
64  orti(m) = 0.0d0
65  scale = 0.0d0
66 C .......... SCALE COLUMN (ALGOL TOL THEN NOT NEEDED) ..........
67  DO 100 i = m, igh
68  100 scale = scale + dabs(ar(i,m-1)) + dabs(ai(i,m-1))
69 C
70  IF (scale .EQ. 0.0d0) goto 200
71  mp = m + igh
72 C .......... FOR I=IGH STEP -1 UNTIL M DO -- ..........
73  DO 110 ii = m, igh
74  i = mp - ii
75  ortr(i) = ar(i,m-1) / scale
76  orti(i) = ai(i,m-1) / scale
77  h = h + ortr(i) * ortr(i) + orti(i) * orti(i)
78  110 CONTINUE
79 C
80  g = dsqrt(h)
81  f = pythag(ortr(m),orti(m))
82  IF (f .EQ. 0.0d0) goto 120
83  h = h + f * g
84  g = g / f
85  ortr(m) = (1.0d0 + g) * ortr(m)
86  orti(m) = (1.0d0 + g) * orti(m)
87  goto 130
88 C
89  120 ortr(m) = g
90  ar(m,m-1) = scale
91 C .......... FORM (I-(U*UT)/H) * A ..........
92  130 DO 160 j = m, n
93  fr = 0.0d0
94  fi = 0.0d0
95 C .......... FOR I=IGH STEP -1 UNTIL M DO -- ..........
96  DO 140 ii = m, igh
97  i = mp - ii
98  fr = fr + ortr(i) * ar(i,j) + orti(i) * ai(i,j)
99  fi = fi + ortr(i) * ai(i,j) - orti(i) * ar(i,j)
100  140 CONTINUE
101 C
102  fr = fr / h
103  fi = fi / h
104 C
105  DO 150 i = m, igh
106  ar(i,j) = ar(i,j) - fr * ortr(i) + fi * orti(i)
107  ai(i,j) = ai(i,j) - fr * orti(i) - fi * ortr(i)
108  150 CONTINUE
109 C
110  160 CONTINUE
111 C .......... FORM (I-(U*UT)/H)*A*(I-(U*UT)/H) ..........
112  DO 190 i = 1, igh
113  fr = 0.0d0
114  fi = 0.0d0
115 C .......... FOR J=IGH STEP -1 UNTIL M DO -- ..........
116  DO 170 jj = m, igh
117  j = mp - jj
118  fr = fr + ortr(j) * ar(i,j) - orti(j) * ai(i,j)
119  fi = fi + ortr(j) * ai(i,j) + orti(j) * ar(i,j)
120  170 CONTINUE
121 C
122  fr = fr / h
123  fi = fi / h
124 C
125  DO 180 j = m, igh
126  ar(i,j) = ar(i,j) - fr * ortr(j) - fi * orti(j)
127  ai(i,j) = ai(i,j) + fr * orti(j) - fi * ortr(j)
128  180 CONTINUE
129 C
130  190 CONTINUE
131 C
132  ortr(m) = scale * ortr(m)
133  orti(m) = scale * orti(m)
134  ar(m,m-1) = -g * ar(m,m-1)
135  ai(m,m-1) = -g * ai(m,m-1)
136  200 CONTINUE
137 C
138  210 RETURN
139  END