SUBROUTINE ZDROTU(N, X, INCX, Y, INCY, F, S, T) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX*16 X(*),Y(*),V,T DOUBLE PRECISION F,S,UR,UC,VR,VC IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN DO I=1,N UR=DBLE(X(I)) UC=DIMAG(X(I)) V=Y(I)*T VR=DBLE(V) VC=DIMAG(V) X(I)=DCMPLX(S*(-F*UR+VR)+UR,S*(-F*UC+VC)+UC) Y(I)=DCMPLX(-S*(F*VR+UR)+VR,-S*(F*VC+UC)+VC) ENDDO ELSE DO I=1,N UR=DBLE(X((I-1)*INCX+1)) UC=DIMAG(X((I-1)*INCX+1)) V=Y((I-1)*INCY+1)*T VR=DBLE(V) VC=DIMAG(V) X((I-1)*INCX+1)=DCMPLX(S*(-F*UR+VR)+UR,S*(-F*UC+VC)+UC) Y((I-1)*INCY+1)=DCMPLX(-S*(F*VR+UR)+VR,-S*(F*VC+UC)+VC) ENDDO ENDIF RETURN END SUBROUTINE ZDROTU SUBROUTINE ZDROTU2(N, X, INCX, Y, INCY, F, S, T) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX*16 X(*),Y(*),V,T DOUBLE PRECISION F,S,UR,UC,VR,VC IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN DO I=1,N UR=DBLE(X(I)) UC=-DIMAG(X(I)) V=Y(I)*T VR=DBLE(V) VC=DIMAG(V) X(I)=DCMPLX(S*(-F*UR+VR)+UR,S*(F*UC-VC)-UC) Y(I)=DCMPLX(-S*(F*VR+UR)+VR,-S*(F*VC+UC)+VC) ENDDO ELSE DO I=1,N UR=DBLE(X((I-1)*INCX+1)) UC=-DIMAG(X((I-1)*INCX+1)) V=Y((I-1)*INCY+1)*T VR=DBLE(V) VC=DIMAG(V) X((I-1)*INCX+1)=DCMPLX(S*(-F*UR+VR)+UR,S*(F*UC-VC)-UC) Y((I-1)*INCY+1)=DCMPLX(-S*(F*VR+UR)+VR,-S*(F*VC+UC)+VC) ENDDO ENDIF RETURN END SUBROUTINE ZDROTU2 SUBROUTINE ZDROTV(N, X, INCX, Y, INCY, C, G, T) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX*16 X(*),Y(*),T,W DOUBLE PRECISION G,C,VR,VC,UR,UC IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN DO I=1,N UR=DBLE(X(I)) UC=DIMAG(X(I)) W=Y(I)*T VR=DBLE(W) VC=DIMAG(W) X(I)=DCMPLX(C*(-G*VR+UR)+VR,C*(-G*VC+UC)+VC) Y(I)=DCMPLX(C*(G*UR+VR)-UR,C*(G*UC+VC)-UC) ENDDO ELSE DO I=1,N UR=DBLE(X((I-1)*INCX+1)) UC=DIMAG(X((I-1)*INCX+1)) W=Y((I-1)*INCY+1)*T VR=DBLE(W) VC=DIMAG(W) X((I-1)*INCX+1)=DCMPLX(C*(-G*VR+UR)+VR,C*(-G*VC+UC)+VC) Y((I-1)*INCY+1)=DCMPLX(C*(G*UR+VR)-UR,C*(G*UC+VC)-UC) ENDDO ENDIF RETURN END SUBROUTINE ZDROTV SUBROUTINE ZDROTV2(N, X, INCX, Y, INCY, C, G, T) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX*16 X(*),Y(*),T,W DOUBLE PRECISION G,C,VR,VC,UR,UC IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN DO I=1,N UR=DBLE(X(I)) UC=-DIMAG(X(I)) W=Y(I)*T VR=DBLE(W) VC=DIMAG(W) X(I)=DCMPLX(C*(-G*VR+UR)+VR,C*(G*VC-UC)-VC) Y(I)=DCMPLX(C*(G*UR+VR)-UR,C*(G*UC+VC)-UC) ENDDO ELSE DO I=1,N UR=DBLE(X((I-1)*INCX+1)) UC=-DIMAG(X((I-1)*INCX+1)) W=Y((I-1)*INCY+1)*T VR=DBLE(W) VC=DIMAG(W) X((I-1)*INCX+1)=DCMPLX(C*(-G*VR+UR)+VR,C*(G*VC-UC)-VC) Y((I-1)*INCY+1)=DCMPLX(C*(G*UR+VR)-UR,C*(G*UC+VC)-UC) ENDDO ENDIF RETURN END SUBROUTINE ZDROTV2 SUBROUTINE ZDROTW(N, X, INCX, Y, INCY, C, G, T) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX*16 X(*),Y(*),T,W DOUBLE PRECISION G,C,VR,VC,UR,UC IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN DO I=1,N UR=DBLE(X(I)) UC=DIMAG(X(I)) W=Y(I)*T VR=DBLE(W) VC=DIMAG(W) X(I)=DCMPLX(C*(G*VR+UR)-VR,C*(G*VC+UC)-VC) Y(I)=DCMPLX(C*(-G*UR+VR)+UR,C*(-G*UC+VC)+UC) ENDDO ELSE DO I=1,N UR=DBLE(X((I-1)*INCX+1)) UC=DIMAG(X((I-1)*INCX+1)) W=Y((I-1)*INCY+1)*T VR=DBLE(W) VC=DIMAG(W) X((I-1)*INCX+1)=DCMPLX(C*(G*VR+UR)-VR,C*(G*VC+UC)-VC) Y((I-1)*INCY+1)=DCMPLX(C*(-G*UR+VR)+UR,C*(-G*UC+VC)+UC) ENDDO ENDIF RETURN END SUBROUTINE ZDROTW SUBROUTINE ZDROTW2(N, X, INCX, Y, INCY, C, G, T) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX*16 X(*),Y(*),T,W DOUBLE PRECISION G,C,VR,VC,UR,UC IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN DO I=1,N UR=DBLE(X(I)) UC=-DIMAG(X(I)) W=Y(I)*T VR=DBLE(W) VC=DIMAG(W) X(I)=DCMPLX(C*(G*VR+UR)-VR,-C*(G*VC+UC)+VC) Y(I)=DCMPLX(C*(-G*UR+VR)+UR,C*(-G*UC+VC)+UC) ENDDO ELSE DO I=1,N UR=DBLE(X((I-1)*INCX+1)) UC=-DIMAG(X((I-1)*INCX+1)) W=Y((I-1)*INCY+1)*T VR=DBLE(W) VC=DIMAG(W) X((I-1)*INCX+1)=DCMPLX(C*(G*VR+UR)-VR,-C*(G*VC+UC)+VC) Y((I-1)*INCY+1)=DCMPLX(C*(-G*UR+VR)+UR,C*(-G*UC+VC)+UC) ENDDO ENDIF RETURN END SUBROUTINE ZDROTW2 SUBROUTINE DLARTGU (F1, G1, C, S, R1) IMPLICIT NONE DOUBLE PRECISION ZERO, ONE, HALF PARAMETER (ONE = 1.0D0, ZERO = 0.0D0, HALF = 0.5D0) DOUBLE PRECISION F1, G1, C, S, R1, U, V EXTERNAL DFMA0 DOUBLE PRECISION DFMA0 IF (F1 .GE. G1) THEN U = G1/F1 R1 = SQRT(DFMA0(U,U,ONE)) C = ONE/R1 S = U/R1 R1 = R1*F1 ELSE V = F1/G1 R1 = SQRT(DFMA0(V,V,ONE)) S = ONE/R1 C = V/R1 R1 = R1*G1 ENDIF RETURN END SUBROUTINE DLARTGU SUBROUTINE ZJAEVD(N,A,LAMBDA,W,INFO) use omp_lib IMPLICIT NONE INTEGER N, INFO, FLAG, LL, MM DOUBLE PRECISION LAMBDA(N),TMP,B(N),Z(N),TMP3,TMP4 COMPLEX*16 A(N,N),W(N,N) COMPLEX*16 TMP2, TMP5 COMPLEX*16 ZZERO, ZONE PARAMETER (ZONE = (1.0D0,0.0D0), ZZERO = (0.0D0,0.0D0)) DOUBLE PRECISION DZERO, DONE, DHALF PARAMETER (DONE = 1.0D0, DZERO = 0.0D0, DHALF = 0.5D0) INTEGER I, J, K, ITER, MAXITER PARAMETER (MAXITER = 100) DOUBLE PRECISION T, U, C, S, F, EPS DOUBLE PRECISION THETA1, THETA2 DOUBLE PRECISION G1, H1, R1 EXTERNAL DFMA0, DLAMCH DOUBLE PRECISION DFMA0, DLAMCH EPS = DLAMCH('P') DO I = 1,N LAMBDA(I) = DBLE(A(I,I)) ENDDO LL = 1 MM = N DO ITER=1, MAXITER WRITE (*,*) 'ITER',ITER,LL,MM DO I=LL,MM-1 DO J=I+1,MM IF (ABS ( A(I,J) ) .GT. EPS*SQRT(ABS(LAMBDA(J)))*SQRT(ABS(LAMBDA(I)))) GO TO 10 ENDDO ENDDO 10 LL = I IF (LL .GE. MM) EXIT DO I=MM,LL+1,-1 DO J=I-1,LL,-1 IF (ABS ( A(J,I) ) .GT. EPS*SQRT(ABS(LAMBDA(J)))*SQRT(ABS(LAMBDA(I)))) GO TO 20 ENDDO ENDDO 20 MM = I IF (LL .GE. MM) EXIT FLAG = 0 DO I = LL,MM B(I) = LAMBDA(I) ENDDO DO I = LL,MM Z(I) = DZERO ENDDO DO J = LL, MM-1 K = J T = ABS(LAMBDA(J)) DO I = J+1, MM IF (ABS(LAMBDA(I)) .GT. T) THEN K = I T = ABS(LAMBDA(I)) ENDIF ENDDO IF (J .NE. K) THEN A(J,K) = DCONJG(A(J,K)) CALL DSWAP(1,LAMBDA(J),1,LAMBDA(K),1) CALL DSWAP(1,B(J),1,B(K),1) CALL DSWAP(1,Z(J),1,Z(K),1) CALL ZSWAP(J - LL, A(LL,J), 1, A(LL,K), 1) DO I = J+1, K-1 TMP2 = DCONJG(A(J,I)) A(J,I) = DCONJG(A(I,K)) A(I,K) = TMP2 ENDDO CALL ZSWAP(MM - K, A(J,K+1), N, A(K,K+1), N) CALL ZSWAP(N, W(1,J), 1, W(1,K), 1) ENDIF DO K = J+1, MM IF (ABS(A(J,K)) .LE. EPS*SQRT(ABS(LAMBDA(J)))*SQRT(ABS(LAMBDA(K)))) THEN A(J,K) = ZZERO ELSE TMP3 = DBLE(A(J,K)) TMP4 = DIMAG(A(J,K)) CALL DLARTG(ABS(TMP3),ABS(TMP4),C,S,TMP) C = SIGN(C,TMP3) S = SIGN(S,TMP4) TMP2 = DCMPLX(C,S) TMP5 = DCONJG(TMP2) H1 = DHALF*LAMBDA(J) - DHALF*LAMBDA(K) G1 = TMP CALL DLARTGU (ABS(H1), G1, C, S, R1) T = G1 / (H1 + SIGN(R1,H1)) R1 = SQRT(DFMA0(T,T,DONE)) C = DONE/R1 S = T/R1 Z(J) = DFMA0 (T, TMP, Z(J)) Z(K) = DFMA0 (-T, TMP, Z(K)) LAMBDA(J) = B(J)+Z(J) LAMBDA(K) = B(K)+Z(K) IF (ABS(LAMBDA(J)) .LT. ABS(LAMBDA(K))) THEN THETA1 = C THETA2 = S IF (THETA2 .GT. DZERO) THEN C = THETA2 S = -THETA1 ELSE C = -THETA2 S = THETA1 ENDIF TMP = LAMBDA(J) LAMBDA(J) = LAMBDA(K) LAMBDA(K) = TMP TMP = B(J) B(J) = B(K) B(K) = TMP TMP = Z(J) Z(J) = Z(K) Z(K) = TMP ENDIF IF (C .GE. ABS(S)) THEN F = S/(DONE+C) C = DFMA0(-S,F,DONE) IF (C .NE. DONE .OR. S .NE. DZERO) THEN FLAG = 1 CALL ZDROTU(J - LL, A(LL,J), 1, A(LL,K), 1, F, S, TMP5) CALL ZDROTU2(K - J - 1, A(J,J+1), N, A(J+1,K), 1, F, S, TMP5) CALL ZDROTU(MM - K, A(J,K+1), N, A(K,K+1), N, F, S, TMP2) CALL ZDROTU(N, W(1,J), 1, W(1,K), 1, F, S, TMP5) ELSE IF(TMP2 .NE. ZONE .OR. TMP5 .NE. ZONE) THEN FLAG = 1 CALL ZSCAL(MM - K, TMP2, A(K, K+1), N) CALL ZSCAL(K - LL, TMP5, A(LL, K), 1) CALL ZSCAL(N, TMP5, W(1,K), 1) ENDIF ELSE IF(S .GE. DZERO) THEN F = C/(DONE+S) S = DFMA0(-C,F,DONE) FLAG = 1 CALL ZDROTV(J - LL, A(LL,J), 1, A(LL,K), 1, C, F, TMP5) CALL ZDROTV2(K - J - 1, A(J,J+1), N, A(J+1,K), 1, C, F, TMP5) CALL ZDROTV(MM - K, A(J,K+1), N, A(K,K+1), N, C, F, TMP2) CALL ZDROTV(N, W(1,J), 1, W(1,K), 1, C, F, TMP5) ELSE IF(S .LE. DZERO) THEN F = C/(DONE-S) S = DFMA0(C,F,-DONE) FLAG = 1 CALL ZDROTW(J - LL, A(LL,J), 1, A(LL,K), 1, C, F, TMP5) CALL ZDROTW2(K - J - 1, A(J,J+1), N, A(J+1,K), 1, C, F, TMP5) CALL ZDROTW(MM - K, A(J,K+1), N, A(K,K+1), N, C, F, TMP2) CALL ZDROTW(N, W(1,J), 1, W(1,K), 1, C, F, TMP5) ENDIF A(J,K) = ZZERO ENDIF ENDDO ENDDO IF (FLAG .EQ. 0) EXIT ENDDO IF (ITER .GT. MAXITER) THEN INFO = -1 RETURN ENDIF DO I = 1, N-1 K = I T = LAMBDA(K) DO J = I+1, N IF (LAMBDA(J) .GT. T) THEN K = J T = LAMBDA(K) ENDIF ENDDO LAMBDA(K) = LAMBDA(I) LAMBDA(I) = T CALL ZSWAP(N, W(1,K), 1, W(1,I), 1) ENDDO INFO = 0 RETURN END SUBROUTINE ZJAEVD PROGRAM TEDTDJE use omp_lib IMPLICIT NONE INTEGER I, J, K, N, INFO, SEED DOUBLE PRECISION TMP1, TMP2 COMPLEX*16 E, S COMPLEX*16 ZZERO, ZONE PARAMETER (ZZERO = (0.0D0,0.0D0), ZONE= (1.0D0,0.0D0)) DOUBLE PRECISION DZERO, DONE, DTWO PARAMETER (DZERO = 0.0D0, DONE= 1.0D0, DTWO= 2.0D0) DOUBLE PRECISION,ALLOCATABLE :: LAMBDA(:) COMPLEX*16,ALLOCATABLE :: A(:,:), B(:,:), W(:,:) DOUBLE PRECISION USERRAND EXTERNAL USERRAND EXTERNAL ZDOTC COMPLEX*16 ZDOTC character*1000 ARG1,ARG2 integer t1, t2, t_rate, t_max, diff CALL GETARG(1,ARG1) CALL GETARG(2,ARG2) read (ARG1,*) N read (ARG2,*) SEED ALLOCATE(LAMBDA(N)) ALLOCATE(A(N,N), B(N,N), W(N,N)) DO I=1,SEED CALL XOR() ENDDO DO I=1,N DO J=1,I-1 TMP1 = USERRAND() TMP2 = USERRAND() A(J,I) = dcmplx(TMP1,TMP2) B(J,I) = A(J,I) B(I,J) = DCONJG(B(J,I)) ENDDO TMP1 = USERRAND() A(I,I) = dcmplx(TMP1,DZERO) B(I,I) = A(I,I) ENDDO DO J=1,N DO K=1,N W(K,J) = ZZERO ENDDO W(J,J) = ZONE ENDDO call system_clock(t1) ! 開始時を記録 CALL ZJAEVD(N,A,LAMBDA,W,INFO) call system_clock(t2, t_rate, t_max) ! 終了時を記録 if ( t2 < t1 ) then diff = (t_max - t1) + t2 + 1 else diff = t2 - t1 endif print "(A, F10.3)", "time it took was:", diff/dble(t_rate) IF (INFO .EQ. -1) THEN WRITE (*,*) 'not converge' RETURN ENDIF E = ZZERO DO I = 1, N DO J = 1, N S = B(I, J) DO K = 1, N S = S - LAMBDA(K)*W(I,K)*DCONJG(W(J,K)) ENDDO E = E+DCONJG(S)*S ENDDO ENDDO WRITE (*,*) 'ERROR',SQRT(DBLE(E)) CALL ZHERK ('U', 'C', N, N, DONE, W, N, DZERO, A, N) DO I=1,N A(I,I) = A(I,I) - ZONE ENDDO E = ZZERO DO I=2,N E = E+ZDOTC(I-1,A(1,I),1,A(1,I),1) ENDDO E = DTWO*E DO I=1,N E = E+DCONJG(A(I,I))*A(I,I) ENDDO WRITE (*,*) 'ERROR',SQRT(DBLE(E)) END PROGRAM TEDTDJE