1C> \ingroup selci
2C> @{
3C>
4      subroutine selci_yacobi(a,u,n,big,jb)
5      implicit real*8(a-h,o-z)
6*
7* $Id$
8*
9      dimension a(1),u(1),big(1),jb(1)
10      integer*4 isum1
11      data zero/0.0d0/
12c 32 bit precision
13c      data eps/1.0e-7/
14c 64 bit precision
15      data eps/1.0d-12/
16      data one/1.0d0/
17      l(i)=(i*(i-1))/2
18c
19      u(1)=one
20      if (n.eq.1) go to 999
21      j=0
22      k=0
23      n2=n*n
24      do 10 i=1,n2
25   10 u(i)=zero
26      do 30 i=1,n
27      j=j+i
28      jj=k+i
29      k=k+n
30      u(jj)=one
31c     section to locate biggest off-diagonal element in row i
32      im1=i-1
33      if (im1.eq.0) go to 30
34      big(i)=zero
35      do 20 jj=1,im1
36      ij=l(i)+jj
37      b=abs(a(ij))
38      if (b.le.big(i)) go to 20
39      big(i)=b
40      jb(i)=jj
41   20 continue
42   30 continue
43c     locate rotation pivot
44   31 continue
45      ibig=2
46      do 35 i=2,n
47      if (big(i).gt.big(ibig)) ibig=i
48   35 continue
49      biggst=big(ibig)
50      jbig=jb(ibig)
51      if (biggst.le.eps) go to 999
52c     begin rotation
53      ii=l(ibig)+ibig
54      ij=l(ibig)+jbig
55      jj=l(jbig)+jbig
56      ab=a(ii)-a(jj)
57      aa=a(ij)+a(ij)
58      d=sqrt(ab*ab+aa*aa)
59      if (ab.lt.zero) d=-d
60      t=aa/(ab+d)
61      tsq=t*t
62      csq=one/(one+tsq)
63      c=sqrt(csq)
64      s=c*t
65      ab=aa*t
66      d=(a(ii)+ab+tsq*a(jj))*csq
67      a(jj)=(a(jj)-ab+tsq*a(ii))*csq
68      a(ii)=d
69      a(ij)=zero
70      iii=4
71      kix=(ibig-1)*n
72      kjx=(jbig-1)*n
73      do 60 i=1,n
74      ki=kix+i
75      kj=kjx+i
76      d=c*u(ki)+s*u(kj)
77      u(kj)=c*u(kj)-s*u(ki)
78      u(ki)=d
79      isum1=i-jbig
80      if (isum1) 45,54,46
81   45 kj=l(jbig)+i
82      ki=l(ibig)+i
83      go to 50
84   46 kj=l(i)+jbig
85      isum1=i-ibig
86      if (isum1)47,55,48
87   47 ki=l(ibig)+i
88      iii=2
89      if (jb(i).eq.jbig) iii=3
90      go to 50
91   48 ki=l(i)+ibig
92      iii=1
93      if (jb(i).eq.ibig.or.jb(i).eq.jbig) iii=3
94   50 d=c*a(ki)+s*a(kj)
95      a(kj)=c*a(kj)-s*a(ki)
96      a(ki)=d
97      go to (51,52,55,60),iii
98   51 b=abs(a(ki))
99      if (b.le.big(i)) go to 52
100      big(i)=b
101      jb(i)=ibig
102   52 b=abs(a(kj))
103      if (b.le.big(i)) go to 60
104      big(i)=b
105      jb(i)=jbig
106      go to 60
107   54 if (i.eq.1) go to 60
108   55 im1=i-1
109      big(i)=zero
110      do 58 j=1,im1
111      ij=l(i)+j
112      b=abs(a(ij))
113      if (b.le.big(i)) go to 58
114      big(i)=b
115      jb(i)=j
116   58 continue
117   60 continue
118c     end of rotation
119      go to 31
120  999 continue
121c lock eigenvectors so that first coefficient is always
122c positive.
123      ij=1
124      do 501 i=1,n
125      if (u(ij).ge.zero) go to 503
126      do 502 j=1,n
127      u(ij+j-1)=-u(ij+j-1)
128  502 continue
129  503 ij=ij+n
130  501 continue
131c put eigen values in big and order them
132      ii = 0
133      do 1231 i = 1,n
134          ii = ii+i
135          big(i) = a(ii)
1361231  continue
137      call selci_order(n,big,u)
138      return
139      end
140      subroutine selci_order(n,e,v)
141      implicit real*8(a-h,o-z)
142      dimension e(n),v(n,n)
143c
144c order eigenvalues and vectors in increasing order.
145c
146      do 10 j=1,n
147          do 20 i=j+1,n
148              if(e(i).lt.e(j)) then
149                t=e(i)
150                e(i)=e(j)
151                e(j)=t
152                do 30 k=1,n
153                    t=v(k,i)
154                    v(k,i)=v(k,j)
15530                  v(k,j)=t
156              endif
15720        continue
15810    continue
159      return
160      end
161C>
162C> @}
163