SUBROUTINE CSROTU(N, X, INCX, Y, INCY, F, S, T) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX X(*),Y(*),T,W REAL F,S,VR,VC,UR,UC IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN DO I=1,N UR=REAL(X(I)) UC=AIMAG(X(I)) W=Y(I)*T VR=REAL(W) VC=AIMAG(W) X(I)=CMPLX(S*(-F*UR+VR)+UR,S*(-F*UC+VC)+UC) Y(I)=CMPLX(-S*(F*VR+UR)+VR,-S*(F*VC+UC)+VC) ENDDO ELSE DO I=1,N UR=REAL(X((I-1)*INCX+1)) UC=AIMAG(X((I-1)*INCX+1)) W=Y((I-1)*INCY+1)*T VR=REAL(W) VC=AIMAG(W) X((I-1)*INCX+1)=CMPLX(S*(-F*UR+VR)+UR,S*(-F*UC+VC)+UC) Y((I-1)*INCY+1)=CMPLX(-S*(F*VR+UR)+VR,-S*(F*VC+UC)+VC) ENDDO ENDIF RETURN END SUBROUTINE CSROTU SUBROUTINE CSROTU2(N, X, INCX, Y, INCY, F, S, T) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX X(*),Y(*),T,W REAL F,S,VR,VC,UR,UC IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN DO I=1,N UR=REAL(X(I)) UC=-AIMAG(X(I)) W=Y(I)*T VR=REAL(W) VC=AIMAG(W) X(I)=CMPLX(S*(-F*UR+VR)+UR,S*(F*UC-VC)-UC) Y(I)=CMPLX(-S*(F*VR+UR)+VR,-S*(F*VC+UC)+VC) ENDDO ELSE DO I=1,N UR=REAL(X((I-1)*INCX+1)) UC=-AIMAG(X((I-1)*INCX+1)) W=Y((I-1)*INCY+1)*T VR=REAL(W) VC=AIMAG(W) X((I-1)*INCX+1)=CMPLX(S*(-F*UR+VR)+UR,S*(F*UC-VC)-UC) Y((I-1)*INCY+1)=CMPLX(-S*(F*VR+UR)+VR,-S*(F*VC+UC)+VC) ENDDO ENDIF RETURN END SUBROUTINE CSROTU2 SUBROUTINE CSROTV(N, X, INCX, Y, INCY, C, G, T) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX X(*),Y(*),T,W REAL G,C,VR,VC,UR,UC IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN DO I=1,N UR=REAL(X(I)) UC=AIMAG(X(I)) W=Y(I)*T VR=REAL(W) VC=AIMAG(W) X(I)=CMPLX(C*(-G*VR+UR)+VR,C*(-G*VC+UC)+VC) Y(I)=CMPLX(C*(G*UR+VR)-UR,C*(G*UC+VC)-UC) ENDDO ELSE DO I=1,N UR=REAL(X((I-1)*INCX+1)) UC=AIMAG(X((I-1)*INCX+1)) W=Y((I-1)*INCY+1)*T VR=REAL(W) VC=AIMAG(W) X((I-1)*INCX+1)=CMPLX(C*(-G*VR+UR)+VR,C*(-G*VC+UC)+VC) Y((I-1)*INCY+1)=CMPLX(C*(G*UR+VR)-UR,C*(G*UC+VC)-UC) ENDDO ENDIF RETURN END SUBROUTINE CSROTV SUBROUTINE CSROTV2(N, X, INCX, Y, INCY, C, G, T) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX X(*),Y(*),T,W REAL G,C,VR,VC,UR,UC IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN DO I=1,N UR=REAL(X(I)) UC=-AIMAG(X(I)) W=Y(I)*T VR=REAL(W) VC=AIMAG(W) X(I)=CMPLX(C*(-G*VR+UR)+VR,C*(G*VC-UC)-VC) Y(I)=CMPLX(C*(G*UR+VR)-UR,C*(G*UC+VC)-UC) ENDDO ELSE DO I=1,N UR=REAL(X((I-1)*INCX+1)) UC=-AIMAG(X((I-1)*INCX+1)) W=Y((I-1)*INCY+1)*T VR=REAL(W) VC=AIMAG(W) X((I-1)*INCX+1)=CMPLX(C*(-G*VR+UR)+VR,C*(G*VC-UC)-VC) Y((I-1)*INCY+1)=CMPLX(C*(G*UR+VR)-UR,C*(G*UC+VC)-UC) ENDDO ENDIF RETURN END SUBROUTINE CSROTV2 SUBROUTINE CSROTW(N, X, INCX, Y, INCY, C, G, T) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX X(*),Y(*),T,W REAL G,C,VR,VC,UR,UC IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN DO I=1,N UR=REAL(X(I)) UC=AIMAG(X(I)) W=Y(I)*T VR=REAL(W) VC=AIMAG(W) X(I)=CMPLX(C*(G*VR+UR)-VR,C*(G*VC+UC)-VC) Y(I)=CMPLX(C*(-G*UR+VR)+UR,C*(-G*UC+VC)+UC) ENDDO ELSE DO I=1,N UR=REAL(X((I-1)*INCX+1)) UC=AIMAG(X((I-1)*INCX+1)) W=Y((I-1)*INCY+1)*T VR=REAL(W) VC=AIMAG(W) X((I-1)*INCX+1)=CMPLX(C*(G*VR+UR)-VR,C*(G*VC+UC)-VC) Y((I-1)*INCY+1)=CMPLX(C*(-G*UR+VR)+UR,C*(-G*UC+VC)+UC) ENDDO ENDIF RETURN END SUBROUTINE CSROTW SUBROUTINE CSROTW2(N, X, INCX, Y, INCY, C, G, T) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX X(*),Y(*),T,W REAL G,C,VR,VC,UR,UC IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN DO I=1,N UR=REAL(X(I)) UC=-AIMAG(X(I)) W=Y(I)*T VR=REAL(W) VC=AIMAG(W) X(I)=CMPLX(C*(G*VR+UR)-VR,-C*(G*VC+UC)+VC) Y(I)=CMPLX(C*(-G*UR+VR)+UR,C*(-G*UC+VC)+UC) ENDDO ELSE DO I=1,N UR=REAL(X((I-1)*INCX+1)) UC=-AIMAG(X((I-1)*INCX+1)) W=Y((I-1)*INCY+1)*T VR=REAL(W) VC=AIMAG(W) X((I-1)*INCX+1)=CMPLX(C*(G*VR+UR)-VR,-C*(G*VC+UC)+VC) Y((I-1)*INCY+1)=CMPLX(C*(-G*UR+VR)+UR,C*(-G*UC+VC)+UC) ENDDO ENDIF RETURN END SUBROUTINE CSROTW2 SUBROUTINE SLARTGU (F1, G1, C, S, R1) IMPLICIT NONE REAL ZERO, ONE, HALF PARAMETER (ONE = 1.0E0, ZERO = 0.0E0, HALF = 0.5E0) REAL F1, G1, C, S, R1, U, V EXTERNAL SFMA0 REAL SFMA0 IF (F1 .GE. G1) THEN U = G1/F1 R1 = SQRT(SFMA0(U,U,ONE)) C = ONE/R1 S = U/R1 R1 = R1*F1 ELSE V = F1/G1 R1 = SQRT(SFMA0(V,V,ONE)) S = ONE/R1 C = V/R1 R1 = R1*G1 ENDIF RETURN END SUBROUTINE SLARTGU SUBROUTINE CJAEVD(N,A,LAMBDA,W,INFO) use omp_lib IMPLICIT NONE INTEGER N, INFO, FLAG, LL, MM REAL LAMBDA(N),TMP,B(N),Z(N),TMP3,TMP4 COMPLEX A(N,N),W(N,N) COMPLEX TMP2,TMP5 COMPLEX CZERO, CONE PARAMETER (CONE = (1.0E0,0.0E0), CZERO = (0.0E0,0.0E0)) REAL SZERO, SONE, SHALF PARAMETER (SONE = 1.0E0, SZERO = 0.0E0, SHALF = 0.5E0) INTEGER I, J, K, ITER, MAXITER PARAMETER (MAXITER = 100) REAL T, U, C, S, F, EPS REAL THETA1, THETA2 REAL G1, H1, R1 EXTERNAL SFMA0, SLAMCH REAL SFMA0, SLAMCH EPS = SLAMCH('P') DO I = 1,N LAMBDA(I) = REAL(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) = SZERO 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) = CONJG(A(J,K)) CALL SSWAP(1,LAMBDA(J),1,LAMBDA(K),1) CALL SSWAP(1,B(J),1,B(K),1) CALL SSWAP(1,Z(J),1,Z(K),1) CALL CSWAP(J - LL, A(LL,J), 1, A(LL,K), 1) DO I = J+1, K-1 TMP2 = CONJG(A(J,I)) A(J,I) = CONJG(A(I,K)) A(I,K) = TMP2 ENDDO CALL CSWAP(MM - K, A(J,K+1), N, A(K,K+1), N) CALL CSWAP(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) = CZERO ELSE TMP3 = REAL(A(J,K)) TMP4 = AIMAG(A(J,K)) CALL SLARTG(ABS(TMP3),ABS(TMP4),C,S,TMP) C = SIGN(C,TMP3) S = SIGN(S,TMP4) TMP2 = CMPLX(C,S) TMP5 = CONJG(TMP2) H1 = SHALF*LAMBDA(J) - SHALF*LAMBDA(K) G1 = TMP CALL SLARTGU (ABS(H1), G1, C, S, R1) T = G1 / (H1 + SIGN(R1,H1)) R1 = SQRT(SFMA0(T,T,SONE)) C = SONE/R1 S = T/R1 Z(J) = SFMA0 (T, TMP, Z(J)) Z(K) = SFMA0 (-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. SZERO) 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/(SONE+C) C = SFMA0(-S,F,SONE) IF (C .NE. SONE .OR. S .NE. SZERO) THEN FLAG = 1 CALL CSROTU(J - LL, A(LL,J), 1, A(LL,K), 1, F, S, TMP5) CALL CSROTU2(K - J - 1, A(J,J+1), N, A(J+1,K), 1, F, S, TMP5) CALL CSROTU(MM - K, A(J,K+1), N, A(K,K+1), N, F, S, TMP2) CALL CSROTU(N, W(1,J), 1, W(1,K), 1, F, S, TMP5) ELSE IF(TMP2 .NE. CONE .OR. TMP5 .NE. CONE) THEN FLAG = 1 CALL CSCAL(MM - K, TMP2, A(K, K+1), N) CALL CSCAL(K - LL, TMP5, A(LL, K), 1) CALL CSCAL(N, TMP5, W(1,K), 1) ENDIF ELSE IF(S .GE. SZERO) THEN F = C/(SONE+S) S = SFMA0(-C,F,SONE) FLAG = 1 CALL CSROTV(J - LL, A(LL,J), 1, A(LL,K), 1, C, F, TMP5) CALL CSROTV2(K - J - 1, A(J,J+1), N, A(J+1,K), 1, C, F, TMP5) CALL CSROTV(MM - K, A(J,K+1), N, A(K,K+1), N, C, F, TMP2) CALL CSROTV(N, W(1,J), 1, W(1,K), 1, C, F, TMP5) ELSE IF(S .LE. SZERO) THEN F = C/(SONE-S) S = SFMA0(C,F,-SONE) FLAG = 1 CALL CSROTW(J - LL, A(LL,J), 1, A(LL,K), 1, C, F, TMP5) CALL CSROTW2(K - J - 1, A(J,J+1), N, A(J+1,K), 1, C, F, TMP5) CALL CSROTW(MM - K, A(J,K+1), N, A(K,K+1), N, C, F, TMP2) CALL CSROTW(N, W(1,J), 1, W(1,K), 1, C, F, TMP5) ENDIF A(J,K) = CZERO 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 CSWAP(N, W(1,K), 1, W(1,I), 1) ENDDO INFO = 0 RETURN END SUBROUTINE CJAEVD PROGRAM TEDTDJE use omp_lib IMPLICIT NONE INTEGER I, J, K, N, INFO, SEED REAL TMP1, TMP2 COMPLEX E, S COMPLEX CZERO, CONE PARAMETER (CZERO = (0.0E0,0.0E0), CONE= (1.0E0,0.0E0)) REAL SZERO, SONE, STWO PARAMETER (SZERO = 0.0E0, SONE= 1.0E0, STWO= 2.0E0) REAL,ALLOCATABLE :: LAMBDA(:) COMPLEX,ALLOCATABLE :: A(:,:), B(:,:), W(:,:) DOUBLE PRECISION USERRAND EXTERNAL USERRAND EXTERNAL CDOTC COMPLEX CDOTC 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 = REAL(USERRAND()) TMP2 = REAL(USERRAND()) A(J,I) = cmplx(TMP1,TMP2) B(J,I) = A(J,I) B(I,J) = CONJG(B(J,I)) ENDDO TMP1 = REAL(USERRAND()) A(I,I) = cmplx(TMP1,SZERO) B(I,I) = A(I,I) ENDDO DO J=1,N DO K=1,N W(K,J) = CZERO ENDDO W(J,J) = CONE ENDDO call system_clock(t1) ! 開始時を記録 CALL CJAEVD(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 = CZERO DO I = 1, N DO J = 1, N S = B(I, J) DO K = 1, N S = S - LAMBDA(K)*W(I,K)*CONJG(W(J,K)) ENDDO E = E+CONJG(S)*S ENDDO ENDDO WRITE (*,*) 'ERROR',SQRT(REAL(E)) CALL CHERK ('U', 'C', N, N, SONE, W, N, SZERO, A, N) DO I=1,N A(I,I) = A(I,I) - CONE ENDDO E = CZERO DO I=2,N E = E+CDOTC(I-1,A(1,I),1,A(1,I),1) ENDDO E = STWO*E DO I=1,N E = E+CONJG(A(I,I))*A(I,I) ENDDO WRITE (*,*) 'ERROR',SQRT(REAL(E)) END PROGRAM TEDTDJE