C Double precision version of CHEBY C PROGRAM CHEBYD IMPLICIT NONE INTEGER size PARAMETER(size=64) INTEGER m,i,j REAL*8 x(0:size),D1(0:size,0:size),D2(0:size,0:size) C eigen-deomposition variables INTEGER info REAL*8 vl,work(10*size) REAL*8 D2int(size,size),wr(size),wi(size),vr(size,size) WRITE(6,*) 'Enter m :' READ(5,*) m IF (size .LT. m) THEN WRITE(6,*) 'SIZE in CHEBY is too small!!' STOP END IF CALL CHEB(x,D1,D2,m,size) WRITE(6,*) 'Back in CHEBY' WRITE(6,*) 'x' DO i=0,m WRITE(6,*) x(i) END DO WRITE(6,*) 'D1' DO i=0,m WRITE(6,*) (D1(i,j),j=0,m) END DO WRITE(6,*) 'D2' DO i=0,m WRITE(6,*) (D2(i,j),j=0,m) END DO DO j=1,m-1 DO i=1,m-1 D2int(i,j)=D2(i,j) END DO END DO CALL DGEEV('N','V',m-1,D2int,size,wr,wi,vl,1,vr,size, & work,10*size,info) WRITE(6,*) ' ' WRITE(6,*) 'DGEEV info : ',info WRITE(6,*) 'evals' DO i=1,m-1 WRITE(6,*) wr(i),wi(i) END DO DO j=1,m-1 WRITE(6,*) 'evector #',j DO i=1,m-1 WRITE(6,*) vr(i,j) END DO END DO STOP END !!!!!!!!!!!!!!!!!!!!! SUBROUTINE CHEB(x,D1,D2,m,arrdim) IMPLICIT NONE C input/output variables INTEGER m,arrdim REAL*8 x(arrdim+1),D1(arrdim+1,arrdim+1),D2(arrdim+1,arrdim+1) C local parameters INTEGER size PARAMETER(size=65) REAL*8 pi PARAMETER(pi=3.14159265358979323846264338327950288D+0) C local variables INTEGER i,j,k,n,n1,n2 REAL*8 th(size),z(size,size),dx(size,size) REAL*8 cvec(size),c(size,size),sum C executable code IF (size .LT. m+1) THEN WRITE(6,*) 'In CHEB: SIZE is too small!!' WRITE(6,*) 'Must be at least ',m+1 STOP END IF n=m+1 j=m DO i=1,n th(i)=DBLE(i-1)*pi/DBLE(m) x(i)=DSIN(DBLE(j)*pi/DBLE(2*m)) j=j-2 END DO DO i=1,n DO j=1,n z(i,j)=2.0D+0*DSIN(th(j)/2.0D+0+th(i)/2.0D+0)* & DSIN(th(j)/2.0D+0-th(i)/2.0D+0) END DO END DO n1=n/2 n2=n-n1 DO j=1,n DO i=1,n1 dx(i,j)=z(i,j) END DO DO i=1,n2 dx(n1+i,j)=-z(n2+1-i,n+1-j) END DO dx(j,j)=1.0D+0 END DO DO j=1,n DO i=1,n z(i,j)=1.0D+0/dx(i,j) END DO z(j,j)=0.0D+0 END DO cvec(1)=1.0D+0 DO i=2,n cvec(i)=-cvec(i-1) END DO DO i=1,n k=1 DO j=i,n c(i,j)=cvec(k) k=k+1 END DO k=2 DO j=i-1,1,-1 c(i,j)=cvec(k) k=k+1 END DO END DO DO i=2,n-1 c(i,1)=c(i,1)/2.0D+0 c(i,n)=c(i,n)/2.0D+0 END DO DO j=2,n-1 c(1,j)=2.0D+0*c(1,j) c(n,j)=2.0D+0*c(n,j) END DO C compute D1 DO j=1,n DO i=1,n D1(i,j)=z(i,j)*c(i,j) END DO D1(j,j)=z(j,j)*(c(j,j)-1.0D+0) END DO DO i=1,n sum=0.0D+0 DO j=1,n sum=sum+D1(i,j) END DO D1(i,i)=-sum END DO C compute D2 DO j=1,n DO i=1,n D2(i,j)=2.0D+0*z(i,j)*(c(i,j)*D1(i,i)-D1(i,j)) END DO END DO DO i=1,n sum=0.0D+0 DO j=1,n sum=sum+D2(i,j) END DO D2(i,i)=-sum END DO RETURN END