SUBROUTINE SROTU(N, X, INCX, Y, INCY, F, S) IMPLICIT NONE INTEGER I,N,INCX,INCY REAL X(*),Y(*),F,S,U,V IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN DO I=1,N U=X(I) V=Y(I) X(I)=S*(-F*U+V)+U Y(I)=-S*(F*V+U)+V ENDDO ELSE DO I=1,N U=X((I-1)*INCX+1) V=Y((I-1)*INCY+1) X((I-1)*INCX+1)=S*(-F*U+V)+U Y((I-1)*INCY+1)=-S*(F*V+U)+V ENDDO ENDIF RETURN END SUBROUTINE SROTU SUBROUTINE SROTV(N, X, INCX, Y, INCY, C, G) IMPLICIT NONE INTEGER I,N,INCX,INCY REAL X(*),Y(*),G,C,U,V IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN DO I=1,N U=X(I) V=Y(I) X(I)=C*(-G*V+U)+V Y(I)=C*(G*U+V)-U ENDDO ELSE DO I=1,N U=X((I-1)*INCX+1) V=Y((I-1)*INCY+1) X((I-1)*INCX+1)=C*(-G*V+U)+V Y((I-1)*INCY+1)=C*(G*U+V)-U ENDDO ENDIF RETURN END SUBROUTINE SROTV SUBROUTINE SROTW(N, X, INCX, Y, INCY, C, G) IMPLICIT NONE INTEGER I,N,INCX,INCY REAL X(*),Y(*),G,C,U,V IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN DO I=1,N U=X(I) V=Y(I) X(I)=C*(G*V+U)-V Y(I)=C*(-G*U+V)+U ENDDO ELSE DO I=1,N U=X((I-1)*INCX+1) V=Y((I-1)*INCY+1) X((I-1)*INCX+1)=C*(G*V+U)-V Y((I-1)*INCY+1)=C*(-G*U+V)+U ENDDO ENDIF RETURN END SUBROUTINE SROTW 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 SJAEVD(N,A,LAMBDA,W,INFO) use omp_lib IMPLICIT NONE INTEGER N, INFO, FLAG, LL, MM REAL A(N,N),W(N,N),LAMBDA(N),B(N),Z(N) REAL ZERO, ONE, HALF PARAMETER (ONE = 1.0E0, ZERO = 0.0E0, HALF = 0.5E0) INTEGER I, J, K, ITER, MAXITER PARAMETER (MAXITER = 100) REAL T, U, C, S, EPS REAL THETA1, THETA2, TMP REAL G1, H1, R1, F EXTERNAL SFMA0, SLAMCH REAL SFMA0, SLAMCH EPS = SLAMCH('P') DO I = 1,N LAMBDA(I) = 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) = ZERO 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 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 SSWAP(J - LL, A(LL,J), 1, A(LL,K), 1) CALL SSWAP(K - J - 1, A(J,J+1), N, A(J+1,K), 1) CALL SSWAP(MM - K, A(J,K+1), N, A(K,K+1), N) CALL SSWAP(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) = ZERO ELSE H1 = HALF * LAMBDA(J) - HALF * LAMBDA(K) G1 = A(J,K) CALL SLARTGU (ABS(H1), ABS(G1), C, S, R1) T = G1 / (H1 + SIGN(R1,H1)) R1 = SQRT(SFMA0(T,T,ONE)) C = ONE/R1 S = T/R1 Z(J) = SFMA0 (T, A(J,K), Z(J)) Z(K) = SFMA0 (-T, A(J,K), 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. ZERO) 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/(ONE+C) C = SFMA0(-S,F,ONE) IF (C .NE. ONE .OR. S .NE. ZERO) THEN FLAG = 1 CALL SROTU(J - LL, A(LL,J), 1, A(LL,K), 1, F, S) CALL SROTU(K - J - 1, A(J,J+1), N, A(J+1,K), 1, F, S) CALL SROTU(MM - K, A(J,K+1), N, A(K,K+1), N, F, S) CALL SROTU(N, W(1,J), 1, W(1,K), 1, F, S) ENDIF ELSE IF(S .GE. ZERO) THEN F = C/(ONE+S) S = SFMA0(-C,F,ONE) FLAG = 1 CALL SROTV(J - LL, A(LL,J), 1, A(LL,K), 1, C, F) CALL SROTV(K - J - 1, A(J,J+1), N, A(J+1,K), 1, C, F) CALL SROTV(MM - K, A(J,K+1), N, A(K,K+1), N, C, F) CALL SROTV(N, W(1,J), 1, W(1,K), 1, C, F) ELSE IF(S .LE. ZERO) THEN F = C/(ONE-S) S = SFMA0(C,F,-ONE) FLAG = 1 CALL SROTW(J - LL, A(LL,J), 1, A(LL,K), 1, C, F) CALL SROTW(K - J - 1, A(J,J+1), N, A(J+1,K), 1, C, F) CALL SROTW(MM - K, A(J,K+1), N, A(K,K+1), N, C, F) CALL SROTW(N, W(1,J), 1, W(1,K), 1, C, F) ENDIF A(J,K) = ZERO 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 SSWAP(N, W(1,K), 1, W(1,I), 1) ENDDO INFO = 0 RETURN END SUBROUTINE SJAEVD PROGRAM TEDTDJE use omp_lib IMPLICIT NONE INTEGER I, J, K, N, INFO, SEED REAL S, E, TMP REAL ZERO, ONE, TWO PARAMETER (ZERO = 0.0E0, ONE=1.0E0, TWO=2.0E0) REAL,ALLOCATABLE :: LAMBDA(:) REAL,ALLOCATABLE :: A(:,:), B(:,:), W(:,:) EXTERNAL SFMA0, SDOT REAL SFMA0, SDOT DOUBLE PRECISION USERRAND EXTERNAL USERRAND 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 IF (MOD(I+2*J,2) .EQ. 0) THEN TMP = REAL(USERRAND()) ELSE TMP = -REAL(USERRAND()) ENDIF A(J,I) = TMP B(J,I) = TMP B(I,J) = TMP ENDDO ENDDO DO J=1,N DO K=1,N W(K,J) = ZERO ENDDO W(J,J) = ONE ENDDO call system_clock(t1) ! 開始時を記録 CALL SJAEVD(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 = ZERO DO I = 1, N DO J = 1, N S = B(I, J) DO K = 1, N S = S - LAMBDA(K)*W(I,K)*W(J,K) ENDDO E = SFMA0(S, S, E) ENDDO ENDDO WRITE (*,*) 'ERROR',SQRT(E) CALL SSYRK ('U', 'T', N, N, ONE, W, N, ZERO, A, N) DO I=1,N A(I,I) = A(I,I) - ONE ENDDO E = ZERO DO I=2,N E = E+SDOT(I-1,A(1,I),1,A(1,I),1) ENDDO E = TWO*E DO I=1,N E = SFMA0(A(I,I),A(I,I),E) ENDDO WRITE (*,*) 'ERROR',SQRT(E) END PROGRAM TEDTDJE