1      subroutine DEVVUN(A, LDA, N, EVALR, EVALI, VEC, IFLAG, WORK)
2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
3c ALL RIGHTS RESERVED.
4c Based on Government Sponsored Research NAS7-03001.
5c     This file contains DEVVUN and CDIV.
6c>> 1996-03-30 DEVUNN Krogh  MIN0 => MIN, added external statement.
7c>> 1994-11-02 DEVVUN Krogh  Changes to use M77CON
8c>> 1992-04-23 DEVVUN CLL   Declaring all variables.
9c>> 1992-04-22 DEVVUN Krogh Removed unused labels 330 and 1220.
10c>> 1992-03-05 DEVVUN Krogh Initial version, converted from EISPACK.
11c
12c     This subroutine uses slight modifcations of EISPACK routines
13c     BALANC, ELMHES, ELTRAN, HQR2 and BALBAK to get the eigenvalues and
14c     eigenvectors of a general real matrix by the QR method.  The first
15c     two are are encapsulated in DEVBH.  The following three are given
16c     inline here.
17c
18c     ELTRAN is a translation of the algol procedure ELMTRANS,
19c     Num. Math. 16, 181-204(1970) by Peters and Wilkinson.
20c     Handbook for Auto. Comp., Vol.ii-Linear Algebra, 372-395(1971).
21c
22c     HQR2 is a translation of the ALGOL procedure HQR2,
23c     Num. Math. 16, 181-204(1970) by Peters and Wilkinson.
24c     Handbook for Auto. Comp., Vol.ii-Linear Algebra, 372-395(1971).
25c
26c     BALBAK is a translation of the ALGOL procedure BALBAK,
27c     Num. Math. 13, 293-304(1969) by Parlett and Reinsch.
28c     Handbook for Auto. Comp., Vol.ii-Linear Algebra, 315-326(1971).
29c
30c     This subroutine finds the eigenvalues and eigenvectors
31c     of a real general matrix by the QR method.
32c
33c     On input
34c
35c     LDA  must be set to the row dimension of two-dimensional
36c          array parameters as declared in the calling program
37c          dimension statement.
38c     N    is the order of the matrix.
39c     A    contains the input matrix whose eigenvalues and eigenvectors
40c          are desired.
41c
42c     On output
43c     A    has been destroyed.
44c     EVALR and EVALI contain real and imaginary parts, respectively, of
45c          the eigenvalues.  The eigenvalues are given in order of
46c          increasing real parts.  When real parts are equal they are
47c          given in order of increasing absolute complex part.  Complex
48c          conjugate pairs of values appear consecutively with
49c          the eigenvalue having the positive imaginary part first.  If
50c          an error exit is made, the eigenvalues should be correct
51c          (but unordered) for indices IFLAG(1)+1,...,N.
52c     VEC  contains the real and imaginary parts of the eigenvectors.
53c          if the I-th eigenvalue is real, the I-th column of VEC
54c          contains its eigenvector.  If the I-th eigenvalue is complex
55c          with positive imaginary part, the I-th and (I+1)-th
56c          columns of VEC contain the real and imaginary parts of its
57c          eigenvector.  The eigenvectors are unnormalized.  If an
58c          error exit is made, none of the eigenvectors has been found.
59c     IFLAG(1) is set to
60c          1   If all eigenvalues are real.
61c          2   If some eigenvalues are complex.
62c          3   If N < 1 on the initial entry.
63c          4   If the limit of 30*N iterations is exhausted.
64c     ------------------------------------------------------------------
65c--D replaces "?": ?EVVUN, ?EVBH, ?NRM2, ?SCAL, ?CDIV
66c     ------------------------------------------------------------------
67      external DNRM2
68      integer I,J,K,L,M,N,EN,NA,LDA,IGH,ITN,ITS,LOW,MP,MP2,ENM2
69      integer II, IFLAG(N), LTYPE
70      double precision A(LDA,N),EVALR(N),EVALI(N),VEC(LDA,N), WORK(N)
71      double precision P,Q,R,S,T,W,X,Y,RA,SA,VI,VR,ZZ,NORM,TST1,TST2
72      double precision DNRM2
73      logical NOTLAS
74c     ------------------------------------------------------------------
75c
76      call DEVBH(A, LDA, N, LOW, IGH, IFLAG, WORK)
77      LTYPE = 1
78c
79c    -------------------------------- ELTRAN ---------------------------
80c
81c     Accumulate the stabilized elementary similarity transformations
82c     used in the preceding reduction of the matrix to upper Hessenberg
83c     form.
84c
85c     .......... initialize VEC to identity matrix ..........
86      do 20 J = 1, N
87         do 10 I = 1, N
88            VEC(I,J) = 0.0D0
89   10    continue
90         VEC(J,J) = 1.0D0
91   20 continue
92      if (IGH .gt. LOW+1) then
93         do 50 MP = IGH-1, LOW+1, -1
94            do 30 I = MP+1, IGH
95               VEC(I,MP) = A(I,MP-1)
96   30       continue
97            I = IFLAG(MP)
98            if (I .ne. MP) then
99               do 40 J = MP, IGH
100                  VEC(MP,J) = VEC(I,J)
101                  VEC(I,J) = 0.0D0
102   40          continue
103               VEC(I,MP) = 1.0D0
104            end if
105   50    continue
106      end if
107c
108c    -------------------------------- HQR2 -----------------------------
109c
110c     Find the eigenvalues of a real upper Hessenberg matrix by the QR
111c     method.
112c
113      IFLAG(1) = 0
114      NORM = 0.0D0
115      K = 1
116c     .......... store roots isolated by balanc
117c                and compute matrix norm ..........
118      do 70 I = 1, N
119         do 60 J = K, N
120            NORM = NORM + abs(A(I,J))
121   60    continue
122         K = I
123         if ((I .lt. LOW) .or. (I .gt. IGH)) then
124            EVALR(I) = A(I,I)
125            EVALI(I) = 0.0D0
126         end if
127   70 continue
128      EN = IGH
129      T = 0.0D0
130      ITN = 30*N
131c     .......... search for next eigenvalues ..........
132   80 if (EN .lt. LOW) go to 340
133      ITS = 0
134      NA = EN - 1
135      ENM2 = NA - 1
136c     .......... look for single small sub-diagonal element
137   90 do 100 L = EN, LOW+1, -1
138         S = abs(A(L-1,L-1)) + abs(A(L,L))
139         if (S .eq. 0.0D0) S = NORM
140         TST1 = S
141         TST2 = TST1 + abs(A(L,L-1))
142         if (TST2 .eq. TST1) go to 110
143  100 continue
144c     .......... form shift ..........
145  110 X = A(EN,EN)
146      if (L .eq. EN) go to 270
147      Y = A(NA,NA)
148      W = A(EN,NA) * A(NA,EN)
149      if (L .eq. NA) go to 280
150      if (ITN .le. 0) then
151c     .......... set error -- all eigenvalues have not
152c                converged after 30*N iterations ..........
153         call ERMSG('DEVVUN', EN, 0,
154     1'ERROR NO. is index of eigenvalue causing convergence failure.',
155     2            '.')
156         IFLAG(1) = 4
157         if (EN .le. 0) IFLAG(1) = 3
158         return
159      end if
160      if (ITS .ne. 10 .and. ITS .ne. 20) go to 130
161c     .......... form exceptional shift ..........
162      T = T + X
163      do 120 I = LOW, EN
164         A(I,I) = A(I,I) - X
165  120 continue
166      S = abs(A(EN,NA)) + abs(A(NA,ENM2))
167      X = 0.75D0 * S
168      Y = X
169      W = -0.4375D0 * S * S
170  130 ITS = ITS + 1
171      ITN = ITN - 1
172c     .......... look for two consecutive small sub-diagonal elements.
173      do 140 M = ENM2, L, -1
174         ZZ = A(M,M)
175         R = X - ZZ
176         S = Y - ZZ
177         P = (R * S - W) / A(M+1,M) + A(M,M+1)
178         Q = A(M+1,M+1) - ZZ - R - S
179         R = A(M+2,M+1)
180         S = abs(P) + abs(Q) + abs(R)
181         P = P / S
182         Q = Q / S
183         R = R / S
184         if (M .eq. L) go to 150
185         TST1 = abs(P)*(abs(A(M-1,M-1)) + abs(ZZ) + abs(A(M+1,M+1)))
186         TST2 = TST1 + abs(A(M,M-1))*(abs(Q) + abs(R))
187         if (TST2 .eq. TST1) go to 150
188  140 continue
189  150 MP2 = M + 2
190      do 160 I = MP2, EN
191         A(I,I-2) = 0.0D0
192         if (I .ne. MP2) A(I,I-3) = 0.0D0
193  160 continue
194c     .......... double QR step involving rows L to EN and
195c                columns M to EN ..........
196      do 260 K = M, NA
197         NOTLAS = K .ne. NA
198         if (K .ne. M) then
199            P = A(K,K-1)
200            Q = A(K+1,K-1)
201            R = 0.0D0
202            if (NOTLAS) R = A(K+2,K-1)
203            X = abs(P) + abs(Q) + abs(R)
204            if (X .eq. 0.0D0) go to 260
205            P = P / X
206            Q = Q / X
207            R = R / X
208         end if
209         S = sign(sqrt(P*P+Q*Q+R*R),P)
210         if (K .ne. M) then
211            A(K,K-1) = -S * X
212         else
213            if (L .ne. M) A(K,K-1) = -A(K,K-1)
214         end if
215         P = P + S
216         X = P / S
217         Y = Q / S
218         ZZ = R / S
219         Q = Q / P
220         R = R / P
221         if (NOTLAS) then
222c        .......... row modification ..........
223            do 200 J = K, N
224               P = A(K,J) + Q * A(K+1,J) + R * A(K+2,J)
225               A(K,J) = A(K,J) - P * X
226               A(K+1,J) = A(K+1,J) - P * Y
227               A(K+2,J) = A(K+2,J) - P * ZZ
228  200       continue
229            J = min(EN,K+3)
230c        .......... column modification ..........
231            do 210 I = 1, J
232               P = X * A(I,K) + Y * A(I,K+1) + ZZ * A(I,K+2)
233               A(I,K) = A(I,K) - P
234               A(I,K+1) = A(I,K+1) - P * Q
235               A(I,K+2) = A(I,K+2) - P * R
236  210       continue
237c        .......... accumulate transformations ..........
238            do 220 I = LOW, IGH
239               P = X * VEC(I,K) + Y * VEC(I,K+1) + ZZ * VEC(I,K+2)
240               VEC(I,K) = VEC(I,K) - P
241               VEC(I,K+1) = VEC(I,K+1) - P * Q
242               VEC(I,K+2) = VEC(I,K+2) - P * R
243  220       continue
244         else
245c        .......... row modification ..........
246            do 230 J = K, N
247               P = A(K,J) + Q * A(K+1,J)
248               A(K,J) = A(K,J) - P * X
249               A(K+1,J) = A(K+1,J) - P * Y
250  230       continue
251c
252            J = min(EN,K+3)
253c        .......... column modification ..........
254            do 240 I = 1, J
255               P = X * A(I,K) + Y * A(I,K+1)
256               A(I,K) = A(I,K) - P
257               A(I,K+1) = A(I,K+1) - P * Q
258  240       continue
259c        .......... accumulate transformations ..........
260            do 250 I = LOW, IGH
261               P = X * VEC(I,K) + Y * VEC(I,K+1)
262               VEC(I,K) = VEC(I,K) - P
263               VEC(I,K+1) = VEC(I,K+1) - P * Q
264  250       continue
265         end if
266  260 continue
267      go to 90
268c     .......... one root found ..........
269  270 A(EN,EN) = X + T
270      EVALR(EN) = A(EN,EN)
271      EVALI(EN) = 0.0D0
272      EN = NA
273      go to 80
274c     .......... two roots found ..........
275  280 P = (Y - X) / 2.0D0
276      Q = P * P + W
277      ZZ = sqrt(abs(Q))
278      A(EN,EN) = X + T
279      X = A(EN,EN)
280      A(NA,NA) = Y + T
281      if (Q .ge. 0.0D0) then
282c        .......... real pair ..........
283         ZZ = P + sign(ZZ,P)
284         EVALR(NA) = X + ZZ
285         EVALR(EN) = EVALR(NA)
286         if (ZZ .ne. 0.0D0) EVALR(EN) = X - W / ZZ
287         EVALI(NA) = 0.0D0
288         EVALI(EN) = 0.0D0
289         X = A(EN,NA)
290         S = abs(X) + abs(ZZ)
291         P = X / S
292         Q = ZZ / S
293         R = sqrt(P*P+Q*Q)
294         P = P / R
295         Q = Q / R
296c        .......... row modification ..........
297         do 290 J = NA, N
298            ZZ = A(NA,J)
299            A(NA,J) = Q * ZZ + P * A(EN,J)
300            A(EN,J) = Q * A(EN,J) - P * ZZ
301  290    continue
302c        .......... column modification ..........
303         do 300 I = 1, EN
304            ZZ = A(I,NA)
305            A(I,NA) = Q * ZZ + P * A(I,EN)
306            A(I,EN) = Q * A(I,EN) - P * ZZ
307  300    continue
308c        .......... accumulate transformations ..........
309         do 310 I = LOW, IGH
310            ZZ = VEC(I,NA)
311            VEC(I,NA) = Q * ZZ + P * VEC(I,EN)
312            VEC(I,EN) = Q * VEC(I,EN) - P * ZZ
313  310    continue
314      else
315c        .......... complex pair ..........
316         LTYPE = 2
317         EVALR(NA) = X + P
318         EVALR(EN) = X + P
319         EVALI(NA) = ZZ
320         EVALI(EN) = -ZZ
321      end if
322      EN = ENM2
323      go to 80
324c     .......... all roots found.  Backsubstitute to find
325c                vectors of upper triangular form ..........
326  340 if (NORM .eq. 0.0D0) go to 1000
327      do 800 EN = N, 1, -1
328         P = EVALR(EN)
329         Q = EVALI(EN)
330         NA = EN - 1
331         if (Q) 710, 600, 800
332c     .......... real vector ..........
333  600    M = EN
334         A(EN,EN) = 1.0D0
335         do 700 I = NA, 1, -1
336            W = A(I,I) - P
337            R = 0.0D0
338            do 610 J = M, EN
339               R = R + A(I,J) * A(J,EN)
340  610       continue
341            if (EVALI(I) .lt. 0.0D0) then
342               ZZ = W
343               S = R
344            else
345               M = I
346               if (EVALI(I) .eq. 0.0D0) then
347                  T = W
348                  if (T .eq. 0.0D0) then
349                     TST1 = NORM
350                     T = TST1
351  640                T = 0.01D0 * T
352                     TST2 = NORM + T
353                     if (TST2 .gt. TST1) go to 640
354                  end if
355                  A(I,EN) = -R / T
356               else
357c           .......... solve real equations ..........
358                  X = A(I,I+1)
359                  Y = A(I+1,I)
360                  Q = (EVALR(I)-P) * (EVALR(I)-P) + EVALI(I) * EVALI(I)
361                  T = (X * S - ZZ * R) / Q
362                  A(I,EN) = T
363                  if (abs(X) .gt. abs(ZZ)) then
364                     A(I+1,EN) = (-R - W * T) / X
365                  else
366                     A(I+1,EN) = (-S - Y * T) / ZZ
367                  end if
368               end if
369c           .......... overflow control ..........
370               T = abs(A(I,EN))
371               if (T .ne. 0.0D0) then
372                  TST1 = T
373                  TST2 = TST1 + 1.0D0/TST1
374                  if (TST2 .le. TST1) then
375                     do 690 J = I, EN
376                        A(J,EN) = A(J,EN)/T
377  690                continue
378                  end if
379               end if
380            end if
381  700    continue
382c     .......... end real vector ..........
383         go to 800
384c
385c     .......... complex vector ..........
386  710    M = NA
387c     .......... last vector component chosen imaginary so that
388c                eigenvector matrix is triangular ..........
389         if (abs(A(EN,NA)) .gt. abs(A(NA,EN))) then
390            A(NA,NA) = Q / A(EN,NA)
391            A(NA,EN) = -(A(EN,EN) - P) / A(EN,NA)
392         else
393            call DCDIV(0.0D0,-A(NA,EN),A(NA,NA)-P,Q,A(NA,NA),A(NA,EN))
394         end if
395         A(EN,NA) = 0.0D0
396         A(EN,EN) = 1.0D0
397         ENM2 = NA - 1
398         do 760 I = ENM2, 1, -1
399            W = A(I,I) - P
400            RA = 0.0D0
401            SA = 0.0D0
402            do 730 J = M, EN
403               RA = RA + A(I,J) * A(J,NA)
404               SA = SA + A(I,J) * A(J,EN)
405  730       continue
406            if (EVALI(I) .lt. 0.0D0) then
407               ZZ = W
408               R = RA
409               S = SA
410            else
411               M = I
412               if (EVALI(I) .eq. 0.0D0) then
413                  call DCDIV(-RA,-SA,W,Q,A(I,NA),A(I,EN))
414               else
415c           .......... solve complex equations ..........
416                  X = A(I,I+1)
417                  Y = A(I+1,I)
418                  VR = (EVALR(I)-P)*(EVALR(I)-P)+EVALI(I)*EVALI(I)-Q*Q
419                  VI = (EVALR(I) - P) * 2.0D0 * Q
420                  if (VR .eq. 0.0D0 .and. VI .eq. 0.0D0) then
421                     TST1 = NORM * (abs(W)+abs(Q)+abs(X)+abs(Y)+abs(ZZ))
422                     VR = TST1
423  740                VR = 0.01D0 * VR
424                     TST2 = TST1 + VR
425                     if (TST2 .gt. TST1) go to 740
426                  end if
427                  call DCDIV(X*R-ZZ*RA+Q*SA,X*S-ZZ*SA-Q*RA,VR,VI,
428     x                   A(I,NA),A(I,EN))
429                  if (abs(X) .gt. abs(ZZ) + abs(Q)) then
430                     A(I+1,NA) = (-RA - W * A(I,NA) + Q * A(I,EN)) / X
431                     A(I+1,EN) = (-SA - W * A(I,EN) - Q * A(I,NA)) / X
432                  else
433                     call DCDIV(-R-Y*A(I,NA),-S-Y*A(I,EN),ZZ,Q,
434     x                   A(I+1,NA),A(I+1,EN))
435                  end if
436               end if
437c           .......... overflow control ..........
438               T = max(abs(A(I,NA)), abs(A(I,EN)))
439               if (T .ne. 0.0D0) then
440                  TST1 = T
441                  TST2 = TST1 + 1.0D0/TST1
442                  if (TST2 .le. TST1) then
443                     do 750 J = I, EN
444                        A(J,NA) = A(J,NA)/T
445                        A(J,EN) = A(J,EN)/T
446  750                continue
447                  end if
448               end if
449            end if
450  760    continue
451c     .......... end complex vector ..........
452  800 continue
453c     .......... end back substitution.
454c                vectors of isolated roots ..........
455      do 840 I = 1, N
456         if (I .ge. LOW .and. I .le. IGH) go to 840
457c
458         do 820 J = I, N
459            VEC(I,J) = A(I,J)
460  820    continue
461c
462  840 continue
463c     .......... multiply by transformation matrix to give
464c                vectors of original full matrix.
465      do 880 J = N, LOW, -1
466         M = min(J,IGH)
467c
468         do 880 I = LOW, IGH
469            ZZ = 0.0D0
470c
471            do 860 K = LOW, M
472               ZZ = ZZ + VEC(I,K) * A(K,J)
473  860       continue
474c
475            VEC(I,J) = ZZ
476  880 continue
477c
478 1000 continue
479c
480c    ------------------------------ BALBAK -----------------------------
481c
482c     Form eigenvectors of a real general matrix by back transforming
483c     those of the corresponding balanced matrix determined by BALANC.
484c
485      if (IGH .ne. LOW) then
486         do 1110 I = LOW, IGH
487            S = WORK(I)
488c        .......... left hand eigenvectors are back transformed
489c                   if the foregoing statement is replaced by
490c                   S=1.0D0/WORK(I). ..........
491            do 1100 J = 1, N
492               VEC(I,J) = VEC(I,J) * S
493 1100       continue
494 1110    continue
495      end if
496c     ......... for I=LOW-1 step -1 until 1,
497c               IGH+1 step 1 until N do -- ..........
498      do 1140 II = 1, N
499         I = II
500         if ((I .lt. LOW) .or. (I .gt. IGH)) then
501            if (I .lt. LOW) I = LOW - II
502            K = WORK(I)
503            if (K .ne. I) then
504               do 1130 J = 1, N
505                  S = VEC(I,J)
506                  VEC(I,J) = VEC(K,J)
507                  VEC(K,J) = S
508 1130          continue
509            end if
510         end if
511 1140 continue
512c                        Normalize the eigenvectors
513      do 1220 I = 1, N
514         P = DNRM2(N, VEC(1, I), 1)
515         if (EVALI(I) .eq. 0.D0) then
516            call DSCAL(N, sign(1.D0, VEC(1,I)) / P, VEC(1,I), 1)
517         else if (EVALI(I) .gt. 0.D0) then
518            Q = DNRM2(N, VEC(1, I+1), 1)
519            if (P .gt. Q) then
520               P = P * sqrt(1.D0 + (Q/P)**2)
521            else if (Q .ne. 0.D0) then
522               P = Q * sqrt(1.D0 + (P/Q)**2)
523            else
524               go to 1220
525            end if
526            if (VEC(1, I+1) .eq. 0.D0) then
527               P = sign(1.D0, VEC(1,I+1)) / P
528               Q = 0.D0
529            else if (abs(VEC(1, I)) .gt. abs(VEC(1, I+1))) then
530               P = 1.D0 / (P*sqrt(1.D0+(VEC(1,I+1)/VEC(1,I))**2))
531               Q = P * (VEC(1,I+1) / abs(VEC(1,I)))
532               P = sign(P, VEC(1, I))
533            else
534               Q = 1.D0 / (P*sqrt(1.D0+(VEC(1,I)/VEC(1,I+1))**2))
535               P = Q * (VEC(1,I) / abs(VEC(1,I+1)))
536               Q = sign(Q, VEC(1, I+1))
537            end if
538            do 1200 J = 1, N
539               R = P * VEC(J, I) + Q * VEC(J, I+1)
540               VEC(J, I+1) = P * VEC(J, I+1) - Q * VEC(J, I)
541               VEC(J, I) = R
542 1200       continue
543         end if
544 1220 continue
545c-- Begin mask code changes
546c                        Set up for Shell sort
547c             Sort so real parts are algebraically increasing
548c             For = real parts, so abs(imag. parts) are increasing
549c             For both =, sort on index -- preserves complex pair order
550c-- End mask code changes
551      do 2000 I = 1, N
552         IFLAG(I) = I
553 2000 continue
554      L = 1
555      do 2010 K = 1, N
556         L = 3*L + 1
557         if (L .ge. N) go to 2020
558 2010 continue
559 2020 L = max(1, (L-4) / 9)
560 2030 do 2100 J = L+1, N
561         K = IFLAG(J)
562         P = EVALR(K)
563         I = J - L
564 2040    if (P - EVALR(IFLAG(I))) 2070, 2050, 2080
565 2050    if (abs(EVALI(K)) - abs(EVALI(IFLAG(I)))) 2070, 2060, 2080
566 2060    if (K .gt. IFLAG(I)) go to 2080
567 2070    IFLAG(I+L) = IFLAG(I)
568         I = I - L
569         if (I .gt. 0) go to 2040
570 2080    IFLAG(I+L) = K
571 2100 continue
572      L = (L-1) / 3
573      if (L .ne. 0) go to 2030
574c              Indices in IFLAG now give the desired order --
575c              Move entries to get this order.
576 2110 do 2150 I = L+1, N
577         if (IFLAG(I) .ne. I) then
578            L = I
579            M = I
580            P = EVALR(I)
581            Q = EVALI(I)
582            do 2115 J = 1, N
583               WORK(J) = VEC(J, I)
584 2115       continue
585 2120       K = IFLAG(M)
586            IFLAG(M) = M
587            if (K .ne. L) then
588               EVALR(M) = EVALR(K)
589               EVALI(M) = EVALI(K)
590               do 2130 J = 1, N
591                  VEC(J, M) = VEC(J, K)
592 2130          continue
593               M = K
594               go to 2120
595            else
596               EVALR(M) = P
597               EVALI(M) = Q
598               do 2140 J = 1, N
599                  VEC(J, M) = WORK(J)
600 2140          continue
601               go to 2110
602            end if
603         end if
604 2150 continue
605      IFLAG(1) = LTYPE
606      return
607      end
608c     ==================================================================
609      subroutine DCDIV(AR,AI,BR,BI,CR,CI)
610c
611c     complex division, (CR,CI) = (AR,AI)/(BR,BI)
612c     ------------------------------------------------------------------
613      double precision AR,AI,BR,BI,CR,CI
614      double precision S,ARS,AIS,BRS,BIS
615      S = abs(BR) + abs(BI)
616      ARS = AR/S
617      AIS = AI/S
618      BRS = BR/S
619      BIS = BI/S
620      S = BRS**2 + BIS**2
621      CR = (ARS*BRS + AIS*BIS)/S
622      CI = (AIS*BRS - ARS*BIS)/S
623      return
624      end
625