SUBROUTINE DROTU(N, X, INCX, Y, INCY, F, S) IMPLICIT NONE INTEGER I,N,INCX,INCY DOUBLE PRECISION 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 DROTU SUBROUTINE DROTV(N, X, INCX, Y, INCY, C, G) IMPLICIT NONE INTEGER I,N,INCX,INCY DOUBLE PRECISION 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 DROTV SUBROUTINE DROTW(N, X, INCX, Y, INCY, C, G) IMPLICIT NONE INTEGER I,N,INCX,INCY DOUBLE PRECISION 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 DROTW 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 DJASVD(N,A,SIGMA,U,V,INFO) use omp_lib IMPLICIT NONE INTEGER N, INFO, FLAG, LL, MM, FLAG2, FLAG3, FLAG4 DOUBLE PRECISION A(N,N),U(N,N),V(N,N),SIGMA(N),B(N),Z(N) DOUBLE PRECISION ZERO, ONE, HALF PARAMETER (ONE = 1.0D0, ZERO = 0.0D0, HALF = 0.5D0) INTEGER I, J, K, ITER, MAXITER PARAMETER (MAXITER = 100) DOUBLE PRECISION T, X DOUBLE PRECISION CSL, SNL, CSR, SNR, T1, T2 DOUBLE PRECISION EPS, R1, R2, TMP DOUBLE PRECISION F1, F2, C1, S1, C2, S2, G1, G2 EXTERNAL DFMA0, DLAMCH DOUBLE PRECISION DFMA0, DLAMCH EPS = DLAMCH('P') DO J=1,N IF (A(J,J) .GE. ZERO) THEN SIGMA(J) = A(J,J) ELSE SIGMA(J) = -A(J,J) CALL DSCAL(J-1,-ONE,A(1,J),1) CALL DSCAL(N,-ONE,V(1,J),1) ENDIF ENDDO IF (ABS(SIGMA(1)) .GE. ABS(SIGMA(N))) THEN FLAG4 = 0 ELSE FLAG4 = 1 ENDIF LL = 1 MM = N DO ITER=1, MAXITER WRITE (*,*) 'ITER',ITER,LL,MM IF (FLAG4 .EQ. 0) THEN DO I=LL,MM-1 DO J=I+1,MM IF (ABS ( A(I,J) ) .GT. EPS*SQRT(ABS(SIGMA(I)))*SQRT(ABS(SIGMA(J)))) 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(SIGMA(I)))*SQRT(ABS(SIGMA(J)))) GO TO 20 ENDDO ENDDO 20 MM = I IF (LL .GE. MM) EXIT FLAG = 0 DO I = LL, MM B(I) = SIGMA(I) ENDDO DO I = LL, MM Z(I) = ZERO ENDDO DO J = LL, MM-1 DO K = J+1, MM IF (ABS ( A(J,K) ) .LE. EPS*SQRT(ABS(SIGMA(J)))*SQRT(ABS(SIGMA(K)))) THEN A(J,K) = ZERO ELSE F1 = SIGMA(J) - SIGMA(K) CALL DLARTGU (ABS(F1), ABS(A(J,K)), C1, S1, R1) F1 = F1 + SIGN(R1, F1) IF (F1 .GE. ZERO) THEN G1 = A(J,K) ELSE G1 = -A(J,K) F1 = -F1 ENDIF F2 = SIGMA(J) + SIGMA(K) CALL DLARTGU (ABS(F2), ABS(A(J,K)), C2, S2, R2) F2 = F2 + SIGN(R2, F2) IF (F2 .GE. ZERO) THEN G2 = -A(J,K) ELSE G2 = A(J,K) F2 = -F2 ENDIF IF (F1 .GE. F2) THEN T1 = G1/F1 C1 = DFMA0(-T1, G2, F2) S1 = DFMA0(T1, F2, G2) C2 = DFMA0(T1, G2, F2) S2 = DFMA0(T1, F2, -G2) ELSE T2 = G2/F2 C1 = DFMA0(-G1, T2, F1) S1 = DFMA0(F1, T2, G1) C2 = DFMA0(G1, T2, F1) S2 = DFMA0(-F1, T2, G1) ENDIF CALL DLARTGU (ABS(C1), ABS(S1), CSL, SNL, R1) CSL = SIGN(CSL, C1) SNL = SIGN(SNL, S1) CALL DLARTGU (ABS(C2), ABS(S2), CSR, SNR, R2) CSR = SIGN(CSR, C2) SNR = SIGN(SNR, S2) T = CSL + CSR R1 = SNR / T R2 = -SNL / T Z(J) = DFMA0( A(J,K), R1, Z(J)) Z(K) = DFMA0( A(J,K), R2, Z(K)) SIGMA(J) = B(J)+Z(J) SIGMA(K) = B(K)+Z(K) IF (ABS(SIGMA(J)) .LT. ABS(SIGMA(K))) THEN F1 = CSL F2 = SNL IF (F2 .GT. ZERO) THEN FLAG2 = 0 CSL = F2 SNL = -F1 ELSE FLAG2 = 1 CSL = -F2 SNL = F1 ENDIF F1 = CSR F2 = SNR IF (F2 .GT. ZERO) THEN FLAG3 = 0 CSR = F2 SNR = -F1 ELSE FLAG3 = 1 CSR = -F2 SNR = F1 ENDIF IF (FLAG2 .EQ. FLAG3) THEN TMP = SIGMA(J) SIGMA(J) = SIGMA(K) SIGMA(K) = TMP TMP = B(J) B(J) = B(K) B(K) = TMP TMP = Z(J) Z(J) = Z(K) Z(K) = TMP ELSE TMP = SIGMA(J) SIGMA(J) = -SIGMA(K) SIGMA(K) = -TMP TMP = B(J) B(J) = -B(K) B(K) = -TMP TMP = Z(J) Z(J) = -Z(K) Z(K) = -TMP ENDIF ENDIF IF (CSL .GE. ABS(SNL)) THEN X = SNL/(ONE+CSL) CSL = DFMA0(-SNL,X,ONE) IF (CSL .NE. ONE .OR. SNL .NE. ZERO) THEN FLAG = 1 CALL DROTU (J - LL, A(LL, J), 1, A(LL, K), 1, X, SNL) CALL DROTU (MM - K, A(J, K+1), N, A(K, K+1), N, X, SNL) CALL DROTU (N, U(1, J), 1, U(1, K), 1, X, SNL) ENDIF ELSE IF(SNL .GE. ZERO) THEN X = CSL/(ONE+SNL) SNL = DFMA0(-CSL,X,ONE) FLAG = 1 CALL DROTV (J - LL, A(LL, J), 1, A(LL, K), 1, CSL, X) CALL DROTV (MM - K, A(J, K+1), N, A(K, K+1), N, CSL, X) CALL DROTV (N, U(1, J), 1, U(1, K), 1, CSL, X) ELSE IF(SNL .LE. ZERO) THEN X = CSL/(ONE-SNL) SNL = DFMA0(CSL,X,-ONE) FLAG = 1 CALL DROTW (J - LL, A(LL, J), 1, A(LL, K), 1, CSL, X) CALL DROTW (MM - K, A(J, K+1), N, A(K, K+1), N, CSL, X) CALL DROTW (N, U(1, J), 1, U(1, K), 1, CSL, X) ENDIF IF (CSR .GE. ABS(SNR)) THEN X = SNR/(ONE+CSR) CSR = DFMA0(-SNR,X,ONE) IF (CSR .NE. ONE .OR. SNR .NE. ZERO) THEN FLAG = 1 CALL DROTU (K - J - 1, A(J,J + 1), N, A(J + 1,K), 1, X, SNR) CALL DROTU (N, V(1, J), 1, V(1, K), 1, X, SNR) ENDIF ELSE IF(SNR .GE. ZERO) THEN X = CSR/(ONE+SNR) SNR = DFMA0(-CSR,X,ONE) FLAG = 1 CALL DROTV (K - J - 1, A(J,J + 1), N, A(J + 1,K), 1, CSR, X) CALL DROTV (N, V(1, J), 1, V(1, K), 1, CSR, X) ELSE IF(SNR .LE. ZERO) THEN X = CSR/(ONE-SNR) SNR = DFMA0(CSR,X,-ONE) FLAG = 1 CALL DROTW (K - J - 1, A(J,J + 1), N, A(J + 1,K), 1, CSR, X) CALL DROTW (N, V(1, J), 1, V(1, K), 1, CSR, X) ENDIF A(J, K) = ZERO ENDIF END DO END DO IF (FLAG .EQ. 0) EXIT DO I=LL,MM-1 DO J=I+1,MM IF (ABS ( A(I,J) ) .GT. EPS*SQRT(ABS(SIGMA(I)))*SQRT(ABS(SIGMA(J)))) GO TO 30 ENDDO ENDDO 30 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(SIGMA(I)))*SQRT(ABS(SIGMA(J)))) GO TO 40 ENDDO ENDDO 40 MM = I IF (LL .GE. MM) EXIT FLAG = 0 DO I = LL, MM B(I) = SIGMA(I) ENDDO DO I = LL, MM Z(I) = ZERO ENDDO DO J = LL, MM-1 DO K = J+1, MM IF (ABS ( A(J,K) ) .LE. EPS*SQRT(ABS(SIGMA(J)))*SQRT(ABS(SIGMA(K)))) THEN A(J,K) = ZERO ELSE F1 = SIGMA(J) - SIGMA(K) CALL DLARTGU (ABS(F1), ABS(A(J,K)), C1, S1, R1) F1 = F1 + SIGN(R1, F1) IF (F1 .GE. ZERO) THEN G1 = A(J,K) ELSE G1 = -A(J,K) F1 = -F1 ENDIF F2 = SIGMA(J) + SIGMA(K) CALL DLARTGU (ABS(F2), ABS(A(J,K)), C2, S2, R2) F2 = F2 + SIGN(R2, F2) IF (F2 .GE. ZERO) THEN G2 = A(J,K) ELSE G2 = -A(J,K) F2 = -F2 ENDIF IF (F1 .GE. F2) THEN T1 = G1/F1 C1 = DFMA0(-T1, G2, F2) S1 = DFMA0(T1, F2, G2) C2 = DFMA0(T1, G2, F2) S2 = DFMA0(T1, F2, -G2) ELSE T2 = G2/F2 C1 = DFMA0(-G1, T2, F1) S1 = DFMA0(F1, T2, G1) C2 = DFMA0(G1, T2, F1) S2 = DFMA0(-F1, T2, G1) ENDIF CALL DLARTGU (ABS(C1), ABS(S1), CSL, SNL, R1) CSL = SIGN(CSL, C1) SNL = SIGN(SNL, S1) CALL DLARTGU (ABS(C2), ABS(S2), CSR, SNR, R2) CSR = SIGN(CSR, C2) SNR = SIGN(SNR, S2) T = CSL + CSR R1 = SNL / T R2 = -SNR / T Z(J) = DFMA0( A(J,K), R1, Z(J)) Z(K) = DFMA0( A(J,K), R2, Z(K)) SIGMA(J) = B(J)+Z(J) SIGMA(K) = B(K)+Z(K) IF (ABS(SIGMA(J)) .LT. ABS(SIGMA(K))) THEN F1 = CSL F2 = SNL IF (F2 .GT. ZERO) THEN FLAG2 = 0 CSL = F2 SNL = -F1 ELSE FLAG2 = 1 CSL = -F2 SNL = F1 ENDIF F1 = CSR F2 = SNR IF (F2 .GT. ZERO) THEN FLAG3 = 0 CSR = F2 SNR = -F1 ELSE FLAG3 = 1 CSR = -F2 SNR = F1 ENDIF IF (FLAG2 .EQ. FLAG3) THEN TMP = SIGMA(J) SIGMA(J) = SIGMA(K) SIGMA(K) = TMP TMP = B(J) B(J) = B(K) B(K) = TMP TMP = Z(J) Z(J) = Z(K) Z(K) = TMP ELSE TMP = SIGMA(J) SIGMA(J) = -SIGMA(K) SIGMA(K) = -TMP TMP = B(J) B(J) = -B(K) B(K) = -TMP TMP = Z(J) Z(J) = -Z(K) Z(K) = -TMP ENDIF ENDIF IF (CSR .GE. ABS(SNR)) THEN X = SNR/(ONE+CSR) CSR = DFMA0(-SNR,X,ONE) IF (CSR .NE. ONE .OR. SNR .NE. ZERO) THEN FLAG = 1 CALL DROTU (J - LL, A(LL,J), 1, A(LL,K), 1, X, SNR) CALL DROTU (MM - K, A(J, K+1), N, A(K, K+1), N, X, SNR) CALL DROTU (N, V(1, J), 1, V(1, K), 1, X, SNR) ENDIF ELSE IF(SNR .GE. ZERO) THEN X = CSR/(ONE+SNR) SNR = DFMA0(-CSR,X,ONE) FLAG = 1 CALL DROTV (J - LL, A(LL,J), 1, A(LL,K), 1, CSR, X) CALL DROTV (MM - K, A(J, K+1), N, A(K, K+1), N, CSR, X) CALL DROTV (N, V(1, J), 1, V(1, K), 1, CSR, X) ELSE IF(SNR .LE. ZERO) THEN X = CSR/(ONE-SNR) SNR = DFMA0(CSR,X,-ONE) FLAG = 1 CALL DROTW (J - LL, A(LL,J), 1, A(LL,K), 1, CSR, X) CALL DROTW (MM - K, A(J, K+1), N, A(K, K+1), N, CSR, X) CALL DROTW (N, V(1, J), 1, V(1, K), 1, CSR, X) ENDIF IF (CSL .GE. ABS(SNL)) THEN X = SNL/(ONE+CSL) CSL = DFMA0(-SNL,X,ONE) IF (CSL .NE. ONE .OR. SNL .NE. ZERO) THEN FLAG = 1 CALL DROTU (K - J - 1, A(J, J+1), N, A(J+1, K), 1, X, SNL) CALL DROTU (N, U(1, J), 1, U(1, K), 1, X, SNL) ENDIF ELSE IF(SNL .GE. ZERO) THEN X = CSL/(ONE+SNL) SNL = DFMA0(-CSL,X,ONE) FLAG = 1 CALL DROTV (K - J - 1, A(J, J+1), N, A(J+1, K), 1, CSL, X) CALL DROTV (N, U(1, J), 1, U(1, K), 1, CSL, X) ELSE IF(SNL .LE. ZERO) THEN X = CSL/(ONE-SNL) SNL = DFMA0(CSL,X,-ONE) FLAG = 1 CALL DROTW (K - J - 1, A(J, J+1), N, A(J+1, K), 1, CSL, X) CALL DROTW (N, U(1, J), 1, U(1, K), 1, CSL, X) ENDIF A(J, K) = ZERO ENDIF ENDDO ENDDO IF (FLAG .EQ. 0) EXIT ELSE DO I=LL,MM-1 DO J=I+1,MM IF (ABS ( A(I,J) ) .GT. EPS*SQRT(ABS(SIGMA(I)))*SQRT(ABS(SIGMA(J)))) GO TO 50 ENDDO ENDDO 50 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(SIGMA(I)))*SQRT(ABS(SIGMA(J)))) GO TO 60 ENDDO ENDDO 60 MM = I IF (LL .GE. MM) EXIT FLAG = 0 DO I = LL, MM B(I) = SIGMA(I) ENDDO DO I = LL, MM Z(I) = ZERO ENDDO DO K = MM, LL+1, -1 DO J = K-1, LL, -1 IF (ABS ( A(J,K) ) .LE. EPS*SQRT(ABS(SIGMA(J)))*SQRT(ABS(SIGMA(K)))) THEN A(J,K) = ZERO ELSE F1 = SIGMA(J) - SIGMA(K) CALL DLARTGU (ABS(F1), ABS(A(J,K)), C1, S1, R1) F1 = F1 + SIGN(R1, F1) IF (F1 .GE. ZERO) THEN G1 = A(J,K) ELSE G1 = -A(J,K) F1 = -F1 ENDIF F2 = SIGMA(J) + SIGMA(K) CALL DLARTGU (ABS(F2), ABS(A(J,K)), C2, S2, R2) F2 = F2 + SIGN(R2, F2) IF (F2 .GE. ZERO) THEN G2 = -A(J,K) ELSE G2 = A(J,K) F2 = -F2 ENDIF IF (F1 .GE. F2) THEN T1 = G1/F1 C1 = DFMA0(-T1, G2, F2) S1 = DFMA0(T1, F2, G2) C2 = DFMA0(T1, G2, F2) S2 = DFMA0(T1, F2, -G2) ELSE T2 = G2/F2 C1 = DFMA0(-G1, T2, F1) S1 = DFMA0(F1, T2, G1) C2 = DFMA0(G1, T2, F1) S2 = DFMA0(-F1, T2, G1) ENDIF CALL DLARTGU (ABS(C1), ABS(S1), CSL, SNL, R1) CSL = SIGN(CSL, C1) SNL = SIGN(SNL, S1) CALL DLARTGU (ABS(C2), ABS(S2), CSR, SNR, R2) CSR = SIGN(CSR, C2) SNR = SIGN(SNR, S2) T = CSL + CSR R1 = SNR / T R2 = -SNL / T Z(J) = DFMA0( A(J,K), R1, Z(J)) Z(K) = DFMA0( A(J,K), R2, Z(K)) SIGMA(J) = B(J)+Z(J) SIGMA(K) = B(K)+Z(K) IF (ABS(SIGMA(J)) .GE. ABS(SIGMA(K))) THEN F1 = CSL F2 = SNL IF (F2 .GT. ZERO) THEN FLAG2 = 0 CSL = F2 SNL = -F1 ELSE FLAG2 = 1 CSL = -F2 SNL = F1 ENDIF F1 = CSR F2 = SNR IF (F2 .GT. ZERO) THEN FLAG3 = 0 CSR = F2 SNR = -F1 ELSE FLAG3 = 1 CSR = -F2 SNR = F1 ENDIF IF (FLAG2 .EQ. FLAG3) THEN TMP = SIGMA(J) SIGMA(J) = SIGMA(K) SIGMA(K) = TMP TMP = B(J) B(J) = B(K) B(K) = TMP TMP = Z(J) Z(J) = Z(K) Z(K) = TMP ELSE TMP = SIGMA(J) SIGMA(J) = -SIGMA(K) SIGMA(K) = -TMP TMP = B(J) B(J) = -B(K) B(K) = -TMP TMP = Z(J) Z(J) = -Z(K) Z(K) = -TMP ENDIF ENDIF IF (CSR .GE. ABS(SNR)) THEN X = SNR/(ONE+CSR) CSR = DFMA0(-SNR,X,ONE) IF (CSR .NE. ONE .OR. SNR .NE. ZERO) THEN FLAG = 1 CALL DROTU (J - LL, A(LL,J), 1, A(LL,K), 1, X, SNR) CALL DROTU (MM - K, A(J, K+1), N, A(K, K+1), N, X, SNR) CALL DROTU (N, V(1, J), 1, V(1, K), 1, X, SNR) ENDIF ELSE IF(SNR .GE. ZERO) THEN X = CSR/(ONE+SNR) SNR = DFMA0(-CSR,X,ONE) FLAG = 1 CALL DROTV (J - LL, A(LL,J), 1, A(LL,K), 1, CSR, X) CALL DROTV (MM - K, A(J, K+1), N, A(K, K+1), N, CSR, X) CALL DROTV (N, V(1, J), 1, V(1, K), 1, CSR, X) ELSE IF(SNR .LE. ZERO) THEN X = CSR/(ONE-SNR) SNR = DFMA0(CSR,X,-ONE) FLAG = 1 CALL DROTW (J - LL, A(LL,J), 1, A(LL,K), 1, CSR, X) CALL DROTW (MM - K, A(J, K+1), N, A(K, K+1), N, CSR, X) CALL DROTW (N, V(1, J), 1, V(1, K), 1, CSR, X) ENDIF IF (CSL .GE. ABS(SNL)) THEN X = SNL/(ONE+CSL) CSL = DFMA0(-SNL,X,ONE) IF (CSL .NE. ONE .OR. SNL .NE. ZERO) THEN FLAG = 1 CALL DROTU (K - J - 1, A(J, J+1), N, A(J+1, K), 1, X, SNL) CALL DROTU (N, U(1, J), 1, U(1, K), 1, X, SNL) ENDIF ELSE IF(SNL .GE. ZERO) THEN X = CSL/(ONE+SNL) SNL = DFMA0(-CSL,X,ONE) FLAG = 1 CALL DROTV (K - J - 1, A(J, J+1), N, A(J+1, K), 1, CSL, X) CALL DROTV (N, U(1, J), 1, U(1, K), 1, CSL, X) ELSE IF(SNL .LE. ZERO) THEN X = CSL/(ONE-SNL) SNL = DFMA0(CSL,X,-ONE) FLAG = 1 CALL DROTW (K - J - 1, A(J, J+1), N, A(J+1, K), 1, CSL, X) CALL DROTW (N, U(1, J), 1, U(1, K), 1, CSL, X) ENDIF A(J, K) = ZERO ENDIF ENDDO ENDDO IF (FLAG .EQ. 0) EXIT DO I=LL,MM-1 DO J=I+1,MM IF (ABS ( A(I,J) ) .GT. EPS*SQRT(ABS(SIGMA(I)))*SQRT(ABS(SIGMA(J)))) GO TO 70 ENDDO ENDDO 70 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(SIGMA(I)))*SQRT(ABS(SIGMA(J)))) GO TO 80 ENDDO ENDDO 80 MM = I IF (LL .GE. MM) EXIT FLAG = 0 DO I = LL, MM B(I) = SIGMA(I) ENDDO DO I = LL, MM Z(I) = ZERO ENDDO DO K = MM, LL+1, -1 DO J = K-1, LL, -1 IF (ABS ( A(J,K) ) .LE. EPS*SQRT(ABS(SIGMA(J)))*SQRT(ABS(SIGMA(K)))) THEN A(J,K) = ZERO ELSE F1 = SIGMA(J) - SIGMA(K) CALL DLARTGU (ABS(F1), ABS(A(J,K)), C1, S1, R1) F1 = F1 + SIGN(R1, F1) IF (F1 .GE. ZERO) THEN G1 = A(J,K) ELSE G1 = -A(J,K) F1 = -F1 ENDIF F2 = SIGMA(J) + SIGMA(K) CALL DLARTGU (ABS(F2), ABS(A(J,K)), C2, S2, R2) F2 = F2 + SIGN(R2, F2) IF (F2 .GE. ZERO) THEN G2 = A(J,K) ELSE G2 = -A(J,K) F2 = -F2 ENDIF IF (F1 .GE. F2) THEN T1 = G1/F1 C1 = DFMA0(-T1, G2, F2) S1 = DFMA0(T1, F2, G2) C2 = DFMA0(T1, G2, F2) S2 = DFMA0(T1, F2, -G2) ELSE T2 = G2/F2 C1 = DFMA0(-G1, T2, F1) S1 = DFMA0(F1, T2, G1) C2 = DFMA0(G1, T2, F1) S2 = DFMA0(-F1, T2, G1) ENDIF CALL DLARTGU (ABS(C1), ABS(S1), CSL, SNL, R1) CSL = SIGN(CSL, C1) SNL = SIGN(SNL, S1) CALL DLARTGU (ABS(C2), ABS(S2), CSR, SNR, R2) CSR = SIGN(CSR, C2) SNR = SIGN(SNR, S2) T = CSL + CSR R1 = SNL / T R2 = -SNR / T Z(J) = DFMA0( A(J,K), R1, Z(J)) Z(K) = DFMA0( A(J,K), R2, Z(K)) SIGMA(J) = B(J)+Z(J) SIGMA(K) = B(K)+Z(K) IF (ABS(SIGMA(J)) .GE. ABS(SIGMA(K))) THEN F1 = CSL F2 = SNL IF (F2 .GT. ZERO) THEN FLAG2 = 0 CSL = F2 SNL = -F1 ELSE FLAG2 = 1 CSL = -F2 SNL = F1 ENDIF F1 = CSR F2 = SNR IF (F2 .GT. ZERO) THEN FLAG3 = 0 CSR = F2 SNR = -F1 ELSE FLAG3 = 1 CSR = -F2 SNR = F1 ENDIF IF (FLAG2 .EQ. FLAG3) THEN TMP = SIGMA(J) SIGMA(J) = SIGMA(K) SIGMA(K) = TMP TMP = B(J) B(J) = B(K) B(K) = TMP TMP = Z(J) Z(J) = Z(K) Z(K) = TMP ELSE TMP = SIGMA(J) SIGMA(J) = -SIGMA(K) SIGMA(K) = -TMP TMP = B(J) B(J) = -B(K) B(K) = -TMP TMP = Z(J) Z(J) = -Z(K) Z(K) = -TMP ENDIF ENDIF IF (CSL .GE. ABS(SNL)) THEN X = SNL/(ONE+CSL) CSL = DFMA0(-SNL,X,ONE) IF (CSL .NE. ONE .OR. SNL .NE. ZERO) THEN FLAG = 1 CALL DROTU (J - LL, A(LL, J), 1, A(LL, K), 1, X, SNL) CALL DROTU (MM - K, A(J, K+1), N, A(K, K+1), N, X, SNL) CALL DROTU (N, U(1, J), 1, U(1, K), 1, X, SNL) ENDIF ELSE IF(SNL .GE. ZERO) THEN X = CSL/(ONE+SNL) SNL = DFMA0(-CSL,X,ONE) FLAG = 1 CALL DROTV (J - LL, A(LL, J), 1, A(LL, K), 1, CSL, X) CALL DROTV (MM - K, A(J, K+1), N, A(K, K+1), N, CSL, X) CALL DROTV (N, U(1, J), 1, U(1, K), 1, CSL, X) ELSE IF(SNL .LE. ZERO) THEN X = CSL/(ONE-SNL) SNL = DFMA0(CSL,X,-ONE) FLAG = 1 CALL DROTW (J - LL, A(LL, J), 1, A(LL, K), 1, CSL, X) CALL DROTW (MM - K, A(J, K+1), N, A(K, K+1), N, CSL, X) CALL DROTW (N, U(1, J), 1, U(1, K), 1, CSL, X) ENDIF IF (CSR .GE. ABS(SNR)) THEN X = SNR/(ONE+CSR) CSR = DFMA0(-SNR,X,ONE) IF (CSR .NE. ONE .OR. SNR .NE. ZERO) THEN FLAG = 1 CALL DROTU (K - J - 1, A(J,J + 1), N, A(J + 1,K), 1, X, SNR) CALL DROTU (N, V(1, J), 1, V(1, K), 1, X, SNR) ENDIF ELSE IF(SNR .GE. ZERO) THEN X = CSR/(ONE+SNR) SNR = DFMA0(-CSR,X,ONE) FLAG = 1 CALL DROTV (K - J - 1, A(J,J + 1), N, A(J + 1,K), 1, CSR, X) CALL DROTV (N, V(1, J), 1, V(1, K), 1, CSR, X) ELSE IF(SNR .LE. ZERO) THEN X = CSR/(ONE-SNR) SNR = DFMA0(CSR,X,-ONE) FLAG = 1 CALL DROTW (K - J - 1, A(J,J + 1), N, A(J + 1,K), 1, CSR, X) CALL DROTW (N, V(1, J), 1, V(1, K), 1, CSR, X) ENDIF A(J, K) = ZERO ENDIF END DO END DO IF (FLAG .EQ. 0) EXIT ENDIF ENDDO IF (ITER .GT. MAXITER) THEN INFO = -1 RETURN ENDIF DO I = 1,N IF (SIGMA(I) .LT. ZERO) THEN SIGMA(I) = -SIGMA(I) CALL DSCAL(N, -ONE, U(1,I), 1) 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 DSWAP(N, U(1,K), 1, U(1,I), 1) CALL DSWAP(N, V(1,K), 1, V(1,I), 1) ENDDO INFO = 0 RETURN END SUBROUTINE DJASVD PROGRAM TESTDJS use omp_lib IMPLICIT NONE INTEGER I, J, K, N, INFO, SEED DOUBLE PRECISION S, E, TMP DOUBLE PRECISION ZERO, ONE, TWO PARAMETER (ZERO = 0.0D0, ONE=1.0D0, TWO=2.0D0) DOUBLE PRECISION,ALLOCATABLE :: SIGMA(:) DOUBLE PRECISION,ALLOCATABLE :: A(:,:), B(:,:), U(:,:), V(:,:) EXTERNAL DFMA0, DDOT DOUBLE PRECISION DFMA0, DDOT EXTERNAL USERRAND DOUBLE PRECISION 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(SIGMA(N)) ALLOCATE(A(N,N), B(N,N), U(N,N), V(N,N)) DO I=1,SEED CALL XOR() ENDDO DO I=1,N DO J=1,I TMP = USERRAND() A(J,I) = TMP B(J,I) = TMP ENDDO DO J=I+1,N A(J,I) = ZERO B(J,I) = ZERO ENDDO ENDDO DO J=1,N DO K=1,N U(K,J) = ZERO ENDDO U(J,J) = ONE ENDDO DO J=1,N DO K=1,N V(K,J) = ZERO ENDDO V(J,J) = ONE ENDDO call system_clock(t1) ! 開始時を記録 CALL DJASVD(N,A,SIGMA,U,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 print "(A, F10.3)", "time it took was:", diff/dble(t_rate) IF (INFO .EQ. -1) THEN WRITE (*,*) 'not converge' RETURN ENDIF E = ZERO !$OMP PARALLEL DO PRIVATE(S,J,K) REDUCTION(+:E) DO I = 1, N DO J = 1, N S = B(I, J) DO K = 1, N S = S - SIGMA(K)*U(I,K)*V(J,K) ENDDO E = DFMA0(S, S, E) ENDDO ENDDO !$OMP END PARALLEL DO WRITE (*,*) 'ERROR',SQRT(E) CALL DSYRK ('U', 'T', N, N, ONE, U, N, ZERO, A, N) DO I=1,N A(I,I) = A(I,I) - ONE ENDDO E = ZERO !$OMP PARALLEL DO REDUCTION(+:E) DO I=2,N E = E+DDOT(I-1,A(1,I),1,A(1,I),1) ENDDO !$OMP END PARALLEL DO E = TWO*E DO I=1,N E = DFMA0(A(I,I),A(I,I),E) ENDDO WRITE (*,*) 'ERROR',SQRT(E) CALL DSYRK ('U', 'T', N, N, ONE, V, N, ZERO, A, N) DO I=1,N A(I,I) = A(I,I) - ONE ENDDO E = ZERO !$OMP PARALLEL DO REDUCTION(+:E) DO I=2,N E = E+DDOT(I-1,A(1,I),1,A(1,I),1) ENDDO !$OMP END PARALLEL DO E = TWO*E DO I=1,N E = DFMA0(A(I,I),A(I,I),E) ENDDO WRITE (*,*) 'ERROR',SQRT(E) END PROGRAM TESTDJS