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 SJASVD(N,A,SIGMA,U,V,INFO) use omp_lib IMPLICIT NONE INTEGER N, INFO, FLAG, LL, MM, FLAG2, FLAG3, FLAG4 REAL A(N,N),U(N,N),V(N,N),SIGMA(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, X REAL CSL, SNL, CSR, SNR, T1, T2 REAL EPS, R1, R2, TMP REAL F1, F2, C1, S1, C2, S2, G1, G2 EXTERNAL SFMA0, SLAMCH REAL SFMA0, SLAMCH EPS = SLAMCH('P') DO J=1,N IF (A(J,J) .GE. ZERO) THEN SIGMA(J) = A(J,J) ELSE SIGMA(J) = -A(J,J) CALL SSCAL(J-1,-ONE,A(1,J),1) CALL SSCAL(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 SLARTGU (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 SLARTGU (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 = SFMA0(-T1, G2, F2) S1 = SFMA0(T1, F2, G2) C2 = SFMA0(T1, G2, F2) S2 = SFMA0(T1, F2, -G2) ELSE T2 = G2/F2 C1 = SFMA0(-G1, T2, F1) S1 = SFMA0(F1, T2, G1) C2 = SFMA0(G1, T2, F1) S2 = SFMA0(-F1, T2, G1) ENDIF CALL SLARTGU (ABS(C1), ABS(S1), CSL, SNL, R1) CSL = SIGN(CSL, C1) SNL = SIGN(SNL, S1) CALL SLARTGU (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) = SFMA0( A(J,K), R1, Z(J)) Z(K) = SFMA0( 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 = SFMA0(-SNL,X,ONE) IF (CSL .NE. ONE .OR. SNL .NE. ZERO) THEN FLAG = 1 CALL SROTU (J - LL, A(LL, J), 1, A(LL, K), 1, X, SNL) CALL SROTU (MM - K, A(J, K+1), N, A(K, K+1), N, X, SNL) CALL SROTU (N, U(1, J), 1, U(1, K), 1, X, SNL) ENDIF ELSE IF(SNL .GE. ZERO) THEN X = CSL/(ONE+SNL) SNL = SFMA0(-CSL,X,ONE) FLAG = 1 CALL SROTV (J - LL, A(LL, J), 1, A(LL, K), 1, CSL, X) CALL SROTV (MM - K, A(J, K+1), N, A(K, K+1), N, CSL, X) CALL SROTV (N, U(1, J), 1, U(1, K), 1, CSL, X) ELSE IF(SNL .LE. ZERO) THEN X = CSL/(ONE-SNL) SNL = SFMA0(CSL,X,-ONE) FLAG = 1 CALL SROTW (J - LL, A(LL, J), 1, A(LL, K), 1, CSL, X) CALL SROTW (MM - K, A(J, K+1), N, A(K, K+1), N, CSL, X) CALL SROTW (N, U(1, J), 1, U(1, K), 1, CSL, X) ENDIF IF (CSR .GE. ABS(SNR)) THEN X = SNR/(ONE+CSR) CSR = SFMA0(-SNR,X,ONE) IF (CSR .NE. ONE .OR. SNR .NE. ZERO) THEN FLAG = 1 CALL SROTU (K - J - 1, A(J,J + 1), N, A(J + 1,K), 1, X, SNR) CALL SROTU (N, V(1, J), 1, V(1, K), 1, X, SNR) ENDIF ELSE IF(SNR .GE. ZERO) THEN X = CSR/(ONE+SNR) SNR = SFMA0(-CSR,X,ONE) FLAG = 1 CALL SROTV (K - J - 1, A(J,J + 1), N, A(J + 1,K), 1, CSR, X) CALL SROTV (N, V(1, J), 1, V(1, K), 1, CSR, X) ELSE IF(SNR .LE. ZERO) THEN X = CSR/(ONE-SNR) SNR = SFMA0(CSR,X,-ONE) FLAG = 1 CALL SROTW (K - J - 1, A(J,J + 1), N, A(J + 1,K), 1, CSR, X) CALL SROTW (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 SLARTGU (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 SLARTGU (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 = SFMA0(-T1, G2, F2) S1 = SFMA0(T1, F2, G2) C2 = SFMA0(T1, G2, F2) S2 = SFMA0(T1, F2, -G2) ELSE T2 = G2/F2 C1 = SFMA0(-G1, T2, F1) S1 = SFMA0(F1, T2, G1) C2 = SFMA0(G1, T2, F1) S2 = SFMA0(-F1, T2, G1) ENDIF CALL SLARTGU (ABS(C1), ABS(S1), CSL, SNL, R1) CSL = SIGN(CSL, C1) SNL = SIGN(SNL, S1) CALL SLARTGU (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) = SFMA0( A(J,K), R1, Z(J)) Z(K) = SFMA0( 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 = SFMA0(-SNR,X,ONE) IF (CSR .NE. ONE .OR. SNR .NE. ZERO) THEN FLAG = 1 CALL SROTU (J - LL, A(LL,J), 1, A(LL,K), 1, X, SNR) CALL SROTU (MM - K, A(J, K+1), N, A(K, K+1), N, X, SNR) CALL SROTU (N, V(1, J), 1, V(1, K), 1, X, SNR) ENDIF ELSE IF(SNR .GE. ZERO) THEN X = CSR/(ONE+SNR) SNR = SFMA0(-CSR,X,ONE) FLAG = 1 CALL SROTV (J - LL, A(LL,J), 1, A(LL,K), 1, CSR, X) CALL SROTV (MM - K, A(J, K+1), N, A(K, K+1), N, CSR, X) CALL SROTV (N, V(1, J), 1, V(1, K), 1, CSR, X) ELSE IF(SNR .LE. ZERO) THEN X = CSR/(ONE-SNR) SNR = SFMA0(CSR,X,-ONE) FLAG = 1 CALL SROTW (J - LL, A(LL,J), 1, A(LL,K), 1, CSR, X) CALL SROTW (MM - K, A(J, K+1), N, A(K, K+1), N, CSR, X) CALL SROTW (N, V(1, J), 1, V(1, K), 1, CSR, X) ENDIF IF (CSL .GE. ABS(SNL)) THEN X = SNL/(ONE+CSL) CSL = SFMA0(-SNL,X,ONE) IF (CSL .NE. ONE .OR. SNL .NE. ZERO) THEN FLAG = 1 CALL SROTU (K - J - 1, A(J, J+1), N, A(J+1, K), 1, X, SNL) CALL SROTU (N, U(1, J), 1, U(1, K), 1, X, SNL) ENDIF ELSE IF(SNL .GE. ZERO) THEN X = CSL/(ONE+SNL) SNL = SFMA0(-CSL,X,ONE) FLAG = 1 CALL SROTV (K - J - 1, A(J, J+1), N, A(J+1, K), 1, CSL, X) CALL SROTV (N, U(1, J), 1, U(1, K), 1, CSL, X) ELSE IF(SNL .LE. ZERO) THEN X = CSL/(ONE-SNL) SNL = SFMA0(CSL,X,-ONE) FLAG = 1 CALL SROTW (K - J - 1, A(J, J+1), N, A(J+1, K), 1, CSL, X) CALL SROTW (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 SLARTGU (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 SLARTGU (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 = SFMA0(-T1, G2, F2) S1 = SFMA0(T1, F2, G2) C2 = SFMA0(T1, G2, F2) S2 = SFMA0(T1, F2, -G2) ELSE T2 = G2/F2 C1 = SFMA0(-G1, T2, F1) S1 = SFMA0(F1, T2, G1) C2 = SFMA0(G1, T2, F1) S2 = SFMA0(-F1, T2, G1) ENDIF CALL SLARTGU (ABS(C1), ABS(S1), CSL, SNL, R1) CSL = SIGN(CSL, C1) SNL = SIGN(SNL, S1) CALL SLARTGU (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) = SFMA0( A(J,K), R1, Z(J)) Z(K) = SFMA0( 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 = SFMA0(-SNR,X,ONE) IF (CSR .NE. ONE .OR. SNR .NE. ZERO) THEN FLAG = 1 CALL SROTU (J - LL, A(LL,J), 1, A(LL,K), 1, X, SNR) CALL SROTU (MM - K, A(J, K+1), N, A(K, K+1), N, X, SNR) CALL SROTU (N, V(1, J), 1, V(1, K), 1, X, SNR) ENDIF ELSE IF(SNR .GE. ZERO) THEN X = CSR/(ONE+SNR) SNR = SFMA0(-CSR,X,ONE) FLAG = 1 CALL SROTV (J - LL, A(LL,J), 1, A(LL,K), 1, CSR, X) CALL SROTV (MM - K, A(J, K+1), N, A(K, K+1), N, CSR, X) CALL SROTV (N, V(1, J), 1, V(1, K), 1, CSR, X) ELSE IF(SNR .LE. ZERO) THEN X = CSR/(ONE-SNR) SNR = SFMA0(CSR,X,-ONE) FLAG = 1 CALL SROTW (J - LL, A(LL,J), 1, A(LL,K), 1, CSR, X) CALL SROTW (MM - K, A(J, K+1), N, A(K, K+1), N, CSR, X) CALL SROTW (N, V(1, J), 1, V(1, K), 1, CSR, X) ENDIF IF (CSL .GE. ABS(SNL)) THEN X = SNL/(ONE+CSL) CSL = SFMA0(-SNL,X,ONE) IF (CSL .NE. ONE .OR. SNL .NE. ZERO) THEN FLAG = 1 CALL SROTU (K - J - 1, A(J, J+1), N, A(J+1, K), 1, X, SNL) CALL SROTU (N, U(1, J), 1, U(1, K), 1, X, SNL) ENDIF ELSE IF(SNL .GE. ZERO) THEN X = CSL/(ONE+SNL) SNL = SFMA0(-CSL,X,ONE) FLAG = 1 CALL SROTV (K - J - 1, A(J, J+1), N, A(J+1, K), 1, CSL, X) CALL SROTV (N, U(1, J), 1, U(1, K), 1, CSL, X) ELSE IF(SNL .LE. ZERO) THEN X = CSL/(ONE-SNL) SNL = SFMA0(CSL,X,-ONE) FLAG = 1 CALL SROTW (K - J - 1, A(J, J+1), N, A(J+1, K), 1, CSL, X) CALL SROTW (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 SLARTGU (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 SLARTGU (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 = SFMA0(-T1, G2, F2) S1 = SFMA0(T1, F2, G2) C2 = SFMA0(T1, G2, F2) S2 = SFMA0(T1, F2, -G2) ELSE T2 = G2/F2 C1 = SFMA0(-G1, T2, F1) S1 = SFMA0(F1, T2, G1) C2 = SFMA0(G1, T2, F1) S2 = SFMA0(-F1, T2, G1) ENDIF CALL SLARTGU (ABS(C1), ABS(S1), CSL, SNL, R1) CSL = SIGN(CSL, C1) SNL = SIGN(SNL, S1) CALL SLARTGU (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) = SFMA0( A(J,K), R1, Z(J)) Z(K) = SFMA0( 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 = SFMA0(-SNL,X,ONE) IF (CSL .NE. ONE .OR. SNL .NE. ZERO) THEN FLAG = 1 CALL SROTU (J - LL, A(LL, J), 1, A(LL, K), 1, X, SNL) CALL SROTU (MM - K, A(J, K+1), N, A(K, K+1), N, X, SNL) CALL SROTU (N, U(1, J), 1, U(1, K), 1, X, SNL) ENDIF ELSE IF(SNL .GE. ZERO) THEN X = CSL/(ONE+SNL) SNL = SFMA0(-CSL,X,ONE) FLAG = 1 CALL SROTV (J - LL, A(LL, J), 1, A(LL, K), 1, CSL, X) CALL SROTV (MM - K, A(J, K+1), N, A(K, K+1), N, CSL, X) CALL SROTV (N, U(1, J), 1, U(1, K), 1, CSL, X) ELSE IF(SNL .LE. ZERO) THEN X = CSL/(ONE-SNL) SNL = SFMA0(CSL,X,-ONE) FLAG = 1 CALL SROTW (J - LL, A(LL, J), 1, A(LL, K), 1, CSL, X) CALL SROTW (MM - K, A(J, K+1), N, A(K, K+1), N, CSL, X) CALL SROTW (N, U(1, J), 1, U(1, K), 1, CSL, X) ENDIF IF (CSR .GE. ABS(SNR)) THEN X = SNR/(ONE+CSR) CSR = SFMA0(-SNR,X,ONE) IF (CSR .NE. ONE .OR. SNR .NE. ZERO) THEN FLAG = 1 CALL SROTU (K - J - 1, A(J,J + 1), N, A(J + 1,K), 1, X, SNR) CALL SROTU (N, V(1, J), 1, V(1, K), 1, X, SNR) ENDIF ELSE IF(SNR .GE. ZERO) THEN X = CSR/(ONE+SNR) SNR = SFMA0(-CSR,X,ONE) FLAG = 1 CALL SROTV (K - J - 1, A(J,J + 1), N, A(J + 1,K), 1, CSR, X) CALL SROTV (N, V(1, J), 1, V(1, K), 1, CSR, X) ELSE IF(SNR .LE. ZERO) THEN X = CSR/(ONE-SNR) SNR = SFMA0(CSR,X,-ONE) FLAG = 1 CALL SROTW (K - J - 1, A(J,J + 1), N, A(J + 1,K), 1, CSR, X) CALL SROTW (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 SSCAL(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 SSWAP(N, U(1,K), 1, U(1,I), 1) CALL SSWAP(N, V(1,K), 1, V(1,I), 1) ENDDO INFO = 0 RETURN END SUBROUTINE SJASVD PROGRAM TESTSJS 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 :: SIGMA(:) REAL,ALLOCATABLE :: A(:,:), B(:,:), U(:,:), V(:,:) EXTERNAL SFMA0, SDOT REAL SFMA0, SDOT 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 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 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 SJASVD(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 = SFMA0(S, S, E) ENDDO ENDDO !$OMP END PARALLEL DO WRITE (*,*) 'ERROR',SQRT(E) CALL SSYRK ('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+SDOT(I-1,A(1,I),1,A(1,I),1) ENDDO !$OMP END PARALLEL DO E = TWO*E DO I=1,N E = SFMA0(A(I,I),A(I,I),E) ENDDO WRITE (*,*) 'ERROR',SQRT(E) CALL SSYRK ('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+SDOT(I-1,A(1,I),1,A(1,I),1) ENDDO !$OMP END PARALLEL DO E = TWO*E DO I=1,N E = SFMA0(A(I,I),A(I,I),E) ENDDO WRITE (*,*) 'ERROR',SQRT(E) END PROGRAM TESTSJS