SUBROUTINE ZCROTUx(N, X, INCX, Y, INCY, F, S, CS, SN) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX*16 X(*),Y(*) DOUBLE PRECISION F,S,T,WR,WC,VR,VC,UR,UC,CS,SN DOUBLE PRECISION ZERO, ONE PARAMETER (ONE = 1.0D0, ZERO = 0.0D0) 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=DBLE(X(I)) WC=DIMAG(X(I)) UR = -SN*(WR*T+WC)+WR UC = SN*(-WC*T+WR)+WC VR=DBLE(Y(I)) VC=AIMAG(Y(I)) X(I)=DCMPLX(S*(-F*UR+VR)+UR,S*(-F*UC+VC)+UC) Y(I)=DCMPLX(-S*(F*VR+UR)+VR,-S*(F*VC+UC)+VC) ENDDO ELSE T = SN/(ONE+CS) DO I=1,N WR=DBLE(X((I-1)*INCX+1)) WC=DIMAG(X((I-1)*INCX+1)) UR = -SN*(WR*T+WC)+WR UC = SN*(-WC*T+WR)+WC VR=DBLE(Y((I-1)*INCY+1)) VC=DIMAG(Y((I-1)*INCY+1)) X((I-1)*INCX+1)=DCMPLX(S*(-F*UR+VR)+UR,S*(-F*UC+VC)+UC) Y((I-1)*INCY+1)=DCMPLX(-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=DBLE(X(I)) WC=DIMAG(X(I)) UR = CS*(WC*T+WR)-WC UC = CS*(-WR*T+WC)+WR VR=DBLE(Y(I)) VC=DIMAG(Y(I)) X(I)=DCMPLX(S*(-F*UR+VR)+UR,S*(-F*UC+VC)+UC) Y(I)=DCMPLX(-S*(F*VR+UR)+VR,-S*(F*VC+UC)+VC) ENDDO ELSE T = CS/(ONE+SN) DO I=1,N WR=DBLE(X((I-1)*INCX+1)) WC=DIMAG(X((I-1)*INCX+1)) UR = CS*(WC*T+WR)-WC UC = CS*(-WR*T+WC)+WR VR=DBLE(Y((I-1)*INCY+1)) VC=DIMAG(Y((I-1)*INCY+1)) X((I-1)*INCX+1)=DCMPLX(S*(-F*UR+VR)+UR,S*(-F*UC+VC)+UC) Y((I-1)*INCY+1)=DCMPLX(-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=DBLE(X(I)) WC=DIMAG(X(I)) UR = SN*(WR*T-WC)-WR UC = SN*(WC*T+WR)-WC VR=DBLE(Y(I)) VC=DIMAG(Y(I)) X(I)=DCMPLX(S*(-F*UR+VR)+UR,S*(-F*UC+VC)+UC) Y(I)=DCMPLX(-S*(F*VR+UR)+VR,-S*(F*VC+UC)+VC) ENDDO ELSE T = SN/(ONE-CS) DO I=1,N WR=DBLE(X((I-1)*INCX+1)) WC=DIMAG(X((I-1)*INCX+1)) UR = SN*(WR*T-WC)-WR UC = SN*(WC*T+WR)-WC VR=DBLE(Y((I-1)*INCY+1)) VC=DIMAG(Y((I-1)*INCY+1)) X((I-1)*INCX+1)=DCMPLX(S*(-F*UR+VR)+UR,S*(-F*UC+VC)+UC) Y((I-1)*INCY+1)=DCMPLX(-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=DBLE(X(I)) WC=DIMAG(X(I)) UR = CS*(-WC*T+WR)+WC UC = CS*(WR*T+WC)-WR VR=DBLE(Y(I)) VC=DIMAG(Y(I)) X(I)=DCMPLX(S*(-F*UR+VR)+UR,S*(-F*UC+VC)+UC) Y(I)=DCMPLX(-S*(F*VR+UR)+VR,-S*(F*VC+UC)+VC) ENDDO ELSE T = CS/(ONE-SN) DO I=1,N WR=DBLE(X((I-1)*INCX+1)) WC=DIMAG(X((I-1)*INCX+1)) UR = CS*(-WC*T+WR)+WC UC = CS*(WR*T+WC)-WR VR=DBLE(Y((I-1)*INCY+1)) VC=DIMAG(Y((I-1)*INCY+1)) X((I-1)*INCX+1)=DCMPLX(S*(-F*UR+VR)+UR,S*(-F*UC+VC)+UC) Y((I-1)*INCY+1)=DCMPLX(-S*(F*VR+UR)+VR,-S*(F*VC+UC)+VC) ENDDO ENDIF ENDIF RETURN END SUBROUTINE ZCROTUX SUBROUTINE ZCROTUY(N, X, INCX, Y, INCY, F, S, CS, SN) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX*16 X(*),Y(*) DOUBLE PRECISION F,S,T,WR,WC,VR,VC,UR,UC DOUBLE PRECISION CS,SN DOUBLE PRECISION ZERO, ONE PARAMETER (ONE = 1.0D0, ZERO = 0.0D0) 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=DBLE(X(I)) UC=DIMAG(X(I)) WR=DBLE(Y(I)) WC=DIMAG(Y(I)) VR = -SN*(WR*T+WC)+WR VC = SN*(-WC*T+WR)+WC X(I)=DCMPLX(S*(-F*UR+VR)+UR,S*(-F*UC+VC)+UC) Y(I)=DCMPLX(-S*(F*VR+UR)+VR,-S*(F*VC+UC)+VC) ENDDO ELSE T = SN/(ONE+CS) DO I=1,N UR=DBLE(X((I-1)*INCX+1)) UC=DIMAG(X((I-1)*INCX+1)) WR=DBLE(Y((I-1)*INCY+1)) WC=DIMAG(Y((I-1)*INCY+1)) VR = -SN*(WR*T+WC)+WR VC = SN*(-WC*T+WR)+WC X((I-1)*INCX+1)=DCMPLX(S*(-F*UR+VR)+UR,S*(-F*UC+VC)+UC) Y((I-1)*INCY+1)=DCMPLX(-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=DBLE(X(I)) UC=DIMAG(X(I)) WR=DBLE(Y(I)) WC=DIMAG(Y(I)) VR = CS*(WC*T+WR)-WC VC = CS*(-WR*T+WC)+WR X(I)=DCMPLX(S*(-F*UR+VR)+UR,S*(-F*UC+VC)+UC) Y(I)=DCMPLX(-S*(F*VR+UR)+VR,-S*(F*VC+UC)+VC) ENDDO ELSE T = CS/(ONE+SN) DO I=1,N UR=DBLE(X((I-1)*INCX+1)) UC=DIMAG(X((I-1)*INCX+1)) WR=DBLE(Y((I-1)*INCY+1)) WC=DIMAG(Y((I-1)*INCY+1)) VR = CS*(WC*T+WR)-WC VC = CS*(-WR*T+WC)+WR X((I-1)*INCX+1)=DCMPLX(S*(-F*UR+VR)+UR,S*(-F*UC+VC)+UC) Y((I-1)*INCY+1)=DCMPLX(-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=DBLE(X(I)) UC=DIMAG(X(I)) WR=DBLE(Y(I)) WC=DIMAG(Y(I)) VR = SN*(WR*T-WC)-WR VC = SN*(WC*T+WR)-WC X(I)=DCMPLX(S*(-F*UR+VR)+UR,S*(-F*UC+VC)+UC) Y(I)=DCMPLX(-S*(F*VR+UR)+VR,-S*(F*VC+UC)+VC) ENDDO ELSE T = SN/(ONE-CS) DO I=1,N UR=DBLE(X((I-1)*INCX+1)) UC=DIMAG(X((I-1)*INCX+1)) WR=DBLE(Y((I-1)*INCY+1)) WC=DIMAG(Y((I-1)*INCY+1)) VR = SN*(WR*T-WC)-WR VC = SN*(WC*T+WR)-WC X((I-1)*INCX+1)=DCMPLX(S*(-F*UR+VR)+UR,S*(-F*UC+VC)+UC) Y((I-1)*INCY+1)=DCMPLX(-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=DBLE(X(I)) UC=DIMAG(X(I)) WR=DBLE(Y(I)) WC=DIMAG(Y(I)) VR = CS*(-WC*T+WR)+WC VC = CS*(WR*T+WC)-WR X(I)=DCMPLX(S*(-F*UR+VR)+UR,S*(-F*UC+VC)+UC) Y(I)=DCMPLX(-S*(F*VR+UR)+VR,-S*(F*VC+UC)+VC) ENDDO ELSE T = CS/(ONE-SN) DO I=1,N UR=DBLE(X((I-1)*INCX+1)) UC=DIMAG(X((I-1)*INCX+1)) WR=DBLE(Y((I-1)*INCY+1)) WC=DIMAG(Y((I-1)*INCY+1)) VR = CS*(-WC*T+WR)+WC VC = CS*(WR*T+WC)-WR X((I-1)*INCX+1)=DCMPLX(S*(-F*UR+VR)+UR,S*(-F*UC+VC)+UC) Y((I-1)*INCY+1)=DCMPLX(-S*(F*VR+UR)+VR,-S*(F*VC+UC)+VC) ENDDO ENDIF ENDIF RETURN END SUBROUTINE ZCROTUY SUBROUTINE ZCROTXX(N, X, INCX, Y, INCY, F, S, CS, SN) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX*16 X(*),Y(*) DOUBLE PRECISION F,S,T,WR,WC,VR,VC,UR,UC DOUBLE PRECISION CS, SN DOUBLE PRECISION ZERO, ONE PARAMETER (ONE = 1.0D0, ZERO = 0.0D0) 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=DBLE(X(I)) WC=DIMAG(X(I)) UR = -SN*(WR*T+WC)+WR UC = SN*(-WC*T+WR)+WC VR=DBLE(Y(I)) VC=DIMAG(Y(I)) X(I)=DCMPLX(S*(F*UR+VR)-UR,S*(F*UC+VC)-UC) Y(I)=DCMPLX(S*(F*VR-UR)-VR,S*(F*VC-UC)-VC) ENDDO ELSE T = SN/(ONE+CS) DO I=1,N WR=DBLE(X((I-1)*INCX+1)) WC=DIMAG(X((I-1)*INCX+1)) UR = -SN*(WR*T+WC)+WR UC = SN*(-WC*T+WR)+WC VR=DBLE(Y((I-1)*INCY+1)) VC=DIMAG(Y((I-1)*INCY+1)) X((I-1)*INCX+1)=DCMPLX(S*(F*UR+VR)-UR,S*(F*UC+VC)-UC) Y((I-1)*INCY+1)=DCMPLX(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=DBLE(X(I)) WC=DIMAG(X(I)) UR = CS*(WC*T+WR)-WC UC = CS*(-WR*T+WC)+WR VR=DBLE(Y(I)) VC=DIMAG(Y(I)) X(I)=DCMPLX(S*(F*UR+VR)-UR,S*(F*UC+VC)-UC) Y(I)=DCMPLX(S*(F*VR-UR)-VR,S*(F*VC-UC)-VC) ENDDO ELSE T = CS/(ONE+SN) DO I=1,N WR=DBLE(X((I-1)*INCX+1)) WC=DIMAG(X((I-1)*INCX+1)) UR = CS*(WC*T+WR)-WC UC = CS*(-WR*T+WC)+WR VR=DBLE(Y((I-1)*INCY+1)) VC=DIMAG(Y((I-1)*INCY+1)) X((I-1)*INCX+1)=DCMPLX(S*(F*UR+VR)-UR,S*(F*UC+VC)-UC) Y((I-1)*INCY+1)=DCMPLX(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=DBLE(X(I)) WC=DIMAG(X(I)) UR = SN*(WR*T-WC)-WR UC = SN*(WC*T+WR)-WC VR=DBLE(Y(I)) VC=DIMAG(Y(I)) X(I)=DCMPLX(S*(F*UR+VR)-UR,S*(F*UC+VC)-UC) Y(I)=DCMPLX(S*(F*VR-UR)-VR,S*(F*VC-UC)-VC) ENDDO ELSE T = SN/(ONE-CS) DO I=1,N WR=DBLE(X((I-1)*INCX+1)) WC=DIMAG(X((I-1)*INCX+1)) UR = SN*(WR*T-WC)-WR UC = SN*(WC*T+WR)-WC VR=DBLE(Y((I-1)*INCY+1)) VC=DIMAG(Y((I-1)*INCY+1)) X((I-1)*INCX+1)=DCMPLX(S*(F*UR+VR)-UR,S*(F*UC+VC)-UC) Y((I-1)*INCY+1)=DCMPLX(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=DBLE(X(I)) WC=DIMAG(X(I)) UR = CS*(-WC*T+WR)+WC UC = CS*(WR*T+WC)-WR VR=DBLE(Y(I)) VC=DIMAG(Y(I)) X(I)=DCMPLX(S*(F*UR+VR)-UR,S*(F*UC+VC)-UC) Y(I)=DCMPLX(S*(F*VR-UR)-VR,S*(F*VC-UC)-VC) ENDDO ELSE T = CS/(ONE-SN) DO I=1,N WR=DBLE(X((I-1)*INCX+1)) WC=DIMAG(X((I-1)*INCX+1)) UR = CS*(-WC*T+WR)+WC UC = CS*(WR*T+WC)-WR VR=DBLE(Y((I-1)*INCY+1)) VC=DIMAG(Y((I-1)*INCY+1)) X((I-1)*INCX+1)=DCMPLX(S*(F*UR+VR)-UR,S*(F*UC+VC)-UC) Y((I-1)*INCY+1)=DCMPLX(S*(F*VR-UR)-VR,S*(F*VC-UC)-VC) ENDDO ENDIF ENDIF RETURN END SUBROUTINE ZCROTXX SUBROUTINE ZCROTXY(N, X, INCX, Y, INCY, F, S, CS, SN) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX*16 X(*),Y(*) DOUBLE PRECISION F,S,T,WR,WC,VR,VC,UR,UC DOUBLE PRECISION CS, SN DOUBLE PRECISION ZERO, ONE PARAMETER (ONE = 1.0D0, ZERO = 0.0D0) 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=DBLE(X(I)) UC=DIMAG(X(I)) WR=DBLE(Y(I)) WC=DIMAG(Y(I)) VR = -SN*(WR*T+WC)+WR VC = SN*(-WC*T+WR)+WC X(I)=DCMPLX(S*(F*UR+VR)-UR,S*(F*UC+VC)-UC) Y(I)=DCMPLX(S*(F*VR-UR)-VR,S*(F*VC-UC)-VC) ENDDO ELSE T = SN/(ONE+CS) DO I=1,N UR=DBLE(X((I-1)*INCX+1)) UC=DIMAG(X((I-1)*INCX+1)) WR=DBLE(Y((I-1)*INCY+1)) WC=DIMAG(Y((I-1)*INCY+1)) VR = -SN*(WR*T+WC)+WR VC = SN*(-WC*T+WR)+WC X((I-1)*INCX+1)=DCMPLX(S*(F*UR+VR)-UR,S*(F*UC+VC)-UC) Y((I-1)*INCY+1)=DCMPLX(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=DBLE(X(I)) UC=DIMAG(X(I)) WR=DBLE(Y(I)) WC=DIMAG(Y(I)) VR = CS*(WC*T+WR)-WC VC = CS*(-WR*T+WC)+WR X(I)=DCMPLX(S*(F*UR+VR)-UR,S*(F*UC+VC)-UC) Y(I)=DCMPLX(S*(F*VR-UR)-VR,S*(F*VC-UC)-VC) ENDDO ELSE T = CS/(ONE+SN) DO I=1,N UR=DBLE(X((I-1)*INCX+1)) UC=DIMAG(X((I-1)*INCX+1)) WR=DBLE(Y((I-1)*INCY+1)) WC=DIMAG(Y((I-1)*INCY+1)) VR = CS*(WC*T+WR)-WC VC = CS*(-WR*T+WC)+WR X((I-1)*INCX+1)=DCMPLX(S*(F*UR+VR)-UR,S*(F*UC+VC)-UC) Y((I-1)*INCY+1)=DCMPLX(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=DBLE(X(I)) UC=DIMAG(X(I)) WR=DBLE(Y(I)) WC=DIMAG(Y(I)) VR = SN*(WR*T-WC)-WR VC = SN*(WC*T+WR)-WC X(I)=DCMPLX(S*(F*UR+VR)-UR,S*(F*UC+VC)-UC) Y(I)=DCMPLX(S*(F*VR-UR)-VR,S*(F*VC-UC)-VC) ENDDO ELSE T = SN/(ONE-CS) DO I=1,N UR=DBLE(X((I-1)*INCX+1)) UC=DIMAG(X((I-1)*INCX+1)) WR=DBLE(Y((I-1)*INCY+1)) WC=DIMAG(Y((I-1)*INCY+1)) VR = SN*(WR*T-WC)-WR VC = SN*(WC*T+WR)-WC X((I-1)*INCX+1)=DCMPLX(S*(F*UR+VR)-UR,S*(F*UC+VC)-UC) Y((I-1)*INCY+1)=DCMPLX(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=DBLE(X(I)) UC=DIMAG(X(I)) WR=DBLE(Y(I)) WC=DIMAG(Y(I)) VR = CS*(-WC*T+WR)+WC VC = CS*(WR*T+WC)-WR X(I)=DCMPLX(S*(F*UR+VR)-UR,S*(F*UC+VC)-UC) Y(I)=DCMPLX(S*(F*VR-UR)-VR,S*(F*VC-UC)-VC) ENDDO ELSE T = CS/(ONE-SN) DO I=1,N UR=DBLE(X((I-1)*INCX+1)) UC=DIMAG(X((I-1)*INCX+1)) WR=DBLE(Y((I-1)*INCY+1)) WC=DIMAG(Y((I-1)*INCY+1)) VR = CS*(-WC*T+WR)+WC VC = CS*(WR*T+WC)-WR X((I-1)*INCX+1)=DCMPLX(S*(F*UR+VR)-UR,S*(F*UC+VC)-UC) Y((I-1)*INCY+1)=DCMPLX(S*(F*VR-UR)-VR,S*(F*VC-UC)-VC) ENDDO ENDIF ENDIF RETURN END SUBROUTINE ZCROTXY SUBROUTINE ZCROTVX(N, X, INCX, Y, INCY, C, G, CS, SN) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX*16 X(*),Y(*) DOUBLE PRECISION G,C,T,WR,WC,UR,UC,VR,VC DOUBLE PRECISION CS, SN DOUBLE PRECISION ZERO, ONE PARAMETER (ONE = 1.0D0, ZERO = 0.0D0) 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=DBLE(X(I)) WC=DIMAG(X(I)) UR = -SN*(WR*T+WC)+WR UC = SN*(-WC*T+WR)+WC VR=DBLE(Y(I)) VC=DIMAG(Y(I)) X(I)=DCMPLX(C*(-G*VR+UR)+VR,C*(-G*VC+UC)+VC) Y(I)=DCMPLX(C*(G*UR+VR)-UR,C*(G*UC+VC)-UC) ENDDO ELSE T = SN/(ONE+CS) DO I=1,N WR=DBLE(X((I-1)*INCX+1)) WC=DIMAG(X((I-1)*INCX+1)) UR = -SN*(WR*T+WC)+WR UC = SN*(-WC*T+WR)+WC VR=DBLE(Y((I-1)*INCY+1)) VC=DIMAG(Y((I-1)*INCY+1)) X((I-1)*INCX+1)=DCMPLX(C*(-G*VR+UR)+VR,C*(-G*VC+UC)+VC) Y((I-1)*INCY+1)=DCMPLX(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=DBLE(X(I)) WC=DIMAG(X(I)) UR = CS*(WC*T+WR)-WC UC = CS*(-WR*T+WC)+WR VR=DBLE(Y(I)) VC=DIMAG(Y(I)) X(I)=DCMPLX(C*(-G*VR+UR)+VR,C*(-G*VC+UC)+VC) Y(I)=DCMPLX(C*(G*UR+VR)-UR,C*(G*UC+VC)-UC) ENDDO ELSE T = CS/(ONE+SN) DO I=1,N WR=DBLE(X((I-1)*INCX+1)) WC=DIMAG(X((I-1)*INCX+1)) UR = CS*(WC*T+WR)-WC UC = CS*(-WR*T+WC)+WR VR=DBLE(Y((I-1)*INCY+1)) VC=DIMAG(Y((I-1)*INCY+1)) X((I-1)*INCX+1)=DCMPLX(C*(-G*VR+UR)+VR,C*(-G*VC+UC)+VC) Y((I-1)*INCY+1)=DCMPLX(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=DBLE(X(I)) WC=DIMAG(X(I)) UR = SN*(WR*T-WC)-WR UC = SN*(WC*T+WR)-WC VR=DBLE(Y(I)) VC=DIMAG(Y(I)) X(I)=DCMPLX(C*(-G*VR+UR)+VR,C*(-G*VC+UC)+VC) Y(I)=DCMPLX(C*(G*UR+VR)-UR,C*(G*UC+VC)-UC) ENDDO ELSE T = SN/(ONE-CS) DO I=1,N WR=DBLE(X((I-1)*INCX+1)) WC=DIMAG(X((I-1)*INCX+1)) UR = SN*(WR*T-WC)-WR UC = SN*(WC*T+WR)-WC VR=DBLE(Y((I-1)*INCY+1)) VC=DIMAG(Y((I-1)*INCY+1)) X((I-1)*INCX+1)=DCMPLX(C*(-G*VR+UR)+VR,C*(-G*VC+UC)+VC) Y((I-1)*INCY+1)=DCMPLX(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=DBLE(X(I)) WC=DIMAG(X(I)) UR = CS*(-WC*T+WR)+WC UC = CS*(WR*T+WC)-WR VR=DBLE(Y(I)) VC=DIMAG(Y(I)) X(I)=DCMPLX(C*(-G*VR+UR)+VR,C*(-G*VC+UC)+VC) Y(I)=DCMPLX(C*(G*UR+VR)-UR,C*(G*UC+VC)-UC) ENDDO ELSE T = CS/(ONE-SN) DO I=1,N WR=DBLE(X((I-1)*INCX+1)) WC=DIMAG(X((I-1)*INCX+1)) UR = CS*(-WC*T+WR)+WC UC = CS*(WR*T+WC)-WR VR=DBLE(Y((I-1)*INCY+1)) VC=DIMAG(Y((I-1)*INCY+1)) X((I-1)*INCX+1)=DCMPLX(C*(-G*VR+UR)+VR,C*(-G*VC+UC)+VC) Y((I-1)*INCY+1)=DCMPLX(C*(G*UR+VR)-UR,C*(G*UC+VC)-UC) ENDDO ENDIF ENDIF RETURN END SUBROUTINE ZCROTVX SUBROUTINE ZCROTVY(N, X, INCX, Y, INCY, C, G, CS, SN) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX*16 X(*),Y(*) DOUBLE PRECISION G,C,T,WR,WC,UR,UC,VR,VC DOUBLE PRECISION CS, SN DOUBLE PRECISION ZERO, ONE PARAMETER (ONE = 1.0D0, ZERO = 0.0D0) 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=DBLE(X(I)) UC=DIMAG(X(I)) WR=DBLE(Y(I)) WC=DIMAG(Y(I)) VR = -SN*(WR*T+WC)+WR VC = SN*(-WC*T+WR)+WC X(I)=DCMPLX(C*(-G*VR+UR)+VR,C*(-G*VC+UC)+VC) Y(I)=DCMPLX(C*(G*UR+VR)-UR,C*(G*UC+VC)-UC) ENDDO ELSE T = SN/(ONE+CS) DO I=1,N UR=DBLE(X((I-1)*INCX+1)) UC=DIMAG(X((I-1)*INCX+1)) WR=DBLE(Y((I-1)*INCY+1)) WC=DIMAG(Y((I-1)*INCY+1)) VR = -SN*(WR*T+WC)+WR VC = SN*(-WC*T+WR)+WC X((I-1)*INCX+1)=DCMPLX(C*(-G*VR+UR)+VR,C*(-G*VC+UC)+VC) Y((I-1)*INCY+1)=DCMPLX(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=DBLE(X(I)) UC=DIMAG(X(I)) WR=DBLE(Y(I)) WC=DIMAG(Y(I)) VR = CS*(WC*T+WR)-WC VC = CS*(-WR*T+WC)+WR X(I)=DCMPLX(C*(-G*VR+UR)+VR,C*(-G*VC+UC)+VC) Y(I)=DCMPLX(C*(G*UR+VR)-UR,C*(G*UC+VC)-UC) ENDDO ELSE T = CS/(ONE+SN) DO I=1,N UR=DBLE(X((I-1)*INCX+1)) UC=DIMAG(X((I-1)*INCX+1)) WR=DBLE(Y((I-1)*INCY+1)) WC=DIMAG(Y((I-1)*INCY+1)) VR = CS*(WC*T+WR)-WC VC = CS*(-WR*T+WC)+WR X((I-1)*INCX+1)=DCMPLX(C*(-G*VR+UR)+VR,C*(-G*VC+UC)+VC) Y((I-1)*INCY+1)=DCMPLX(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=DBLE(X(I)) UC=DIMAG(X(I)) WR=DBLE(Y(I)) WC=DIMAG(Y(I)) VR = SN*(WR*T-WC)-WR VC = SN*(WC*T+WR)-WC X(I)=DCMPLX(C*(-G*VR+UR)+VR,C*(-G*VC+UC)+VC) Y(I)=DCMPLX(C*(G*UR+VR)-UR,C*(G*UC+VC)-UC) ENDDO ELSE T = SN/(ONE-CS) DO I=1,N UR=DBLE(X((I-1)*INCX+1)) UC=DIMAG(X((I-1)*INCX+1)) WR=DBLE(Y((I-1)*INCY+1)) WC=DIMAG(Y((I-1)*INCY+1)) VR = SN*(WR*T-WC)-WR VC = SN*(WC*T+WR)-WC X((I-1)*INCX+1)=DCMPLX(C*(-G*VR+UR)+VR,C*(-G*VC+UC)+VC) Y((I-1)*INCY+1)=DCMPLX(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=DBLE(X(I)) UC=DIMAG(X(I)) WR=DBLE(Y(I)) WC=DIMAG(Y(I)) VR = CS*(-WC*T+WR)+WC VC = CS*(WR*T+WC)-WR X(I)=DCMPLX(C*(-G*VR+UR)+VR,C*(-G*VC+UC)+VC) Y(I)=DCMPLX(C*(G*UR+VR)-UR,C*(G*UC+VC)-UC) ENDDO ELSE T = CS/(ONE-SN) DO I=1,N UR=DBLE(X((I-1)*INCX+1)) UC=DIMAG(X((I-1)*INCX+1)) WR=DBLE(Y((I-1)*INCY+1)) WC=DIMAG(Y((I-1)*INCY+1)) VR = CS*(-WC*T+WR)+WC VC = CS*(WR*T+WC)-WR X((I-1)*INCX+1)=DCMPLX(C*(-G*VR+UR)+VR,C*(-G*VC+UC)+VC) Y((I-1)*INCY+1)=DCMPLX(C*(G*UR+VR)-UR,C*(G*UC+VC)-UC) ENDDO ENDIF ENDIF RETURN END SUBROUTINE ZCROTVY SUBROUTINE ZCROTWX(N, X, INCX, Y, INCY, C, G, CS, SN) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX*16 X(*),Y(*) DOUBLE PRECISION G,C,T,WR,WC,UR,UC,VR,VC DOUBLE PRECISION CS, SN DOUBLE PRECISION ZERO, ONE PARAMETER (ONE = 1.0D0, ZERO = 0.0D0) 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=DBLE(X(I)) WC=DIMAG(X(I)) UR = -SN*(WR*T+WC)+WR UC = SN*(-WC*T+WR)+WC VR=DBLE(Y(I)) VC=DIMAG(Y(I)) X(I)=DCMPLX(C*(G*VR+UR)-VR,C*(G*VC+UC)-VC) Y(I)=DCMPLX(C*(-G*UR+VR)+UR,C*(-G*UC+VC)+UC) ENDDO ELSE T = SN/(ONE+CS) DO I=1,N WR=DBLE(X((I-1)*INCX+1)) WC=DIMAG(X((I-1)*INCX+1)) UR = -SN*(WR*T+WC)+WR UC = SN*(-WC*T+WR)+WC VR=DBLE(Y((I-1)*INCY+1)) VC=DIMAG(Y((I-1)*INCY+1)) X((I-1)*INCX+1)=DCMPLX(C*(G*VR+UR)-VR,C*(G*VC+UC)-VC) Y((I-1)*INCY+1)=DCMPLX(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=DBLE(X(I)) WC=DIMAG(X(I)) UR = CS*(WC*T+WR)-WC UC = CS*(-WR*T+WC)+WR VR=DBLE(Y(I)) VC=DIMAG(Y(I)) X(I)=DCMPLX(C*(G*VR+UR)-VR,C*(G*VC+UC)-VC) Y(I)=DCMPLX(C*(-G*UR+VR)+UR,C*(-G*UC+VC)+UC) ENDDO ELSE T = CS/(ONE+SN) DO I=1,N WR=DBLE(X((I-1)*INCX+1)) WC=DIMAG(X((I-1)*INCX+1)) UR = CS*(WC*T+WR)-WC UC = CS*(-WR*T+WC)+WR VR=DBLE(Y((I-1)*INCY+1)) VC=DIMAG(Y((I-1)*INCY+1)) X((I-1)*INCX+1)=DCMPLX(C*(G*VR+UR)-VR,C*(G*VC+UC)-VC) Y((I-1)*INCY+1)=DCMPLX(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=DBLE(X(I)) WC=DIMAG(X(I)) UR = SN*(WR*T-WC)-WR UC = SN*(WC*T+WR)-WC VR=DBLE(Y(I)) VC=DIMAG(Y(I)) X(I)=DCMPLX(C*(G*VR+UR)-VR,C*(G*VC+UC)-VC) Y(I)=DCMPLX(C*(-G*UR+VR)+UR,C*(-G*UC+VC)+UC) ENDDO ELSE T = SN/(ONE-CS) DO I=1,N WR=DBLE(X((I-1)*INCX+1)) WC=DIMAG(X((I-1)*INCX+1)) UR = SN*(WR*T-WC)-WR UC = SN*(WC*T+WR)-WC VR=DBLE(Y((I-1)*INCY+1)) VC=DIMAG(Y((I-1)*INCY+1)) X((I-1)*INCX+1)=DCMPLX(C*(G*VR+UR)-VR,C*(G*VC+UC)-VC) Y((I-1)*INCY+1)=DCMPLX(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=DBLE(X(I)) WC=DIMAG(X(I)) UR = CS*(-WC*T+WR)+WC UC = CS*(WR*T+WC)-WR VR=DBLE(Y(I)) VC=DIMAG(Y(I)) X(I)=DCMPLX(C*(G*VR+UR)-VR,C*(G*VC+UC)-VC) Y(I)=DCMPLX(C*(-G*UR+VR)+UR,C*(-G*UC+VC)+UC) ENDDO ELSE T = CS/(ONE-SN) DO I=1,N WR=DBLE(X((I-1)*INCX+1)) WC=DIMAG(X((I-1)*INCX+1)) UR = CS*(-WC*T+WR)+WC UC = CS*(WR*T+WC)-WR VR=DBLE(Y((I-1)*INCY+1)) VC=DIMAG(Y((I-1)*INCY+1)) X((I-1)*INCX+1)=DCMPLX(C*(G*VR+UR)-VR,C*(G*VC+UC)-VC) Y((I-1)*INCY+1)=DCMPLX(C*(-G*UR+VR)+UR,C*(-G*UC+VC)+UC) ENDDO ENDIF ENDIF RETURN END SUBROUTINE ZCROTWX SUBROUTINE ZCROTWY(N, X, INCX, Y, INCY, C, G, CS, SN) IMPLICIT NONE INTEGER I,N,INCX,INCY COMPLEX*16 X(*),Y(*) DOUBLE PRECISION G,C,T,WR,WC,UR,UC,VR,VC DOUBLE PRECISION CS, SN DOUBLE PRECISION ZERO, ONE PARAMETER (ONE = 1.0D0, ZERO = 0.0D0) 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=DBLE(X(I)) UC=DIMAG(X(I)) WR=DBLE(Y(I)) WC=DIMAG(Y(I)) VR = -SN*(WR*T+WC)+WR VC = SN*(-WC*T+WR)+WC X(I)=DCMPLX(C*(G*VR+UR)-VR,C*(G*VC+UC)-VC) Y(I)=DCMPLX(C*(-G*UR+VR)+UR,C*(-G*UC+VC)+UC) ENDDO ELSE T = SN/(ONE+CS) DO I=1,N UR=DBLE(X((I-1)*INCX+1)) UC=DIMAG(X((I-1)*INCX+1)) WR=DBLE(Y((I-1)*INCY+1)) WC=DIMAG(Y((I-1)*INCY+1)) VR = -SN*(WR*T+WC)+WR VC = SN*(-WC*T+WR)+WC X((I-1)*INCX+1)=DCMPLX(C*(G*VR+UR)-VR,C*(G*VC+UC)-VC) Y((I-1)*INCY+1)=DCMPLX(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=DBLE(X(I)) UC=DIMAG(X(I)) WR=DBLE(Y(I)) WC=DIMAG(Y(I)) VR = CS*(WC*T+WR)-WC VC = CS*(-WR*T+WC)+WR X(I)=DCMPLX(C*(G*VR+UR)-VR,C*(G*VC+UC)-VC) Y(I)=DCMPLX(C*(-G*UR+VR)+UR,C*(-G*UC+VC)+UC) ENDDO ELSE T = CS/(ONE+SN) DO I=1,N UR=DBLE(X((I-1)*INCX+1)) UC=DIMAG(X((I-1)*INCX+1)) WR=DBLE(Y((I-1)*INCY+1)) WC=DIMAG(Y((I-1)*INCY+1)) VR = CS*(WC*T+WR)-WC VC = CS*(-WR*T+WC)+WR X((I-1)*INCX+1)=DCMPLX(C*(G*VR+UR)-VR,C*(G*VC+UC)-VC) Y((I-1)*INCY+1)=DCMPLX(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=DBLE(X(I)) UC=DIMAG(X(I)) WR=DBLE(Y(I)) WC=DIMAG(Y(I)) VR = SN*(WR*T-WC)-WR VC = SN*(WC*T+WR)-WC X(I)=DCMPLX(C*(G*VR+UR)-VR,C*(G*VC+UC)-VC) Y(I)=DCMPLX(C*(-G*UR+VR)+UR,C*(-G*UC+VC)+UC) ENDDO ELSE T = SN/(ONE-CS) DO I=1,N UR=DBLE(X((I-1)*INCX+1)) UC=DIMAG(X((I-1)*INCX+1)) WR=DBLE(Y((I-1)*INCY+1)) WC=DIMAG(Y((I-1)*INCY+1)) VR = SN*(WR*T-WC)-WR VC = SN*(WC*T+WR)-WC X((I-1)*INCX+1)=DCMPLX(C*(G*VR+UR)-VR,C*(G*VC+UC)-VC) Y((I-1)*INCY+1)=DCMPLX(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=DBLE(X(I)) UC=DIMAG(X(I)) WR=DBLE(Y(I)) WC=DIMAG(Y(I)) VR = CS*(-WC*T+WR)+WC VC = CS*(WR*T+WC)-WR X(I)=DCMPLX(C*(G*VR+UR)-VR,C*(G*VC+UC)-VC) Y(I)=DCMPLX(C*(-G*UR+VR)+UR,C*(-G*UC+VC)+UC) ENDDO ELSE T = CS/(ONE-SN) DO I=1,N UR=DBLE(X((I-1)*INCX+1)) UC=DIMAG(X((I-1)*INCX+1)) WR=DBLE(Y((I-1)*INCY+1)) WC=DIMAG(Y((I-1)*INCY+1)) VR = CS*(-WC*T+WR)+WC VC = CS*(WR*T+WC)-WR X((I-1)*INCX+1)=DCMPLX(C*(G*VR+UR)-VR,C*(G*VC+UC)-VC) Y((I-1)*INCY+1)=DCMPLX(C*(-G*UR+VR)+UR,C*(-G*UC+VC)+UC) ENDDO ENDIF ENDIF RETURN END SUBROUTINE ZCROTWY SUBROUTINE ZSCALG(N, X, INCX, CS, SN) IMPLICIT NONE INTEGER N, I, INCX COMPLEX*16 X(*) DOUBLE PRECISION CS, SN, T, XR, XC DOUBLE PRECISION ZERO, ONE PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) IF(CS .GE. ZERO .AND. CS .GE. ABS(SN)) THEN IF (INCX .EQ. 1) THEN T = SN/(ONE+CS) DO I=1,N XR=DBLE(X(I)) XC=DIMAG(X(I)) X(I) = DCMPLX(-SN*(XR*T+XC)+XR, SN*(-XC*T+XR)+XC) ENDDO ELSE T = SN/(ONE+CS) DO I=1,N XR=DBLE(X((I-1)*INCX+1)) XC=DIMAG(X((I-1)*INCX+1)) X((I-1)*INCX+1) = DCMPLX(-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=DBLE(X(I)) XC=DIMAG(X(I)) X(I) = DCMPLX(CS*(XC*T+XR)-XC, CS*(-XR*T+XC)+XR) ENDDO ELSE T = CS/(ONE+SN) DO I=1,N XR=DBLE(X((I-1)*INCX+1)) XC=DIMAG(X((I-1)*INCX+1)) X((I-1)*INCX+1) = DCMPLX(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=DBLE(X(I)) XC=DIMAG(X(I)) X(I) = DCMPLX(SN*(XR*T-XC)-XR, SN*(XC*T+XR)-XC) ENDDO ELSE T = SN/(ONE-CS) DO I=1,N XR=DBLE(X((I-1)*INCX+1)) XC=DIMAG(X((I-1)*INCX+1)) X((I-1)*INCX+1) = DCMPLX(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=DBLE(X(I)) XC=DIMAG(X(I)) X(I) = DCMPLX(CS*(-XC*T+XR)+XC, CS*(XR*T+XC)-XR) ENDDO ELSE T = CS/(ONE-SN) DO I=1,N XR=DBLE(X((I-1)*INCX+1)) XC=DIMAG(X((I-1)*INCX+1)) X((I-1)*INCX+1) = DCMPLX(CS*(-XC*T+XR)+XC, CS*(XR*T+XC)-XR) ENDDO ENDIF ENDIF RETURN END SUBROUTINE ZSCALG SUBROUTINE DLARTGU (F1, G1, C, S, R1) IMPLICIT NONE DOUBLE PRECISION ZERO, ONE, HALF PARAMETER (ONE = 1.0D0, ZERO = 0.0D0, HALF = 0.5E0) DOUBLE PRECISION F1, G1, C, S, R1, U, V EXTERNAL DFMA0 DOUBLE PRECISION DFMA0 IF (F1 .GE. G1) THEN U = G1/F1 R1 = DSQRT(DFMA0(U,U,ONE)) C = ONE/R1 S = U/R1 R1 = R1*F1 ELSE V = F1/G1 R1 = DSQRT(DFMA0(V,V,ONE)) S = ONE/R1 C = V/R1 R1 = R1*G1 ENDIF RETURN END SUBROUTINE DLARTGU SUBROUTINE ZJASVD(N,A,SIGMA,U,V,INFO) use omp_lib IMPLICIT NONE INTEGER N, INFO, FLAG, LL, MM, FLAG2, FLAG3, FLAG4 DOUBLE PRECISION SIGMA(N),B(N),Z(N),TMP,TMP3,TMP4,TMP5 COMPLEX*16 A(N,N),U(N,N),V(N,N) COMPLEX*16 TMP2,TMP6 COMPLEX*16 ZZERO, ZONE PARAMETER (ZONE = (1.0D0,0.0D0), ZZERO = (0.0D0,0.0D0)) DOUBLE PRECISION DZERO, DONE, DHALF PARAMETER (DONE = 1.0D0, DZERO = 0.0D0, DHALF = 0.5E0) INTEGER I, J, K, ITER, MAXITER PARAMETER (MAXITER = 100) DOUBLE PRECISION T, X, ALPHA, BETA DOUBLE PRECISION THETA1, THETA2, CS, SN, CSL, SNL, CSR, SNR DOUBLE PRECISION EPS, R1, R2 DOUBLE PRECISION T1, G1, F1, T2, G2, F2, C1, S1, C2, S2 EXTERNAL DFMA0, DLAMCH DOUBLE PRECISION DFMA0, DLAMCH EPS = DLAMCH('P') DO J=1,N TMP4 = DBLE(A(J,J)) TMP5 = DIMAG(A(J,J)) CALL DLARTG(ABS(TMP4),ABS(TMP5),CSL,SNL,SIGMA(J)) IF (SIGMA(J) .GT. DZERO) THEN CSL = DSIGN(CSL,TMP4) SNL = DSIGN(SNL,TMP5) TMP2 = DCMPLX(CSL,SNL) TMP2 = DCONJG(TMP2) CALL ZSCAL(J-1,TMP2,A(1,J),1) CALL ZSCAL(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*DSQRT(ABS(SIGMA(I)))*DSQRT(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*DSQRT(ABS(SIGMA(I)))*DSQRT(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) = DZERO ENDDO DO J = LL, MM-1 DO K = J+1, MM IF (ABS(A(J,K)) .LE. EPS*DSQRT(ABS(SIGMA(J)))*DSQRT(ABS(SIGMA(K)))) THEN A(J,K) = ZZERO ELSE TMP4 = DBLE(A(J,K)) TMP5 = DIMAG(A(J,K)) CALL DLARTG(ABS(TMP4),ABS(TMP5),CSL,SNL,TMP) CS = DSIGN(CSL,TMP4) SN = DSIGN(SNL,TMP5) F1 = SIGMA(J) - SIGMA(K) CALL DLARTGU (ABS(F1), ABS(TMP), C1, S1, R1) F1 = F1 + DSIGN(R1,F1) IF (F1 .GE. DZERO) THEN G1 = TMP ELSE G1 = -TMP F1 = -F1 ENDIF F2 = SIGMA(J) + SIGMA(K) CALL DLARTGU (ABS(F2), ABS(TMP), C2, S2, R2) F2 = F2 + DSIGN(R2,F2) IF (F2 .GE. DZERO) THEN G2 = -TMP ELSE G2 = TMP 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 = DSIGN(CSL, C1) SNL = DSIGN(SNL, S1) CALL DLARTGU (ABS(C2), ABS(S2), CSR, SNR, R2) CSR = DSIGN(CSR, C2) SNR = DSIGN(SNR, S2) T = CSL + CSR R1 = SNR / T R2 = -SNL / T Z(J) = DFMA0( TMP, R1, Z(J)) Z(K) = DFMA0( 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. DZERO) THEN FLAG2 = 0 CSL = THETA2 SNL = -THETA1 ELSE FLAG2 = 1 CSL = -THETA2 SNL = THETA1 ENDIF THETA1 = CSR THETA2 = SNR IF (THETA2 .GT. DZERO) 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/(DONE+CSL) CSL = DFMA0(-SNL,X,DONE) FLAG = 1 IF (CSL .NE. DONE .OR. SNL .NE. DZERO) THEN CALL ZCROTUY (J-LL+1, A(LL, J), 1, A(LL, K), 1, X, SNL, CS, SN) CALL ZCROTUY (MM-K+1, A(J, K), N, A(K, K), N, X, SNL, CS, SN) CALL ZCROTUY (N, U(1, J), 1, U(1, K), 1, X, SNL, CS, -SN) ELSE CALL ZSCALG (J-LL+1, A(LL, K), 1, CS, SN) CALL ZSCALG (MM-K+1, A(K, K), N, CS, SN) CALL ZSCALG (N, U(1, K), 1, CS, -SN) ENDIF ELSE IF(SNL .GE. DZERO) THEN X = CSL/(DONE+SNL) SNL = DFMA0(-CSL,X,DONE) FLAG = 1 CALL ZCROTVY (J-LL+1, A(LL, J), 1, A(LL, K), 1, CSL, X, CS, SN) CALL ZCROTVY (MM-K+1, A(J, K), N, A(K, K), N, CSL, X, CS, SN) CALL ZCROTVY (N, U(1, J), 1, U(1, K), 1, CSL, X, CS, -SN) ELSE IF(SNL .LE. DZERO) THEN X = CSL/(DONE-SNL) SNL = DFMA0(CSL,X,-DONE) FLAG = 1 CALL ZCROTWY (J-LL+1, A(LL, J), 1, A(LL, K), 1, CSL, X, CS, SN) CALL ZCROTWY (MM-K+1, A(J, K), N, A(K, K), N, CSL, X, CS, SN) CALL ZCROTWY (N, U(1, J), 1, U(1, K), 1, CSL, X, CS, -SN) ENDIF IF (CSR .GE. ABS(SNR)) THEN X = SNR/(DONE+CSR) CSR = DFMA0(-SNR,X,DONE) FLAG = 1 IF (CSR .NE. DONE .OR. SNR .NE. DZERO) THEN CALL ZCROTUY (K-J+1, A(J,J), N, A(J,K), 1, X, SNR, CS, -SN) CALL ZCROTUY (N, V(1, J), 1, V(1, K), 1, X, SNR, CS, -SN) ELSE CALL ZSCALG (K-J+1, A(J, K), 1, CS, -SN) CALL ZSCALG (N, V(1, K), 1, CS, -SN) ENDIF ELSE IF(SNR .GE. DZERO) THEN X = CSR/(DONE+SNR) SNR = DFMA0(-CSR,X,DONE) FLAG = 1 CALL ZCROTVY (K-J+1, A(J,J), N, A(J,K), 1, CSR, X, CS, -SN) CALL ZCROTVY (N, V(1, J), 1, V(1, K), 1, CSR, X, CS, -SN) ELSE IF(SNR .LE. DZERO) THEN X = CSR/(DONE-SNR) SNR = DFMA0(CSR,X,-DONE) FLAG = 1 CALL ZCROTWY (K-J+1, A(J,J), N, A(J,K), 1, CSR, X, CS, -SN) CALL ZCROTWY (N, V(1, J), 1, V(1, K), 1, CSR, X, CS, -SN) ENDIF A(J, K) = ZZERO 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*DSQRT(ABS(SIGMA(I)))*DSQRT(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*DSQRT(ABS(SIGMA(I)))*DSQRT(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) = DZERO ENDDO DO J = LL, MM-1 DO K = J+1, MM IF (ABS(A(J,K)) .LE. EPS*DSQRT(ABS(SIGMA(J)))*DSQRT(ABS(SIGMA(K)))) THEN A(J,K) = ZZERO ELSE TMP4 = DBLE(A(J,K)) TMP5 = DIMAG(A(J,K)) CALL DLARTG(ABS(TMP4),ABS(TMP5),CSL,SNL,TMP) CS = DSIGN(CSL,TMP4) SN = DSIGN(SNL,TMP5) F1 = SIGMA(J) - SIGMA(K) CALL DLARTGU (ABS(F1), ABS(TMP), C1, S1, R1) F1 = F1 + DSIGN(R1,F1) IF (F1 .GE. DZERO) THEN G1 = TMP ELSE G1 = -TMP F1 = -F1 ENDIF F2 = SIGMA(J) + SIGMA(K) CALL DLARTGU (ABS(F2), ABS(TMP), C2, S2, R2) F2 = F2 + DSIGN(R2,F2) IF (F2 .GE. DZERO) THEN G2 = TMP ELSE G2 = -TMP 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 = DSIGN(CSL, C1) SNL = DSIGN(SNL, S1) CALL DLARTGU (ABS(C2), ABS(S2), CSR, SNR, R2) CSR = DSIGN(CSR, C2) SNR = DSIGN(SNR, S2) T = CSL + CSR R1 = SNL / T R2 = -SNR / T Z(J) = DFMA0( TMP, R1, Z(J)) Z(K) = DFMA0( 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. DZERO) THEN FLAG2 = 0 CSL = THETA2 SNL = -THETA1 ELSE FLAG2 = 1 CSL = -THETA2 SNL = THETA1 ENDIF THETA1 = CSR THETA2 = SNR IF (THETA2 .GT. DZERO) 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/(DONE+CSR) CSR = DFMA0(-SNR,X,DONE) FLAG = 1 IF (CSR .NE. DONE .OR. SNR .NE. DZERO) THEN CALL ZCROTUY (J-LL+1, A(LL,J), 1, A(LL,K), 1, X, SNR, CS, SN) CALL ZCROTUY (MM-k+1, A(J,K), N, A(K,K), N, X, SNR, CS, SN) CALL ZCROTUY (N, V(1, J), 1, V(1, K), 1, X, SNR, CS, SN) ELSE CALL ZSCALG (J-LL+1, A(LL, K), 1, CS, SN) CALL ZSCALG (MM-K+1, A(K, K), N, CS, SN) CALL ZSCALG (N, V(1, K), 1, CS, SN) ENDIF ELSE IF(SNR .GE. DZERO) THEN X = CSR/(DONE+SNR) SNR = DFMA0(-CSR,X,DONE) FLAG = 1 CALL ZCROTVY (J-LL+1, A(LL,J), 1, A(LL,K), 1, CSR, X, CS, SN) CALL ZCROTVY (MM-k+1, A(J,K), N, A(K,K), N, CSR, X, CS, SN) CALL ZCROTVY (N, V(1, J), 1, V(1, K), 1, CSR, X, CS, SN) ELSE IF(SNR .LE. DZERO) THEN X = CSR/(DONE-SNR) SNR = DFMA0(CSR,X,-DONE) FLAG = 1 CALL ZCROTWY (J-LL+1, A(LL,J), 1, A(LL,K), 1, CSR, X, CS, SN) CALL ZCROTWY (MM-K+1, A(J,K), N, A(K,K), N, CSR, X, CS, SN) CALL ZCROTWY (N, V(1, J), 1, V(1, K), 1, CSR, X, CS, SN) ENDIF IF (CSL .GE. ABS(SNL)) THEN X = SNL/(DONE+CSL) CSL = DFMA0(-SNL,X,DONE) FLAG = 1 IF (CSL .NE. DONE .OR. SNL .NE. DZERO) THEN CALL ZCROTUY (K-J+1, A(J, J), N, A(J, K), 1, X, SNL, CS, -SN) CALL ZCROTUY (N, U(1, J), 1, U(1, K), 1, X, SNL, CS, SN) ELSE CALL ZSCALG (K-J+1, A(J, K), 1, CS, -SN) CALL ZSCALG (N, U(1, K), 1, CS, SN) ENDIF ELSE IF(SNL .GE. DZERO) THEN X = CSL/(DONE+SNL) SNL = DFMA0(-CSL,X,DONE) FLAG = 1 CALL ZCROTVY (K-J+1, A(j, J), N, A(J, K), 1, CSL, X, CS, -SN) CALL ZCROTVY (N, U(1, J), 1, U(1, K), 1, CSL, X, CS, SN) ELSE IF(SNL .LE. DZERO) THEN X = CSL/(DONE-SNL) SNL = DFMA0(CSL,X,-DONE) FLAG = 1 CALL ZCROTWY (K-J+1, A(J, J), N, A(J, K), 1, CSL, X, CS, -SN) CALL ZCROTWY (N, U(1, J), 1, U(1, K), 1, CSL, X, CS, SN) ENDIF A(J, K) = ZZERO 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*DSQRT(ABS(SIGMA(I)))*DSQRT(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*DSQRT(ABS(SIGMA(I)))*DSQRT(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) = DZERO ENDDO DO K = MM, LL+1, -1 DO J = K-1, LL, -1 IF (ABS(A(J,K)) .LE. EPS*DSQRT(ABS(SIGMA(J)))*DSQRT(ABS(SIGMA(K)))) THEN A(J,K) = ZZERO ELSE TMP4 = DBLE(A(J,K)) TMP5 = DIMAG(A(J,K)) CALL DLARTG(ABS(TMP4),ABS(TMP5),CSL,SNL,TMP) CS = DSIGN(CSL,TMP4) SN = DSIGN(SNL,TMP5) F1 = SIGMA(J) - SIGMA(K) CALL DLARTGU (ABS(F1), ABS(TMP), C1, S1, R1) F1 = F1 + DSIGN(R1,F1) IF (F1 .GE. DZERO) THEN G1 = TMP ELSE G1 = -TMP F1 = -F1 ENDIF F2 = SIGMA(J) + SIGMA(K) CALL DLARTGU (ABS(F2), ABS(TMP), C2, S2, R2) F2 = F2 + DSIGN(R2,F2) IF (F2 .GE. DZERO) THEN G2 = -TMP ELSE G2 = TMP 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 = DSIGN(CSL, C1) SNL = DSIGN(SNL, S1) CALL DLARTGU (ABS(C2), ABS(S2), CSR, SNR, R2) CSR = DSIGN(CSR, C2) SNR = DSIGN(SNR, S2) T = CSL + CSR R1 = SNR / T R2 = -SNL / T Z(J) = DFMA0( TMP, R1, Z(J)) Z(K) = DFMA0( 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. DZERO) THEN FLAG2 = 0 CSL = THETA2 SNL = -THETA1 ELSE FLAG2 = 1 CSL = -THETA2 SNL = THETA1 ENDIF THETA1 = CSR THETA2 = SNR IF (THETA2 .GT. DZERO) 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/(DONE+CSL) CSL = DFMA0(-SNL,X,DONE) FLAG = 1 IF (CSL .NE. DONE .OR. SNL .NE. DZERO) THEN CALL ZCROTUX (K-J+1, A(J, J), N, A(J, K), 1, X, SNL, CS, -SN) CALL ZCROTUX (N, U(1, J), 1, U(1, K), 1, X, SNL, CS, SN) ELSE CALL ZSCALG (K-J+1, A(J, J), N, CS, -SN) CALL ZSCALG (N, U(1, J), 1, CS, SN) ENDIF ELSE IF(SNL .GE. DZERO) THEN X = CSL/(DONE+SNL) SNL = DFMA0(-CSL,X,DONE) FLAG = 1 CALL ZCROTVX (K-J+1, A(J, J), N, A(J, K), 1, CSL, X, CS, -SN) CALL ZCROTVX (N, U(1, J), 1, U(1, K), 1, CSL, X, CS, SN) ELSE IF(SNL .LE. DZERO) THEN X = CSL/(DONE-SNL) SNL = DFMA0(CSL,X,-DONE) FLAG = 1 CALL ZCROTWX (K-J+1, A(J, J), N, A(J, K), 1, CSL, X, CS, -SN) CALL ZCROTWX (N, U(1, J), 1, U(1, K), 1, CSL, X, CS, SN) ENDIF IF (CSR .GE. ABS(SNR)) THEN X = SNR/(DONE+CSR) CSR = DFMA0(-SNR,X,DONE) FLAG = 1 IF (CSR .NE. DONE .OR. SNR .NE. DZERO) THEN CALL ZCROTUX (J-LL+1, A(LL,J), 1, A(LL,K), 1, X, SNR, CS, SN) CALL ZCROTUX (MM-K+1, A(J,K), N, A(K,K), N, X, SNR, CS, SN) CALL ZCROTUX (N, V(1, J), 1, V(1, K), 1, X, SNR, CS, SN) ELSE CALL ZSCALG (J-LL+1, A(LL, J), 1, CS, SN) CALL ZSCALG (MM-K+1, A(J, K), N, CS, SN) CALL ZSCALG (N, V(1, J), 1, CS, SN) ENDIF ELSE IF(SNR .GE. DZERO) THEN X = CSR/(DONE+SNR) SNR = DFMA0(-CSR,X,DONE) FLAG = 1 CALL ZCROTVX (J-LL+1, A(LL,J), 1, A(LL,K), 1, CSR, X, CS, SN) CALL ZCROTVX (MM-K+1, A(J,K), N, A(K,K), N, CSR, X, CS, SN) CALL ZCROTVX (N, V(1, J), 1, V(1, K), 1, CSR, X, CS, SN) ELSE IF(SNR .LE. DZERO) THEN X = CSR/(DONE-SNR) SNR = DFMA0(CSR,X,-DONE) FLAG = 1 CALL ZCROTWX (J-LL+1, A(LL,J), 1, A(LL,K), 1, CSR, X, CS, SN) CALL ZCROTWX (MM-K+1, A(J,K), N, A(K,K), N, CSR, X, CS, SN) CALL ZCROTWX (N, V(1, J), 1, V(1, K), 1, CSR, X, CS, SN) ENDIF A(J, K) = ZZERO 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*DSQRT(ABS(SIGMA(I)))*DSQRT(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*DSQRT(ABS(SIGMA(I)))*DSQRT(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) = DZERO ENDDO DO K = MM, LL+1, -1 DO J = K-1, LL, -1 IF (ABS(A(J,K)) .LE. EPS*DSQRT(ABS(SIGMA(J)))*DSQRT(ABS(SIGMA(K)))) THEN A(J,K) = ZZERO ELSE TMP4 = DBLE(A(J,K)) TMP5 = DIMAG(A(J,K)) CALL DLARTG(ABS(TMP4),ABS(TMP5),CSL,SNL,TMP) CS = DSIGN(CSL,TMP4) SN = DSIGN(SNL,TMP5) F1 = SIGMA(J) - SIGMA(K) CALL DLARTGU (ABS(F1), ABS(TMP), C1, S1, R1) F1 = F1 + DSIGN(R1,F1) IF (F1 .GE. DZERO) THEN G1 = TMP ELSE G1 = -TMP F1 = -F1 ENDIF F2 = SIGMA(J) + SIGMA(K) CALL DLARTGU (ABS(F2), ABS(TMP), C2, S2, R2) F2 = F2 + DSIGN(R2,F2) IF (F2 .GE. DZERO) THEN G2 = TMP ELSE G2 = -TMP 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 = DSIGN(CSL, C1) SNL = DSIGN(SNL, S1) CALL DLARTGU (ABS(C2), ABS(S2), CSR, SNR, R2) CSR = DSIGN(CSR, C2) SNR = DSIGN(SNR, S2) T = CSL + CSR R1 = SNL / T R2 = -SNR / T Z(J) = DFMA0( TMP, R1, Z(J)) Z(K) = DFMA0( 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. DZERO) THEN FLAG2 = 0 CSL = THETA2 SNL = -THETA1 ELSE FLAG2 = 1 CSL = -THETA2 SNL = THETA1 ENDIF THETA1 = CSR THETA2 = SNR IF (THETA2 .GT. DZERO) 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/(DONE+CSR) CSR = DFMA0(-SNR,X,DONE) FLAG = 1 IF (CSR .NE. DONE .OR. SNR .NE. DZERO) THEN CALL ZCROTUX (K-J+1, A(J,J), N, A(J,K), 1, X, SNR, CS, -SN) CALL ZCROTUX (N, V(1, J), 1, V(1, K), 1, X, SNR, CS, -SN) ELSE CALL ZSCALG (K-J+1, A(J, J), N, CS, -SN) CALL ZSCALG (N, V(1, J), 1, CS, -SN) ENDIF ELSE IF(SNR .GE. DZERO) THEN X = CSR/(DONE+SNR) SNR = DFMA0(-CSR,X,DONE) FLAG = 1 CALL ZCROTVX (K-J+1, A(J,J), N, A(J,K), 1, CSR, X, CS, -SN) CALL ZCROTVX (N, V(1, J), 1, V(1, K), 1, CSR, X, CS, -SN) ELSE IF(SNR .LE. DZERO) THEN X = CSR/(DONE-SNR) SNR = DFMA0(CSR,X,-DONE) FLAG = 1 CALL ZCROTWX (K-J+1, A(J,J), N, A(J,K), 1, CSR, X, CS, -SN) CALL ZCROTWX (N, V(1, J), 1, V(1, K), 1, CSR, X, CS, -SN) ENDIF IF (CSL .GE. ABS(SNL)) THEN X = SNL/(DONE+CSL) CSL = DFMA0(-SNL,X,DONE) FLAG = 1 IF (CSL .NE. DONE .OR. SNL .NE. DZERO) THEN CALL ZCROTUX (J-LL+1, A(LL, J), 1, A(LL, K), 1, X, SNL, CS, SN) CALL ZCROTUX (MM-K+1, A(J, K), N, A(K, K), N, X, SNL, CS, SN) CALL ZCROTUX (N, U(1, J), 1, U(1, K), 1, X, SNL, CS, -SN) ELSE CALL ZSCALG (J-LL+1, A(LL, J), 1, CS, SN) CALL ZSCALG (MM-K+1, A(J, K), N, CS, SN) CALL ZSCALG (N, U(1, J), 1, CS, -SN) ENDIF ELSE IF(SNL .GE. DZERO) THEN X = CSL/(DONE+SNL) SNL = DFMA0(-CSL,X,DONE) FLAG = 1 CALL ZCROTVX (J-LL+1, A(LL, J), 1, A(LL, K), 1, CSL, X, CS, SN) CALL ZCROTVX (MM-K+1, A(J, K), N, A(K, K), N, CSL, X, CS, SN) CALL ZCROTVX (N, U(1, J), 1, U(1, K), 1, CSL, X, CS, -SN) ELSE IF(SNL .LE. DZERO) THEN X = CSL/(DONE-SNL) SNL = DFMA0(CSL,X,-DONE) FLAG = 1 CALL ZCROTWX (J-LL+1, A(LL, J), 1, A(LL, K), 1, CSL, X, CS, SN) CALL ZCROTWX (MM-K+1, A(J, K), N, A(K, K), N, CSL, X, CS, SN) CALL ZCROTWX (N, U(1, J), 1, U(1, K), 1, CSL, X, CS, -SN) ENDIF A(J, K) = ZZERO 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. DZERO) THEN SIGMA(I) = -SIGMA(I) CALL ZSCAL(N, -ZONE, 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 ZSWAP(N, U(1,K), 1, U(1,I), 1) CALL ZSWAP(N, V(1,K), 1, V(1,I), 1) ENDDO INFO = 0 RETURN END SUBROUTINE ZJASVD PROGRAM TESTDJS use omp_lib IMPLICIT NONE INTEGER I, J, K, N, INFO, SEED DOUBLE PRECISION TMP1, TMP2, PI COMPLEX*16 E, S COMPLEX*16 ZZERO, ZONE PARAMETER (ZZERO = (0.0D0,0.0D0), ZONE=(1.0D0,0.0D0)) DOUBLE PRECISION DZERO, DONE, DTWO PARAMETER (DZERO = 0.0D0, DONE=1.0D0, DTWO=2.0D0) DOUBLE PRECISION,ALLOCATABLE :: SIGMA(:) COMPLEX*16,ALLOCATABLE :: A(:,:), B(:,:), U(:,:), V(:,:) EXTERNAL USERRAND DOUBLE PRECISION USERRAND EXTERNAL ZDOTC COMPLEX*16 ZDOTC character*1000 ARG1,ARG2 integer t1, t2, t_rate, t_max, diff PI = DACOS(-DONE) 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) = DCMPLX(TMP1,TMP2) B(J,I) = A(J,I) ENDDO DO J=I+1,N A(J,I) = ZZERO B(J,I) = ZZERO ENDDO ENDDO DO J=1,N DO K=1,N U(K,J) = ZZERO ENDDO U(J,J) = ZONE ENDDO DO J=1,N DO K=1,N V(K,J) = ZZERO ENDDO V(J,J) = ZONE ENDDO call system_clock(t1) !開始時を記録 CALL ZJASVD(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 = ZZERO !$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',DSQRT(DBLE(E)) CALL ZHERK ('U', 'C', N, N, DONE, U, N, DZERO, A, N) DO I=1,N A(I,I) = A(I,I) - ZONE ENDDO E = ZZERO !$OMP PARALLEL DO REDUCTION(+:E) DO I=2,N E = E+ZDOTC(I-1,A(1,I),1,A(1,I),1) ENDDO !$OMP END PARALLEL DO E = DTWO*E DO I=1,N E = E+CONJG(A(I,I))*A(I,I) ENDDO WRITE (*,*) 'ERROR',DSQRT(DBLE(E)) CALL ZHERK ('U', 'C', N, N, DONE, V, N, DZERO, A, N) DO I=1,N A(I,I) = A(I,I) - ZONE ENDDO E = ZZERO !$OMP PARALLEL DO REDUCTION(+:E) DO I=2,N E = E+ZDOTC(I-1,A(1,I),1,A(1,I),1) ENDDO !$OMP END PARALLEL DO E = DTWO*E DO I=1,N E = E+CONJG(A(I,I))*A(I,I) ENDDO WRITE (*,*) 'ERROR',DSQRT(DBLE(E)) END PROGRAM TESTDJS