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, T, 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 CJASVD(M,N,A,RANK,SIGMA,V,INFO) use omp_lib IMPLICIT NONE INTEGER M, N, INFO, IINFO, FLAG, RANK INTEGER IPIV1(N),IPIV2(N) REAL SIGMA(N) COMPLEX A(M,N),V(N,N),WORK(M),SC(N) REAL CS(N),SN(N) REAL SZERO, SONE, SHALF PARAMETER (SONE = 1.0E0, SZERO = 0.0E0, SHALF = 0.5E0) COMPLEX CZERO, CONE, CHALF PARAMETER (CONE = (1.0E0,0.0E0), CZERO = (0.0E0,0.0E0), & CHALF = (0.5E0,0.0E0)) INTEGER I, J, K, L, P, Q, R, ITER, MAXITER PARAMETER (MAXITER = 100) REAL T REAL TMP1, TMP2, TMP4, TMP5 REAL TMP10, TMP11, TMP12, TMP13, TMP14, TMP15 COMPLEX TMP6, TMP7, TMP8, TMP13C, TMP14C, W, B REAL H1, G1, R1, X, YR, YC, ZR, ZC, S1, S2 REAL C, S, C0, S0 REAL EPS, SAFMIN, SCALE1, SIGMX EXTERNAL SFMA0, SLAMCH, CDOTC, SCNRM2 REAL SFMA0, SLAMCH, SCNRM2 COMPLEX CDOTC EPS = SLAMCH('P') SAFMIN = SLAMCH('S') SCALE1 = SONE/REAL(M)/SQRT(SAFMIN) SIGMX = MAXVAL(ABS(A(1:M,1:N))) IF (SIGMX .GT. SZERO) THEN CALL CLASCL( 'G', 0, 0, SIGMX, SCALE1, M, N, & A, M, IINFO ) ENDIF DO I = 1,N SIGMA(I) = SCNRM2(M,A(1,I),1) ENDDO DO J = 1, N IPIV1(J) = 1 ENDDO DO ITER=1, MAXITER ! WRITE (*,*) 'ITER',ITER FLAG = 0 DO J = 1, N IPIV2(J) = 0 ENDDO DO J = 1, N-1 IF (IPIV1(J) .EQ. 0 .OR. SIGMA(J) .EQ. SZERO) CYCLE R = J DO L = J+1, N IF (IPIV1(L) .EQ. 1 .AND. SIGMA(L) .GT. SIGMA(R)) R = L ENDDO IF (R .NE. J) THEN CALL CSWAP(M,A(1,J),1,A(1,R),1) CALL SSWAP(1,SIGMA(J),1,SIGMA(R),1) L = IPIV1(R) IPIV1(R) = IPIV1(J) IPIV1(J) = L L = IPIV2(R) IPIV2(R) = IPIV2(J) IPIV2(J) = L ENDIF IF (IPIV1(J) .EQ. 0 .OR. SIGMA(J) .EQ. SZERO) THEN Q = 0 GO TO 40 ENDIF IF (SIGMA(J) .LT. SAFMIN) THEN TMP1 = SONE/SAFMIN ELSE TMP1 = SONE/SIGMA(J) ENDIF S1 = SIGMA(J)*TMP1 DO L = 1,M WORK(L) = A(L,J)*TMP1 ENDDO K = J 10 DO P = K+1, N IF (IPIV1(P) .EQ. 0 .OR. SIGMA(P) .EQ. SZERO) THEN SN(P)=SZERO SC(P)=CONE CYCLE ENDIF TMP8 = CDOTC(M,WORK(1),1,A(1,P),1) GO TO 20 ENDDO Q = N GO TO 40 20 K = P TMP4 = REAL(TMP8) TMP5 = AIMAG(TMP8) CALL SLARTGU (ABS(TMP4),ABS(TMP5),C,S,TMP2) IF (TMP2 .LT. SAFMIN) THEN IF (TMP2 .LE. SQRT(REAL(M))*EPS*SIGMA(K)*S1) THEN SN(K)=SZERO SC(K)=CONE GO TO 10 ENDIF ELSE IF (TMP2 .LE. EPS*SIGMA(K)*S1) THEN SN(K)=SZERO SC(K)=CONE GO TO 10 ENDIF ENDIF C = SIGN(C,TMP4) S = SIGN(S,TMP5) TMP6 = CMPLX(C,-S) S2 = SIGMA(K)*TMP1 H1 = (SIGMA(J) - SIGMA(K))*((S1+S2)*SHALF) G1 = TMP2 CALL SLARTGU(ABS(H1), ABS(G1), C, S, R1) T = G1 / (H1 + SIGN(R1,H1)) R1 = SQRT(SFMA0(T,T,SONE)) C = SONE/R1 S = T/R1 TMP10 = SIGMA(J)*SQRT(SFMA0(T,TMP2/(SIGMA(J)*S1),SONE)) IF (TMP10 .LT. SAFMIN) THEN TMP10 = SAFMIN ENDIF TMP1 = SONE/TMP10 TMP15 = SFMA0(-T,TMP2/(SIGMA(K)*S2),SONE) IF (TMP15 .GT. SZERO) THEN TMP11 = SIGMA(K)*SQRT(TMP15) IF (TMP11 .LT. SAFMIN) THEN TMP11 = SAFMIN ENDIF TMP12 = SONE/TMP11 ELSE TMP11 = SZERO ENDIF CS(K)=C SN(K)=S SC(K)=TMP6 IF (S .NE. SZERO) THEN X = S/(SONE+C) IPIV2(J) = 1 IPIV2(K) = 1 DO P = K+1, N IF (IPIV1(P) .EQ. 0 .OR. SIGMA(P) .EQ. SZERO) THEN SN(P)=SZERO SC(P)=CONE CYCLE ENDIF IF (TMP11 .GT. SZERO) THEN TMP13C = CZERO TMP8 = CZERO TMP14C = CZERO DO L = 1,M YR = REAL(A(L,J)) YC = AIMAG(A(L,J)) W = TMP6*A(L,K) ZR = REAL(W) ZC = AIMAG(W) A(L,J) = cmplx(S*(-X*YR+ZR)+YR,S*(-X*YC+ZC)+YC) WORK(L) = A(L,J)*TMP1 TMP13C = TMP13C+CONJG(WORK(L))*WORK(L) TMP8 = TMP8 + CONJG(WORK(L))*A(L,P) A(L,K) = cmplx(-S*(X*ZR+YR)+ZR,-S*(X*ZC+YC)+ZC) B = A(L,K)*TMP12 TMP14C = TMP14C+CONJG(B)*B ENDDO TMP13 = REAL(TMP13C) TMP14 = REAL(TMP14C) SIGMA(J) = TMP10*SQRT(TMP13) S1 = SIGMA(J)*TMP1 SIGMA(K) = TMP11*SQRT(TMP14) ELSE TMP13C = CZERO TMP8 = CZERO DO L = 1,M YR = REAL(A(L,J)) YC = AIMAG(A(L,J)) W = TMP6*A(L,K) ZR = REAL(W) ZC = AIMAG(W) A(L,J) = cmplx(S*(-X*YR+ZR)+YR,S*(-X*YC+ZC)+YC) WORK(L) = A(L,J)*TMP1 TMP13C = TMP13C+CONJG(WORK(L))*WORK(L) TMP8 = TMP8 + CONJG(WORK(L))*A(L,P) A(L,K) = cmplx(-S*(X*ZR+YR)+ZR,-S*(X*ZC+YC)+ZC) ENDDO TMP13 = REAL(TMP13C) SIGMA(J) = TMP10*SQRT(TMP13) S1 = SIGMA(J)*TMP1 SIGMA(K) = SCNRM2(M,A(1,K),1) ENDIF GO TO 30 ENDDO IF (TMP11 .GT. SZERO) THEN TMP13C = CZERO TMP14C = CZERO DO L = 1,M YR = REAL(A(L,J)) YC = AIMAG(A(L,J)) W = TMP6*A(L,K) ZR = REAL(W) ZC = AIMAG(W) A(L,J) = cmplx(S*(-X*YR+ZR)+YR,S*(-X*YC+ZC)+YC) WORK(L) = A(L,J)*TMP1 TMP13C = TMP13C+CONJG(WORK(L))*WORK(L) A(L,K) = cmplx(-S*(X*ZR+YR)+ZR,-S*(X*ZC+YC)+ZC) B = A(L,K)*TMP12 TMP14C = TMP14C+CONJG(B)*B ENDDO TMP13 = REAL(TMP13C) TMP14 = REAL(TMP14C) SIGMA(J) = TMP10*SQRT(TMP13) S1 = SIGMA(J)*TMP1 SIGMA(K) = TMP11*SQRT(TMP14) ELSE TMP13C = CZERO DO L = 1,M YR = REAL(A(L,J)) YC = AIMAG(A(L,J)) W = TMP6*A(L,K) ZR = REAL(W) ZC = AIMAG(W) A(L,J) = cmplx(S*(-X*YR+ZR)+YR,S*(-X*YC+ZC)+YC) WORK(L) = A(L,J)*TMP1 TMP13C = TMP13C+CONJG(WORK(L))*WORK(L) A(L,K) = cmplx(-S*(X*ZR+YR)+ZR,-S*(X*ZC+YC)+ZC) ENDDO TMP13 = REAL(TMP13C) SIGMA(J) = TMP10*SQRT(TMP13) S1 = SIGMA(J)*TMP1 SIGMA(K) = SCNRM2(M,A(1,K),1) ENDIF 30 IF (FLAG .EQ. 0) THEN TMP7 = CDOTC(M,WORK(1),1,A(1,K),1) TMP13 = REAL(TMP7) TMP14 = AIMAG(TMP7) CALL SLARTGU(ABS(TMP13),ABS(TMP14),C0,S0,TMP4) IF (TMP4 .LT. SAFMIN) THEN IF (TMP4 .LE. SQRT(REAL(M))*EPS*SIGMA(K)*S1) FLAG = 1 ELSE IF (TMP4 .LE. EPS*SIGMA(K)*S1) FLAG = 1 ENDIF ENDIF IF (SIGMA(J) .GT. SZERO .AND. P .LE. N) GO TO 20 Q = K ELSE IF(TMP6 .NE. CONE) THEN IPIV2(K) = 1 IF (TMP11 .GT. SZERO) THEN TMP14C = CZERO DO L = 1,M A(L,K) = TMP6*A(L,K) B = A(L,K)*TMP12 TMP14C = TMP14C+CONJG(B)*B ENDDO TMP14 = REAL(TMP14C) SIGMA(K) = TMP11*SQRT(TMP14) ELSE CALL CSCAL(M,TMP6,A(1,K),1) SIGMA(K) = SCNRM2(M,A(1,K),1) ENDIF IF (FLAG .EQ. 0) THEN TMP7 = CDOTC(M,WORK(1),1,A(1,K),1) TMP13 = REAL(TMP7) TMP14 = AIMAG(TMP7) CALL SLARTGU(ABS(TMP13),ABS(TMP14),C0,S0,TMP4) IF (TMP4 .LT. SAFMIN) THEN IF (TMP4 .LE. SQRT(REAL(M))*EPS*SIGMA(K)*S1) FLAG = 1 ELSE IF (TMP4 .LE. EPS*SIGMA(K)*S1) FLAG = 1 ENDIF ENDIF GO TO 10 ELSE GO TO 10 ENDIF 40 IF (R .NE. J) THEN CALL CSWAP(N,V(1,J),1,V(1,R),1) ENDIF DO P=J+1,Q S=SN(P) TMP6=SC(P) IF (S .NE. SZERO) THEN C=CS(P) X = S/(SONE+C) DO L = 1,N YR = REAL(V(L,J)) YC = AIMAG(V(L,J)) W = TMP6*V(L,P) ZR = REAL(W) ZC = AIMAG(W) V(L,J) = cmplx(S*(-X*YR+ZR)+YR,S*(-X*YC+ZC)+YC) V(L,P) = cmplx(-S*(X*ZR+YR)+ZR,-S*(X*ZC+YC)+ZC) ENDDO ELSE IF(TMP6 .NE. CONE) THEN CALL CSCAL(N,TMP6,V(1,P),1) ENDIF ENDDO ENDDO IF (FLAG .EQ. 0) EXIT DO J = 1, N IPIV1(J) = IPIV2(J) ENDDO ENDDO IF (ITER .GT. MAXITER) THEN INFO = -1 RETURN ENDIF DO I = 1,N IF (SIGMA(I) .GT. SZERO) THEN DO L = 1, M A(L,I) = A(L,I)/SIGMA(I) ENDDO ENDIF ENDDO DO I = 1, N-1 K = I T = SIGMA(K) DO J = I+1, N IF (SIGMA(J) .GT. T) THEN K = J T = SIGMA(K) ENDIF ENDDO SIGMA(K) = SIGMA(I) SIGMA(I) = T CALL CSWAP(M, A(1,K), 1, A(1,I), 1) CALL CSWAP(N, V(1,K), 1, V(1,I), 1) ENDDO DO I = 1,N IF(SIGMA(I) .LT. SAFMIN) THEN RANK = I-1 GO TO 50 ENDIF ENDDO RANK = N 50 IF (SIGMX .GT. SZERO) THEN CALL SLASCL( 'G', 0, 0, SCALE1, SIGMX, N, 1, SIGMA, N, IINFO ) ENDIF INFO = 0 RETURN END SUBROUTINE CJASVD PROGRAM TESTSJS use omp_lib IMPLICIT NONE INTEGER I, J, K, M, N, INFO, SEED, RANK REAL TMP1, TMP2 COMPLEX E COMPLEX CZERO, CONE PARAMETER (CZERO = (0.0E0,0.0E0), CONE=(1.0E0,0.0E0)) REAL SZERO, SONE, STWO, SAFMIN PARAMETER (SZERO = 0.0E0, SONE=1.0E0, STWO=2.0E0) REAL,ALLOCATABLE :: SIGMA(:) COMPLEX,ALLOCATABLE :: A(:,:), B(:,:), WORK(:,:), V(:,:) REAL SLAMCH EXTERNAL USERRAND, SLAMCH DOUBLE PRECISION USERRAND EXTERNAL CDOTC COMPLEX CDOTC character*1000 ARG1,ARG2,ARG3 integer t1, t2, t_rate, t_max, diff SAFMIN = SLAMCH('S') CALL GETARG(1,ARG1) CALL GETARG(2,ARG2) CALL GETARG(3,ARG3) read (ARG1,*) M read (ARG2,*) N read (ARG3,*) SEED ALLOCATE(SIGMA(N)) ALLOCATE(A(M,N), B(M,N), WORK(N,N), V(N,N)) DO I=1,SEED CALL XOR() ENDDO DO I=1,N DO J=1,M TMP1 = REAL(USERRAND()) TMP2 = REAL(USERRAND()) A(J,I) = CMPLX(TMP1,TMP2) B(J,I) = A(J,I) ENDDO ENDDO DO J=1,N DO K=1,N V(K,J) = CZERO ENDDO V(J,J) = CONE ENDDO call system_clock(t1) ! 開始時を記録 CALL CJASVD(M,N,A,RANK,SIGMA,V,INFO) call system_clock(t2, t_rate, t_max) ! 終了時を記録 if ( t2 < t1 ) then diff = (t_max - t1) + t2 + 1 else diff = t2 - t1 endif WRITE (*,*) 'RANK',RANK print "(A, F10.3)", "time it took was:", diff/dble(t_rate) IF (INFO .EQ. -1) THEN WRITE (*,*) 'not converge' RETURN ENDIF !$OMP PARALLEL DO PRIVATE(K,I) DO J = 1, N DO K = 1, RANK DO I = 1, M B(I, J) = B(I, J) - SIGMA(K)*A(I,K)*CONJG(V(J,K)) ENDDO ENDDO ENDDO !$OMP END PARALLEL DO E = CZERO !$OMP PARALLEL DO PRIVATE(I) REDUCTION(+:E) DO J = 1, N DO I = 1, M E = E+CONJG(B(I, J))*B(I, J) ENDDO ENDDO !$OMP END PARALLEL DO WRITE (*,*) 'ERROR',SQRT(REAL(E)) CALL CHERK ('U', 'C', RANK, M, SONE, A, M, SZERO, WORK, N) DO I=1,RANK WORK(I,I) = WORK(I,I) - CONE ENDDO E = CZERO !$OMP PARALLEL DO REDUCTION(+:E) DO I=2,RANK E = E+CDOTC(I-1,WORK(1,I),1,WORK(1,I),1) ENDDO !$OMP END PARALLEL DO E = STWO*E DO I=1,RANK E = E+CONJG(WORK(I,I))*WORK(I,I) ENDDO WRITE (*,*) 'ERROR',SQRT(REAL(E)) CALL CHERK ('U', 'C', RANK, N, SONE, V, N, SZERO, WORK, N) DO I=1,RANK WORK(I,I) = WORK(I,I) - CONE ENDDO E = CZERO !$OMP PARALLEL DO REDUCTION(+:E) DO I=2,RANK E = E+CDOTC(I-1,WORK(1,I),1,WORK(1,I),1) ENDDO !$OMP END PARALLEL DO E = STWO*E DO I=1,RANK E = E+CONJG(WORK(I,I))*WORK(I,I) ENDDO WRITE (*,*) 'ERROR',SQRT(REAL(E)) END PROGRAM TESTSJS