SUBROUTINE CSROTUx(N, X, INCX, Y, INCY, F, S, CS, SN) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX X(*),Y(*) REAL F,S,T,WR,WC,VR,VC,UR,UC REAL CS,SN REAL ZERO, ONE PARAMETER (ONE = 1.0E0, ZERO = 0.0E0) IF(CS .GE. ZERO .AND. CS .GE. ABS(SN)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = SN/(ONE+CS) DO I=1,N WR=REAL(X(I)) WC=AIMAG(X(I)) UR = -SN*(WR*T+WC)+WR UC = SN*(-WC*T+WR)+WC VR=REAL(Y(I)) VC=AIMAG(Y(I)) 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 T = SN/(ONE+CS) DO I=1,N WR=REAL(X((I-1)*INCX+1)) WC=AIMAG(X((I-1)*INCX+1)) UR = -SN*(WR*T+WC)+WR UC = SN*(-WC*T+WR)+WC VR=REAL(Y((I-1)*INCY+1)) VC=AIMAG(Y((I-1)*INCY+1)) 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 ELSE IF(SN .GE. ZERO .AND. SN .GE. ABS(CS)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = CS/(ONE+SN) DO I=1,N WR=REAL(X(I)) WC=AIMAG(X(I)) UR = CS*(WC*T+WR)-WC UC = CS*(-WR*T+WC)+WR VR=REAL(Y(I)) VC=AIMAG(Y(I)) 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 T = CS/(ONE+SN) DO I=1,N WR=REAL(X((I-1)*INCX+1)) WC=AIMAG(X((I-1)*INCX+1)) UR = CS*(WC*T+WR)-WC UC = CS*(-WR*T+WC)+WR VR=REAL(Y((I-1)*INCY+1)) VC=AIMAG(Y((I-1)*INCY+1)) 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 ELSE IF(CS .LE. ZERO .AND. ABS(CS) .GE. ABS(SN)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = SN/(ONE-CS) DO I=1,N WR=REAL(X(I)) WC=AIMAG(X(I)) UR = SN*(WR*T-WC)-WR UC = SN*(WC*T+WR)-WC VR=REAL(Y(I)) VC=AIMAG(Y(I)) 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 T = SN/(ONE-CS) DO I=1,N WR=REAL(X((I-1)*INCX+1)) WC=AIMAG(X((I-1)*INCX+1)) UR = SN*(WR*T-WC)-WR UC = SN*(WC*T+WR)-WC VR=REAL(Y((I-1)*INCY+1)) VC=AIMAG(Y((I-1)*INCY+1)) 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 ELSE IF(SN .LE. ZERO .AND. ABS(SN) .GE. ABS(CS)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = CS/(ONE-SN) DO I=1,N WR=REAL(X(I)) WC=AIMAG(X(I)) UR = CS*(-WC*T+WR)+WC UC = CS*(WR*T+WC)-WR VR=REAL(Y(I)) VC=AIMAG(Y(I)) 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 T = CS/(ONE-SN) DO I=1,N WR=REAL(X((I-1)*INCX+1)) WC=AIMAG(X((I-1)*INCX+1)) UR = CS*(-WC*T+WR)+WC UC = CS*(WR*T+WC)-WR VR=REAL(Y((I-1)*INCY+1)) VC=AIMAG(Y((I-1)*INCY+1)) 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 ENDIF RETURN END SUBROUTINE CSROTUX SUBROUTINE CSROTUY(N, X, INCX, Y, INCY, F, S, CS, SN) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX X(*),Y(*) REAL F,S,T,WR,WC,VR,VC,UR,UC REAL CS,SN REAL ZERO, ONE PARAMETER (ONE = 1.0E0, ZERO = 0.0E0) IF(CS .GE. ZERO .AND. CS .GE. ABS(SN)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = SN/(ONE+CS) DO I=1,N UR=REAL(X(I)) UC=AIMAG(X(I)) WR=REAL(Y(I)) WC=AIMAG(Y(I)) VR = -SN*(WR*T+WC)+WR VC = SN*(-WC*T+WR)+WC 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 T = SN/(ONE+CS) DO I=1,N UR=REAL(X((I-1)*INCX+1)) UC=AIMAG(X((I-1)*INCX+1)) WR=REAL(Y((I-1)*INCY+1)) WC=AIMAG(Y((I-1)*INCY+1)) VR = -SN*(WR*T+WC)+WR VC = SN*(-WC*T+WR)+WC 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 ELSE IF(SN .GE. ZERO .AND. SN .GE. ABS(CS)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = CS/(ONE+SN) DO I=1,N UR=REAL(X(I)) UC=AIMAG(X(I)) WR=REAL(Y(I)) WC=AIMAG(Y(I)) VR = CS*(WC*T+WR)-WC VC = CS*(-WR*T+WC)+WR 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 T = CS/(ONE+SN) DO I=1,N UR=REAL(X((I-1)*INCX+1)) UC=AIMAG(X((I-1)*INCX+1)) WR=REAL(Y((I-1)*INCY+1)) WC=AIMAG(Y((I-1)*INCY+1)) VR = CS*(WC*T+WR)-WC VC = CS*(-WR*T+WC)+WR 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 ELSE IF(CS .LE. ZERO .AND. ABS(CS) .GE. ABS(SN)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = SN/(ONE-CS) DO I=1,N UR=REAL(X(I)) UC=AIMAG(X(I)) WR=REAL(Y(I)) WC=AIMAG(Y(I)) VR = SN*(WR*T-WC)-WR VC = SN*(WC*T+WR)-WC 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 T = SN/(ONE-CS) DO I=1,N UR=REAL(X((I-1)*INCX+1)) UC=AIMAG(X((I-1)*INCX+1)) WR=REAL(Y((I-1)*INCY+1)) WC=AIMAG(Y((I-1)*INCY+1)) VR = SN*(WR*T-WC)-WR VC = SN*(WC*T+WR)-WC 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 ELSE IF(SN .LE. ZERO .AND. ABS(SN) .GE. ABS(CS)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = CS/(ONE-SN) DO I=1,N UR=REAL(X(I)) UC=AIMAG(X(I)) WR=REAL(Y(I)) WC=AIMAG(Y(I)) VR = CS*(-WC*T+WR)+WC VC = CS*(WR*T+WC)-WR 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 T = CS/(ONE-SN) DO I=1,N UR=REAL(X((I-1)*INCX+1)) UC=AIMAG(X((I-1)*INCX+1)) WR=REAL(Y((I-1)*INCY+1)) WC=AIMAG(Y((I-1)*INCY+1)) VR = CS*(-WC*T+WR)+WC VC = CS*(WR*T+WC)-WR 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 ENDIF RETURN END SUBROUTINE CSROTUY SUBROUTINE CSROTXX(N, X, INCX, Y, INCY, F, S, CS, SN) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX X(*),Y(*) REAL F,S,T,WR,WC,VR,VC,UR,UC REAL CS, SN REAL ZERO, ONE PARAMETER (ONE = 1.0E0, ZERO = 0.0E0) IF(CS .GE. ZERO .AND. CS .GE. ABS(SN)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = SN/(ONE+CS) DO I=1,N WR=REAL(X(I)) WC=AIMAG(X(I)) UR = -SN*(WR*T+WC)+WR UC = SN*(-WC*T+WR)+WC VR=REAL(Y(I)) VC=AIMAG(Y(I)) 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 T = SN/(ONE+CS) DO I=1,N WR=REAL(X((I-1)*INCX+1)) WC=AIMAG(X((I-1)*INCX+1)) UR = -SN*(WR*T+WC)+WR UC = SN*(-WC*T+WR)+WC VR=REAL(Y((I-1)*INCY+1)) VC=AIMAG(Y((I-1)*INCY+1)) 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 ELSE IF(SN .GE. ZERO .AND. SN .GE. ABS(CS)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = CS/(ONE+SN) DO I=1,N WR=REAL(X(I)) WC=AIMAG(X(I)) UR = CS*(WC*T+WR)-WC UC = CS*(-WR*T+WC)+WR VR=REAL(Y(I)) VC=AIMAG(Y(I)) 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 T = CS/(ONE+SN) DO I=1,N WR=REAL(X((I-1)*INCX+1)) WC=AIMAG(X((I-1)*INCX+1)) UR = CS*(WC*T+WR)-WC UC = CS*(-WR*T+WC)+WR VR=REAL(Y((I-1)*INCY+1)) VC=AIMAG(Y((I-1)*INCY+1)) 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 ELSE IF(CS .LE. ZERO .AND. ABS(CS) .GE. ABS(SN)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = SN/(ONE-CS) DO I=1,N WR=REAL(X(I)) WC=AIMAG(X(I)) UR = SN*(WR*T-WC)-WR UC = SN*(WC*T+WR)-WC VR=REAL(Y(I)) VC=AIMAG(Y(I)) 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 T = SN/(ONE-CS) DO I=1,N WR=REAL(X((I-1)*INCX+1)) WC=AIMAG(X((I-1)*INCX+1)) UR = SN*(WR*T-WC)-WR UC = SN*(WC*T+WR)-WC VR=REAL(Y((I-1)*INCY+1)) VC=AIMAG(Y((I-1)*INCY+1)) 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 ELSE IF(SN .LE. ZERO .AND. ABS(SN) .GE. ABS(CS)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = CS/(ONE-SN) DO I=1,N WR=REAL(X(I)) WC=AIMAG(X(I)) UR = CS*(-WC*T+WR)+WC UC = CS*(WR*T+WC)-WR VR=REAL(Y(I)) VC=AIMAG(Y(I)) 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 T = CS/(ONE-SN) DO I=1,N WR=REAL(X((I-1)*INCX+1)) WC=AIMAG(X((I-1)*INCX+1)) UR = CS*(-WC*T+WR)+WC UC = CS*(WR*T+WC)-WR VR=REAL(Y((I-1)*INCY+1)) VC=AIMAG(Y((I-1)*INCY+1)) 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 ENDIF RETURN END SUBROUTINE CSROTXX SUBROUTINE CSROTXY(N, X, INCX, Y, INCY, F, S, CS, SN) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX X(*),Y(*) REAL F,S,T,WR,WC,VR,VC,UR,UC REAL CS, SN REAL ZERO, ONE PARAMETER (ONE = 1.0E0, ZERO = 0.0E0) IF(CS .GE. ZERO .AND. CS .GE. ABS(SN)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = SN/(ONE+CS) DO I=1,N UR=REAL(X(I)) UC=AIMAG(X(I)) WR=REAL(Y(I)) WC=AIMAG(Y(I)) VR = -SN*(WR*T+WC)+WR VC = SN*(-WC*T+WR)+WC 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 T = SN/(ONE+CS) DO I=1,N UR=REAL(X((I-1)*INCX+1)) UC=AIMAG(X((I-1)*INCX+1)) WR=REAL(Y((I-1)*INCY+1)) WC=AIMAG(Y((I-1)*INCY+1)) VR = -SN*(WR*T+WC)+WR VC = SN*(-WC*T+WR)+WC 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 ELSE IF(SN .GE. ZERO .AND. SN .GE. ABS(CS)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = CS/(ONE+SN) DO I=1,N UR=REAL(X(I)) UC=AIMAG(X(I)) WR=REAL(Y(I)) WC=AIMAG(Y(I)) VR = CS*(WC*T+WR)-WC VC = CS*(-WR*T+WC)+WR 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 T = CS/(ONE+SN) DO I=1,N UR=REAL(X((I-1)*INCX+1)) UC=AIMAG(X((I-1)*INCX+1)) WR=REAL(Y((I-1)*INCY+1)) WC=AIMAG(Y((I-1)*INCY+1)) VR = CS*(WC*T+WR)-WC VC = CS*(-WR*T+WC)+WR 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 ELSE IF(CS .LE. ZERO .AND. ABS(CS) .GE. ABS(SN)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = SN/(ONE-CS) DO I=1,N UR=REAL(X(I)) UC=AIMAG(X(I)) WR=REAL(Y(I)) WC=AIMAG(Y(I)) VR = SN*(WR*T-WC)-WR VC = SN*(WC*T+WR)-WC 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 T = SN/(ONE-CS) DO I=1,N UR=REAL(X((I-1)*INCX+1)) UC=AIMAG(X((I-1)*INCX+1)) WR=REAL(Y((I-1)*INCY+1)) WC=AIMAG(Y((I-1)*INCY+1)) VR = SN*(WR*T-WC)-WR VC = SN*(WC*T+WR)-WC 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 ELSE IF(SN .LE. ZERO .AND. ABS(SN) .GE. ABS(CS)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = CS/(ONE-SN) DO I=1,N UR=REAL(X(I)) UC=AIMAG(X(I)) WR=REAL(Y(I)) WC=AIMAG(Y(I)) VR = CS*(-WC*T+WR)+WC VC = CS*(WR*T+WC)-WR 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 T = CS/(ONE-SN) DO I=1,N UR=REAL(X((I-1)*INCX+1)) UC=AIMAG(X((I-1)*INCX+1)) WR=REAL(Y((I-1)*INCY+1)) WC=AIMAG(Y((I-1)*INCY+1)) VR = CS*(-WC*T+WR)+WC VC = CS*(WR*T+WC)-WR 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 ENDIF RETURN END SUBROUTINE CSROTXY SUBROUTINE CSROTVX(N, X, INCX, Y, INCY, C, G, CS, SN) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX X(*),Y(*) REAL G,C,T,WR,WC,UR,UC,VR,VC REAL CS, SN REAL ZERO, ONE PARAMETER (ONE = 1.0E0, ZERO = 0.0E0) IF(CS .GE. ZERO .AND. CS .GE. ABS(SN)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = SN/(ONE+CS) DO I=1,N WR=REAL(X(I)) WC=AIMAG(X(I)) UR = -SN*(WR*T+WC)+WR UC = SN*(-WC*T+WR)+WC VR=REAL(Y(I)) VC=AIMAG(Y(I)) 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 T = SN/(ONE+CS) DO I=1,N WR=REAL(X((I-1)*INCX+1)) WC=AIMAG(X((I-1)*INCX+1)) UR = -SN*(WR*T+WC)+WR UC = SN*(-WC*T+WR)+WC VR=REAL(Y((I-1)*INCY+1)) VC=AIMAG(Y((I-1)*INCY+1)) 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 ELSE IF(SN .GE. ZERO .AND. SN .GE. ABS(CS)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = CS/(ONE+SN) DO I=1,N WR=REAL(X(I)) WC=AIMAG(X(I)) UR = CS*(WC*T+WR)-WC UC = CS*(-WR*T+WC)+WR VR=REAL(Y(I)) VC=AIMAG(Y(I)) 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 T = CS/(ONE+SN) DO I=1,N WR=REAL(X((I-1)*INCX+1)) WC=AIMAG(X((I-1)*INCX+1)) UR = CS*(WC*T+WR)-WC UC = CS*(-WR*T+WC)+WR VR=REAL(Y((I-1)*INCY+1)) VC=AIMAG(Y((I-1)*INCY+1)) 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 ELSE IF(CS .LE. ZERO .AND. ABS(CS) .GE. ABS(SN)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = SN/(ONE-CS) DO I=1,N WR=REAL(X(I)) WC=AIMAG(X(I)) UR = SN*(WR*T-WC)-WR UC = SN*(WC*T+WR)-WC VR=REAL(Y(I)) VC=AIMAG(Y(I)) 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 T = SN/(ONE-CS) DO I=1,N WR=REAL(X((I-1)*INCX+1)) WC=AIMAG(X((I-1)*INCX+1)) UR = SN*(WR*T-WC)-WR UC = SN*(WC*T+WR)-WC VR=REAL(Y((I-1)*INCY+1)) VC=AIMAG(Y((I-1)*INCY+1)) 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 ELSE IF(SN .LE. ZERO .AND. ABS(SN) .GE. ABS(CS)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = CS/(ONE-SN) DO I=1,N WR=REAL(X(I)) WC=AIMAG(X(I)) UR = CS*(-WC*T+WR)+WC UC = CS*(WR*T+WC)-WR VR=REAL(Y(I)) VC=AIMAG(Y(I)) 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 T = CS/(ONE-SN) DO I=1,N WR=REAL(X((I-1)*INCX+1)) WC=AIMAG(X((I-1)*INCX+1)) UR = CS*(-WC*T+WR)+WC UC = CS*(WR*T+WC)-WR VR=REAL(Y((I-1)*INCY+1)) VC=AIMAG(Y((I-1)*INCY+1)) 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 ENDIF RETURN END SUBROUTINE CSROTVX SUBROUTINE CSROTVY(N, X, INCX, Y, INCY, C, G, CS, SN) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX X(*),Y(*) REAL G,C,T,WR,WC,UR,UC,VR,VC REAL CS, SN REAL ZERO, ONE PARAMETER (ONE = 1.0E0, ZERO = 0.0E0) IF(CS .GE. ZERO .AND. CS .GE. ABS(SN)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = SN/(ONE+CS) DO I=1,N UR=REAL(X(I)) UC=AIMAG(X(I)) WR=REAL(Y(I)) WC=AIMAG(Y(I)) VR = -SN*(WR*T+WC)+WR VC = SN*(-WC*T+WR)+WC 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 T = SN/(ONE+CS) DO I=1,N UR=REAL(X((I-1)*INCX+1)) UC=AIMAG(X((I-1)*INCX+1)) WR=REAL(Y((I-1)*INCY+1)) WC=AIMAG(Y((I-1)*INCY+1)) VR = -SN*(WR*T+WC)+WR VC = SN*(-WC*T+WR)+WC 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 ELSE IF(SN .GE. ZERO .AND. SN .GE. ABS(CS)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = CS/(ONE+SN) DO I=1,N UR=REAL(X(I)) UC=AIMAG(X(I)) WR=REAL(Y(I)) WC=AIMAG(Y(I)) VR = CS*(WC*T+WR)-WC VC = CS*(-WR*T+WC)+WR 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 T = CS/(ONE+SN) DO I=1,N UR=REAL(X((I-1)*INCX+1)) UC=AIMAG(X((I-1)*INCX+1)) WR=REAL(Y((I-1)*INCY+1)) WC=AIMAG(Y((I-1)*INCY+1)) VR = CS*(WC*T+WR)-WC VC = CS*(-WR*T+WC)+WR 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 ELSE IF(CS .LE. ZERO .AND. ABS(CS) .GE. ABS(SN)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = SN/(ONE-CS) DO I=1,N UR=REAL(X(I)) UC=AIMAG(X(I)) WR=REAL(Y(I)) WC=AIMAG(Y(I)) VR = SN*(WR*T-WC)-WR VC = SN*(WC*T+WR)-WC 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 T = SN/(ONE-CS) DO I=1,N UR=REAL(X((I-1)*INCX+1)) UC=AIMAG(X((I-1)*INCX+1)) WR=REAL(Y((I-1)*INCY+1)) WC=AIMAG(Y((I-1)*INCY+1)) VR = SN*(WR*T-WC)-WR VC = SN*(WC*T+WR)-WC 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 ELSE IF(SN .LE. ZERO .AND. ABS(SN) .GE. ABS(CS)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = CS/(ONE-SN) DO I=1,N UR=REAL(X(I)) UC=AIMAG(X(I)) WR=REAL(Y(I)) WC=AIMAG(Y(I)) VR = CS*(-WC*T+WR)+WC VC = CS*(WR*T+WC)-WR 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 T = CS/(ONE-SN) DO I=1,N UR=REAL(X((I-1)*INCX+1)) UC=AIMAG(X((I-1)*INCX+1)) WR=REAL(Y((I-1)*INCY+1)) WC=AIMAG(Y((I-1)*INCY+1)) VR = CS*(-WC*T+WR)+WC VC = CS*(WR*T+WC)-WR 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 ENDIF RETURN END SUBROUTINE CSROTVY SUBROUTINE CSROTWX(N, X, INCX, Y, INCY, C, G, CS, SN) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX X(*),Y(*) REAL G,C,T,WR,WC,UR,UC,VR,VC REAL CS, SN REAL ZERO, ONE PARAMETER (ONE = 1.0E0, ZERO = 0.0E0) IF(CS .GE. ZERO .AND. CS .GE. ABS(SN)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = SN/(ONE+CS) DO I=1,N WR=REAL(X(I)) WC=AIMAG(X(I)) UR = -SN*(WR*T+WC)+WR UC = SN*(-WC*T+WR)+WC VR=REAL(Y(I)) VC=AIMAG(Y(I)) 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 T = SN/(ONE+CS) DO I=1,N WR=REAL(X((I-1)*INCX+1)) WC=AIMAG(X((I-1)*INCX+1)) UR = -SN*(WR*T+WC)+WR UC = SN*(-WC*T+WR)+WC VR=REAL(Y((I-1)*INCY+1)) VC=AIMAG(Y((I-1)*INCY+1)) 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 ELSE IF(SN .GE. ZERO .AND. SN .GE. ABS(CS)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = CS/(ONE+SN) DO I=1,N WR=REAL(X(I)) WC=AIMAG(X(I)) UR = CS*(WC*T+WR)-WC UC = CS*(-WR*T+WC)+WR VR=REAL(Y(I)) VC=AIMAG(Y(I)) 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 T = CS/(ONE+SN) DO I=1,N WR=REAL(X((I-1)*INCX+1)) WC=AIMAG(X((I-1)*INCX+1)) UR = CS*(WC*T+WR)-WC UC = CS*(-WR*T+WC)+WR VR=REAL(Y((I-1)*INCY+1)) VC=AIMAG(Y((I-1)*INCY+1)) 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 ELSE IF(CS .LE. ZERO .AND. ABS(CS) .GE. ABS(SN)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = SN/(ONE-CS) DO I=1,N WR=REAL(X(I)) WC=AIMAG(X(I)) UR = SN*(WR*T-WC)-WR UC = SN*(WC*T+WR)-WC VR=REAL(Y(I)) VC=AIMAG(Y(I)) 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 T = SN/(ONE-CS) DO I=1,N WR=REAL(X((I-1)*INCX+1)) WC=AIMAG(X((I-1)*INCX+1)) UR = SN*(WR*T-WC)-WR UC = SN*(WC*T+WR)-WC VR=REAL(Y((I-1)*INCY+1)) VC=AIMAG(Y((I-1)*INCY+1)) 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 ELSE IF(SN .LE. ZERO .AND. ABS(SN) .GE. ABS(CS)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = CS/(ONE-SN) DO I=1,N WR=REAL(X(I)) WC=AIMAG(X(I)) UR = CS*(-WC*T+WR)+WC UC = CS*(WR*T+WC)-WR VR=REAL(Y(I)) VC=AIMAG(Y(I)) 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 T = CS/(ONE-SN) DO I=1,N WR=REAL(X((I-1)*INCX+1)) WC=AIMAG(X((I-1)*INCX+1)) UR = CS*(-WC*T+WR)+WC UC = CS*(WR*T+WC)-WR VR=REAL(Y((I-1)*INCY+1)) VC=AIMAG(Y((I-1)*INCY+1)) 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 ENDIF RETURN END SUBROUTINE CSROTWX SUBROUTINE CSROTWY(N, X, INCX, Y, INCY, C, G, CS, SN) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX X(*),Y(*) REAL G,C,T,WR,WC,UR,UC,VR,VC REAL CS, SN REAL ZERO, ONE PARAMETER (ONE = 1.0E0, ZERO = 0.0E0) IF(CS .GE. ZERO .AND. CS .GE. ABS(SN)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = SN/(ONE+CS) DO I=1,N UR=REAL(X(I)) UC=AIMAG(X(I)) WR=REAL(Y(I)) WC=AIMAG(Y(I)) VR = -SN*(WR*T+WC)+WR VC = SN*(-WC*T+WR)+WC 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 T = SN/(ONE+CS) DO I=1,N UR=REAL(X((I-1)*INCX+1)) UC=AIMAG(X((I-1)*INCX+1)) WR=REAL(Y((I-1)*INCY+1)) WC=AIMAG(Y((I-1)*INCY+1)) VR = -SN*(WR*T+WC)+WR VC = SN*(-WC*T+WR)+WC 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 ELSE IF(SN .GE. ZERO .AND. SN .GE. ABS(CS)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = CS/(ONE+SN) DO I=1,N UR=REAL(X(I)) UC=AIMAG(X(I)) WR=REAL(Y(I)) WC=AIMAG(Y(I)) VR = CS*(WC*T+WR)-WC VC = CS*(-WR*T+WC)+WR 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 T = CS/(ONE+SN) DO I=1,N UR=REAL(X((I-1)*INCX+1)) UC=AIMAG(X((I-1)*INCX+1)) WR=REAL(Y((I-1)*INCY+1)) WC=AIMAG(Y((I-1)*INCY+1)) VR = CS*(WC*T+WR)-WC VC = CS*(-WR*T+WC)+WR 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 ELSE IF(CS .LE. ZERO .AND. ABS(CS) .GE. ABS(SN)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = SN/(ONE-CS) DO I=1,N UR=REAL(X(I)) UC=AIMAG(X(I)) WR=REAL(Y(I)) WC=AIMAG(Y(I)) VR = SN*(WR*T-WC)-WR VC = SN*(WC*T+WR)-WC 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 T = SN/(ONE-CS) DO I=1,N UR=REAL(X((I-1)*INCX+1)) UC=AIMAG(X((I-1)*INCX+1)) WR=REAL(Y((I-1)*INCY+1)) WC=AIMAG(Y((I-1)*INCY+1)) VR = SN*(WR*T-WC)-WR VC = SN*(WC*T+WR)-WC 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 ELSE IF(SN .LE. ZERO .AND. ABS(SN) .GE. ABS(CS)) THEN IF (INCX .EQ. 1 .AND. INCY .EQ. 1) THEN T = CS/(ONE-SN) DO I=1,N UR=REAL(X(I)) UC=AIMAG(X(I)) WR=REAL(Y(I)) WC=AIMAG(Y(I)) VR = CS*(-WC*T+WR)+WC VC = CS*(WR*T+WC)-WR 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 T = CS/(ONE-SN) DO I=1,N UR=REAL(X((I-1)*INCX+1)) UC=AIMAG(X((I-1)*INCX+1)) WR=REAL(Y((I-1)*INCY+1)) WC=AIMAG(Y((I-1)*INCY+1)) VR = CS*(-WC*T+WR)+WC VC = CS*(WR*T+WC)-WR 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 ENDIF RETURN END SUBROUTINE CSROTWY SUBROUTINE CSCALG(N, X, INCX, CS, SN) IMPLICIT NONE INTEGER N, I, INCX COMPLEX X(*) REAL CS, SN, T, XR, XC REAL ZERO, ONE PARAMETER (ZERO = 0.0E0, ONE = 1.0E0) IF(CS .GE. ZERO .AND. CS .GE. ABS(SN)) THEN IF (INCX .EQ. 1) THEN T = SN/(ONE+CS) DO I=1,N XR=REAL(X(I)) XC=AIMAG(X(I)) X(I) = CMPLX(-SN*(XR*T+XC)+XR, SN*(-XC*T+XR)+XC) ENDDO ELSE T = SN/(ONE+CS) DO I=1,N XR=REAL(X((I-1)*INCX+1)) XC=AIMAG(X((I-1)*INCX+1)) X((I-1)*INCX+1) = CMPLX(-SN*(XR*T+XC)+XR, SN*(-XC*T+XR)+XC) ENDDO ENDIF ELSE IF(SN .GE. ZERO .AND. SN .GE. ABS(CS)) THEN IF (INCX .EQ. 1) THEN T = CS/(ONE+SN) DO I=1,N XR=REAL(X(I)) XC=AIMAG(X(I)) X(I) = CMPLX(CS*(XC*T+XR)-XC, CS*(-XR*T+XC)+XR) ENDDO ELSE T = CS/(ONE+SN) DO I=1,N XR=REAL(X((I-1)*INCX+1)) XC=AIMAG(X((I-1)*INCX+1)) X((I-1)*INCX+1) = CMPLX(CS*(XC*T+XR)-XC, CS*(-XR*T+XC)+XR) ENDDO ENDIF ELSE IF(CS .LE. ZERO .AND. ABS(CS) .GE. ABS(SN)) THEN IF (INCX .EQ. 1) THEN T = SN/(ONE-CS) DO I=1,N XR=REAL(X(I)) XC=AIMAG(X(I)) X(I) = CMPLX(SN*(XR*T-XC)-XR, SN*(XC*T+XR)-XC) ENDDO ELSE T = SN/(ONE-CS) DO I=1,N XR=REAL(X((I-1)*INCX+1)) XC=AIMAG(X((I-1)*INCX+1)) X((I-1)*INCX+1) = CMPLX(SN*(XR*T-XC)-XR, SN*(XC*T+XR)-XC) ENDDO ENDIF ELSE IF(SN .LE. ZERO .AND. ABS(SN) .GE. ABS(CS)) THEN IF (INCX .EQ. 1) THEN T = CS/(ONE-SN) DO I=1,N XR=REAL(X(I)) XC=AIMAG(X(I)) X(I) = CMPLX(CS*(-XC*T+XR)+XC, CS*(XR*T+XC)-XR) ENDDO ELSE T = CS/(ONE-SN) DO I=1,N XR=REAL(X((I-1)*INCX+1)) XC=AIMAG(X((I-1)*INCX+1)) X((I-1)*INCX+1) = CMPLX(CS*(-XC*T+XR)+XC, CS*(XR*T+XC)-XR) ENDDO ENDIF ENDIF RETURN END SUBROUTINE CSCALG 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 CJASVD(N,A,SIGMA,U,V,INFO) use omp_lib IMPLICIT NONE INTEGER N, INFO, FLAG, LL, MM, FLAG2, FLAG3, FLAG4 REAL SIGMA(N),B(N),Z(N),TMP,TMP3,TMP4,TMP5 COMPLEX A(N,N),U(N,N),V(N,N) COMPLEX TMP2,TMP6 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, X, ALPHA, BETA REAL THETA1, THETA2, CS, SN, CSL, SNL, CSR, SNR REAL EPS, R1, R2 REAL T1, G1, F1, T2, G2, F2, C1, S1, C2, S2 EXTERNAL SFMA0, SLAMCH REAL SFMA0, SLAMCH EPS = SLAMCH('P') DO J=1,N TMP4 = REAL(A(J,J)) TMP5 = AIMAG(A(J,J)) CALL SLARTG(ABS(TMP4),ABS(TMP5),CSL,SNL,SIGMA(J)) IF (SIGMA(J) .GT. SZERO) THEN CSL = SIGN(CSL,TMP4) SNL = SIGN(SNL,TMP5) TMP2 = CMPLX(CSL,SNL) TMP2 = CONJG(TMP2) CALL CSCAL(J-1,TMP2,A(1,J),1) CALL CSCAL(N,TMP2,V(1,J),1) ENDIF A(J,J) = SIGMA(J) ENDDO IF (SIGMA(1) .GE. SIGMA(N)) THEN FLAG4 = 0 ELSE FLAG4 = 1 ENDIF LL = 1 MM = N DO ITER=1, MAXITER 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) = SZERO 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) = CZERO ELSE TMP4 = REAL(A(J,K)) TMP5 = AIMAG(A(J,K)) CALL SLARTG(ABS(TMP4),ABS(TMP5),CSL,SNL,TMP) CS = SIGN(CSL,TMP4) SN = SIGN(SNL,TMP5) F1 = SIGMA(J) - SIGMA(K) CALL SLARTGU (ABS(F1), ABS(TMP), C1, S1, R1) F1 = F1 + SIGN(R1,F1) IF (F1 .GE. SZERO) THEN G1 = TMP ELSE G1 = -TMP F1 = -F1 ENDIF F2 = SIGMA(J) + SIGMA(K) CALL SLARTGU (ABS(F2), ABS(TMP), C2, S2, R2) F2 = F2 + SIGN(R2,F2) IF (F2 .GE. SZERO) THEN G2 = -TMP ELSE G2 = TMP 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( TMP, R1, Z(J)) Z(K) = SFMA0( TMP, R2, Z(K)) SIGMA(J) = B(J)+Z(J) SIGMA(K) = B(K)+Z(K) IF (ABS(SIGMA(J)) .LT. ABS(SIGMA(K))) THEN THETA1 = CSL THETA2 = SNL IF (THETA2 .GT. SZERO) THEN FLAG2 = 0 CSL = THETA2 SNL = -THETA1 ELSE FLAG2 = 1 CSL = -THETA2 SNL = THETA1 ENDIF THETA1 = CSR THETA2 = SNR IF (THETA2 .GT. SZERO) THEN FLAG3 = 0 CSR = THETA2 SNR = -THETA1 ELSE FLAG3 = 1 CSR = -THETA2 SNR = THETA1 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/(SONE+CSL) CSL = SFMA0(-SNL,X,SONE) FLAG = 1 IF (CSL .NE. SONE .OR. SNL .NE. SZERO) THEN CALL CSROTUY (J-LL+1, A(LL, J), 1, A(LL, K), 1, X, SNL, CS, SN) CALL CSROTUY (MM-K+1, A(J, K), N, A(K, K), N, X, SNL, CS, SN) CALL CSROTUY (N, U(1, J), 1, U(1, K), 1, X, SNL, CS, -SN) ELSE CALL CSCALG (J-LL+1, A(LL, K), 1, CS, SN) CALL CSCALG (MM-K+1, A(K, K), N, CS, SN) CALL CSCALG (N, U(1, K), 1, CS, -SN) ENDIF ELSE IF(SNL .GE. SZERO) THEN X = CSL/(SONE+SNL) SNL = SFMA0(-CSL,X,SONE) FLAG = 1 CALL CSROTVY (J-LL+1, A(LL, J), 1, A(LL, K), 1, CSL, X, CS, SN) CALL CSROTVY (MM-K+1, A(J, K), N, A(K, K), N, CSL, X, CS, SN) CALL CSROTVY (N, U(1, J), 1, U(1, K), 1, CSL, X, CS, -SN) ELSE IF(SNL .LE. SZERO) THEN X = CSL/(SONE-SNL) SNL = SFMA0(CSL,X,-SONE) FLAG = 1 CALL CSROTWY (J-LL+1, A(LL, J), 1, A(LL, K), 1, CSL, X, CS, SN) CALL CSROTWY (MM-K+1, A(J, K), N, A(K, K), N, CSL, X, CS, SN) CALL CSROTWY (N, U(1, J), 1, U(1, K), 1, CSL, X, CS, -SN) ENDIF IF (CSR .GE. ABS(SNR)) THEN X = SNR/(SONE+CSR) CSR = SFMA0(-SNR,X,SONE) FLAG = 1 IF (CSR .NE. SONE .OR. SNR .NE. SZERO) THEN CALL CSROTUY (K-J+1, A(J,J), N, A(J,K), 1, X, SNR, CS, -SN) CALL CSROTUY (N, V(1, J), 1, V(1, K), 1, X, SNR, CS, -SN) ELSE CALL CSCALG (K-J+1, A(J, K), 1, CS, -SN) CALL CSCALG (N, V(1, K), 1, CS, -SN) ENDIF ELSE IF(SNR .GE. SZERO) THEN X = CSR/(SONE+SNR) SNR = SFMA0(-CSR,X,SONE) FLAG = 1 CALL CSROTVY (K-J+1, A(J,J), N, A(J,K), 1, CSR, X, CS, -SN) CALL CSROTVY (N, V(1, J), 1, V(1, K), 1, CSR, X, CS, -SN) ELSE IF(SNR .LE. SZERO) THEN X = CSR/(SONE-SNR) SNR = SFMA0(CSR,X,-SONE) FLAG = 1 CALL CSROTWY (K-J+1, A(J,J), N, A(J,K), 1, CSR, X, CS, -SN) CALL CSROTWY (N, V(1, J), 1, V(1, K), 1, CSR, X, CS, -SN) ENDIF A(J, K) = CZERO 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 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) = SZERO 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) = CZERO ELSE TMP4 = REAL(A(J,K)) TMP5 = AIMAG(A(J,K)) CALL SLARTG(ABS(TMP4),ABS(TMP5),CSL,SNL,TMP) CS = SIGN(CSL,TMP4) SN = SIGN(SNL,TMP5) F1 = SIGMA(J) - SIGMA(K) CALL SLARTGU (ABS(F1), ABS(TMP), C1, S1, R1) F1 = F1 + SIGN(R1,F1) IF (F1 .GE. SZERO) THEN G1 = TMP ELSE G1 = -TMP F1 = -F1 ENDIF F2 = SIGMA(J) + SIGMA(K) CALL SLARTGU (ABS(F2), ABS(TMP), C2, S2, R2) F2 = F2 + SIGN(R2,F2) IF (F2 .GE. SZERO) THEN G2 = TMP ELSE G2 = -TMP 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( TMP, R1, Z(J)) Z(K) = SFMA0( TMP, R2, Z(K)) SIGMA(J) = B(J)+Z(J) SIGMA(K) = B(K)+Z(K) IF (ABS(SIGMA(J)) .LT. ABS(SIGMA(K))) THEN THETA1 = CSL THETA2 = SNL IF (THETA2 .GT. SZERO) THEN FLAG2 = 0 CSL = THETA2 SNL = -THETA1 ELSE FLAG2 = 1 CSL = -THETA2 SNL = THETA1 ENDIF THETA1 = CSR THETA2 = SNR IF (THETA2 .GT. SZERO) THEN FLAG3 = 0 CSR = THETA2 SNR = -THETA1 ELSE FLAG3 = 1 CSR = -THETA2 SNR = THETA1 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/(SONE+CSR) CSR = SFMA0(-SNR,X,SONE) FLAG = 1 IF (CSR .NE. SONE .OR. SNR .NE. SZERO) THEN CALL CSROTUY (J-LL+1, A(LL,J), 1, A(LL,K), 1, X, SNR, CS, SN) CALL CSROTUY (MM-k+1, A(J,K), N, A(K,K), N, X, SNR, CS, SN) CALL CSROTUY (N, V(1, J), 1, V(1, K), 1, X, SNR, CS, SN) ELSE CALL CSCALG (J-LL+1, A(LL, K), 1, CS, SN) CALL CSCALG (MM-K+1, A(K, K), N, CS, SN) CALL CSCALG (N, V(1, K), 1, CS, SN) ENDIF ELSE IF(SNR .GE. SZERO) THEN X = CSR/(SONE+SNR) SNR = SFMA0(-CSR,X,SONE) FLAG = 1 CALL CSROTVY (J-LL+1, A(LL,J), 1, A(LL,K), 1, CSR, X, CS, SN) CALL CSROTVY (MM-k+1, A(J,K), N, A(K,K), N, CSR, X, CS, SN) CALL CSROTVY (N, V(1, J), 1, V(1, K), 1, CSR, X, CS, SN) ELSE IF(SNR .LE. SZERO) THEN X = CSR/(SONE-SNR) SNR = SFMA0(CSR,X,-SONE) FLAG = 1 CALL CSROTWY (J-LL+1, A(LL,J), 1, A(LL,K), 1, CSR, X, CS, SN) CALL CSROTWY (MM-K+1, A(J,K), N, A(K,K), N, CSR, X, CS, SN) CALL CSROTWY (N, V(1, J), 1, V(1, K), 1, CSR, X, CS, SN) ENDIF IF (CSL .GE. ABS(SNL)) THEN X = SNL/(SONE+CSL) CSL = SFMA0(-SNL,X,SONE) FLAG = 1 IF (CSL .NE. SONE .OR. SNL .NE. SZERO) THEN CALL CSROTUY (K-J+1, A(J, J), N, A(J, K), 1, X, SNL, CS, -SN) CALL CSROTUY (N, U(1, J), 1, U(1, K), 1, X, SNL, CS, SN) ELSE CALL CSCALG (K-J+1, A(J, K), 1, CS, -SN) CALL CSCALG (N, U(1, K), 1, CS, SN) ENDIF ELSE IF(SNL .GE. SZERO) THEN X = CSL/(SONE+SNL) SNL = SFMA0(-CSL,X,SONE) FLAG = 1 CALL CSROTVY (K-J+1, A(j, J), N, A(J, K), 1, CSL, X, CS, -SN) CALL CSROTVY (N, U(1, J), 1, U(1, K), 1, CSL, X, CS, SN) ELSE IF(SNL .LE. SZERO) THEN X = CSL/(SONE-SNL) SNL = SFMA0(CSL,X,-SONE) FLAG = 1 CALL CSROTWY (K-J+1, A(J, J), N, A(J, K), 1, CSL, X, CS, -SN) CALL CSROTWY (N, U(1, J), 1, U(1, K), 1, CSL, X, CS, SN) ENDIF A(J, K) = CZERO 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) = SZERO 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) = CZERO ELSE TMP4 = REAL(A(J,K)) TMP5 = AIMAG(A(J,K)) CALL SLARTG(ABS(TMP4),ABS(TMP5),CSL,SNL,TMP) CS = SIGN(CSL,TMP4) SN = SIGN(SNL,TMP5) F1 = SIGMA(J) - SIGMA(K) CALL SLARTGU (ABS(F1), ABS(TMP), C1, S1, R1) F1 = F1 + SIGN(R1,F1) IF (F1 .GE. SZERO) THEN G1 = TMP ELSE G1 = -TMP F1 = -F1 ENDIF F2 = SIGMA(J) + SIGMA(K) CALL SLARTGU (ABS(F2), ABS(TMP), C2, S2, R2) F2 = F2 + SIGN(R2,F2) IF (F2 .GE. SZERO) THEN G2 = -TMP ELSE G2 = TMP 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( TMP, R1, Z(J)) Z(K) = SFMA0( TMP, R2, Z(K)) SIGMA(J) = B(J)+Z(J) SIGMA(K) = B(K)+Z(K) IF (ABS(SIGMA(J)) .GE. ABS(SIGMA(K))) THEN THETA1 = CSL THETA2 = SNL IF (THETA2 .GT. SZERO) THEN FLAG2 = 0 CSL = THETA2 SNL = -THETA1 ELSE FLAG2 = 1 CSL = -THETA2 SNL = THETA1 ENDIF THETA1 = CSR THETA2 = SNR IF (THETA2 .GT. SZERO) THEN FLAG3 = 0 CSR = THETA2 SNR = -THETA1 ELSE FLAG3 = 1 CSR = -THETA2 SNR = THETA1 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/(SONE+CSL) CSL = SFMA0(-SNL,X,SONE) FLAG = 1 IF (CSL .NE. SONE .OR. SNL .NE. SZERO) THEN CALL CSROTUX (K-J+1, A(J, J), N, A(J, K), 1, X, SNL, CS, -SN) CALL CSROTUX (N, U(1, J), 1, U(1, K), 1, X, SNL, CS, SN) ELSE CALL CSCALG (K-J+1, A(J, J), N, CS, -SN) CALL CSCALG (N, U(1, J), 1, CS, SN) ENDIF ELSE IF(SNL .GE. SZERO) THEN X = CSL/(SONE+SNL) SNL = SFMA0(-CSL,X,SONE) FLAG = 1 CALL CSROTVX (K-J+1, A(J, J), N, A(J, K), 1, CSL, X, CS, -SN) CALL CSROTVX (N, U(1, J), 1, U(1, K), 1, CSL, X, CS, SN) ELSE IF(SNL .LE. SZERO) THEN X = CSL/(SONE-SNL) SNL = SFMA0(CSL,X,-SONE) FLAG = 1 CALL CSROTWX (K-J+1, A(J, J), N, A(J, K), 1, CSL, X, CS, -SN) CALL CSROTWX (N, U(1, J), 1, U(1, K), 1, CSL, X, CS, SN) ENDIF IF (CSR .GE. ABS(SNR)) THEN X = SNR/(SONE+CSR) CSR = SFMA0(-SNR,X,SONE) FLAG = 1 IF (CSR .NE. SONE .OR. SNR .NE. SZERO) THEN CALL CSROTUX (J-LL+1, A(LL,J), 1, A(LL,K), 1, X, SNR, CS, SN) CALL CSROTUX (MM-K+1, A(J,K), N, A(K,K), N, X, SNR, CS, SN) CALL CSROTUX (N, V(1, J), 1, V(1, K), 1, X, SNR, CS, SN) ELSE CALL CSCALG (J-LL+1, A(LL, J), 1, CS, SN) CALL CSCALG (MM-K+1, A(J, K), N, CS, SN) CALL CSCALG (N, V(1, J), 1, CS, SN) ENDIF ELSE IF(SNR .GE. SZERO) THEN X = CSR/(SONE+SNR) SNR = SFMA0(-CSR,X,SONE) FLAG = 1 CALL CSROTVX (J-LL+1, A(LL,J), 1, A(LL,K), 1, CSR, X, CS, SN) CALL CSROTVX (MM-K+1, A(J,K), N, A(K,K), N, CSR, X, CS, SN) CALL CSROTVX (N, V(1, J), 1, V(1, K), 1, CSR, X, CS, SN) ELSE IF(SNR .LE. SZERO) THEN X = CSR/(SONE-SNR) SNR = SFMA0(CSR,X,-SONE) FLAG = 1 CALL CSROTWX (J-LL+1, A(LL,J), 1, A(LL,K), 1, CSR, X, CS, SN) CALL CSROTWX (MM-K+1, A(J,K), N, A(K,K), N, CSR, X, CS, SN) CALL CSROTWX (N, V(1, J), 1, V(1, K), 1, CSR, X, CS, SN) ENDIF A(J, K) = CZERO 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 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) = SZERO 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) = CZERO ELSE TMP4 = REAL(A(J,K)) TMP5 = AIMAG(A(J,K)) CALL SLARTG(ABS(TMP4),ABS(TMP5),CSL,SNL,TMP) CS = SIGN(CSL,TMP4) SN = SIGN(SNL,TMP5) F1 = SIGMA(J) - SIGMA(K) CALL SLARTGU (ABS(F1), ABS(TMP), C1, S1, R1) F1 = F1 + SIGN(R1,F1) IF (F1 .GE. SZERO) THEN G1 = TMP ELSE G1 = -TMP F1 = -F1 ENDIF F2 = SIGMA(J) + SIGMA(K) CALL SLARTGU (ABS(F2), ABS(TMP), C2, S2, R2) F2 = F2 + SIGN(R2,F2) IF (F2 .GE. SZERO) THEN G2 = TMP ELSE G2 = -TMP 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( TMP, R1, Z(J)) Z(K) = SFMA0( TMP, R2, Z(K)) SIGMA(J) = B(J)+Z(J) SIGMA(K) = B(K)+Z(K) IF (ABS(SIGMA(J)) .GE. ABS(SIGMA(K))) THEN THETA1 = CSL THETA2 = SNL IF (THETA2 .GT. SZERO) THEN FLAG2 = 0 CSL = THETA2 SNL = -THETA1 ELSE FLAG2 = 1 CSL = -THETA2 SNL = THETA1 ENDIF THETA1 = CSR THETA2 = SNR IF (THETA2 .GT. SZERO) THEN FLAG3 = 0 CSR = THETA2 SNR = -THETA1 ELSE FLAG3 = 1 CSR = -THETA2 SNR = THETA1 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/(SONE+CSR) CSR = SFMA0(-SNR,X,SONE) FLAG = 1 IF (CSR .NE. SONE .OR. SNR .NE. SZERO) THEN CALL CSROTUX (K-J+1, A(J,J), N, A(J,K), 1, X, SNR, CS, -SN) CALL CSROTUX (N, V(1, J), 1, V(1, K), 1, X, SNR, CS, -SN) ELSE CALL CSCALG (K-J+1, A(J, J), N, CS, -SN) CALL CSCALG (N, V(1, J), 1, CS, -SN) ENDIF ELSE IF(SNR .GE. SZERO) THEN X = CSR/(SONE+SNR) SNR = SFMA0(-CSR,X,SONE) FLAG = 1 CALL CSROTVX (K-J+1, A(J,J), N, A(J,K), 1, CSR, X, CS, -SN) CALL CSROTVX (N, V(1, J), 1, V(1, K), 1, CSR, X, CS, -SN) ELSE IF(SNR .LE. SZERO) THEN X = CSR/(SONE-SNR) SNR = SFMA0(CSR,X,-SONE) FLAG = 1 CALL CSROTWX (K-J+1, A(J,J), N, A(J,K), 1, CSR, X, CS, -SN) CALL CSROTWX (N, V(1, J), 1, V(1, K), 1, CSR, X, CS, -SN) ENDIF IF (CSL .GE. ABS(SNL)) THEN X = SNL/(SONE+CSL) CSL = SFMA0(-SNL,X,SONE) FLAG = 1 IF (CSL .NE. SONE .OR. SNL .NE. SZERO) THEN CALL CSROTUX (J-LL+1, A(LL, J), 1, A(LL, K), 1, X, SNL, CS, SN) CALL CSROTUX (MM-K+1, A(J, K), N, A(K, K), N, X, SNL, CS, SN) CALL CSROTUX (N, U(1, J), 1, U(1, K), 1, X, SNL, CS, -SN) ELSE CALL CSCALG (J-LL+1, A(LL, J), 1, CS, SN) CALL CSCALG (MM-K+1, A(J, K), N, CS, SN) CALL CSCALG (N, U(1, J), 1, CS, -SN) ENDIF ELSE IF(SNL .GE. SZERO) THEN X = CSL/(SONE+SNL) SNL = SFMA0(-CSL,X,SONE) FLAG = 1 CALL CSROTVX (J-LL+1, A(LL, J), 1, A(LL, K), 1, CSL, X, CS, SN) CALL CSROTVX (MM-K+1, A(J, K), N, A(K, K), N, CSL, X, CS, SN) CALL CSROTVX (N, U(1, J), 1, U(1, K), 1, CSL, X, CS, -SN) ELSE IF(SNL .LE. SZERO) THEN X = CSL/(SONE-SNL) SNL = SFMA0(CSL,X,-SONE) FLAG = 1 CALL CSROTWX (J-LL+1, A(LL, J), 1, A(LL, K), 1, CSL, X, CS, SN) CALL CSROTWX (MM-K+1, A(J, K), N, A(K, K), N, CSL, X, CS, SN) CALL CSROTWX (N, U(1, J), 1, U(1, K), 1, CSL, X, CS, -SN) ENDIF A(J, K) = CZERO ENDIF ENDDO ENDDO IF(FLAG .EQ. 0) EXIT ENDIF ENDDO IF (ITER .GT. MAXITER) THEN INFO = -1 RETURN ENDIF DO I = 1,N IF (SIGMA(I) .LE. SZERO) THEN SIGMA(I) = -SIGMA(I) CALL CSCAL(N, -CONE, 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 CSWAP(N, U(1,K), 1, U(1,I), 1) CALL CSWAP(N, V(1,K), 1, V(1,I), 1) ENDDO INFO = 0 RETURN END SUBROUTINE CJASVD PROGRAM TESTDJS use omp_lib IMPLICIT NONE INTEGER I, J, K, N, INFO, SEED REAL TMP1, TMP2, PI 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 :: SIGMA(:) COMPLEX,ALLOCATABLE :: A(:,:), B(:,:), U(:,:), V(:,:) EXTERNAL USERRAND DOUBLE PRECISION USERRAND EXTERNAL CDOTC COMPLEX CDOTC character*1000 ARG1,ARG2 integer t1, t2, t_rate, t_max, diff PI = ACOS(-SONE) 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 TMP1 = USERRAND() TMP2 = USERRAND() A(J,I) = CMPLX(TMP1,TMP2) B(J,I) = A(J,I) ENDDO DO J=I+1,N A(J,I) = CZERO B(J,I) = CZERO ENDDO ENDDO DO J=1,N DO K=1,N U(K,J) = CZERO ENDDO U(J,J) = CONE 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(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 = CZERO !$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)*CONJG(V(J,K)) ENDDO E = E+CONJG(S)*S ENDDO ENDDO !$OMP END PARALLEL DO WRITE (*,*) 'ERROR',SQRT(REAL(E)) CALL CHERK ('U', 'C', N, N, SONE, U, N, SZERO, A, N) DO I=1,N A(I,I) = A(I,I) - CONE ENDDO E = CZERO !$OMP PARALLEL DO REDUCTION(+:E) DO I=2,N E = E+CDOTC(I-1,A(1,I),1,A(1,I),1) ENDDO !$OMP END PARALLEL DO E = STWO*E DO I=1,N E = E+CONJG(A(I,I))*A(I,I) ENDDO WRITE (*,*) 'ERROR',SQRT(REAL(E)) CALL CHERK ('U', 'C', N, N, SONE, V, N, SZERO, A, N) DO I=1,N A(I,I) = A(I,I) - CONE ENDDO E = CZERO !$OMP PARALLEL DO REDUCTION(+:E) DO I=2,N E = E+CDOTC(I-1,A(1,I),1,A(1,I),1) ENDDO !$OMP END PARALLEL DO E = STWO*E DO I=1,N E = E+CONJG(A(I,I))*A(I,I) ENDDO WRITE (*,*) 'ERROR',SQRT(REAL(E)) END PROGRAM TESTDJS