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