mirror of https://gitlab.com/QEF/q-e.git
667 lines
19 KiB
Fortran
667 lines
19 KiB
Fortran
!
|
|
! Copyright (C) 2002-2005 FPMD-CPV groups
|
|
! This file is distributed under the terms of the
|
|
! GNU General Public License. See the file `License'
|
|
! in the root directory of the present distribution,
|
|
! or http://www.gnu.org/copyleft/gpl.txt .
|
|
!
|
|
! AB INITIO COSTANT PRESSURE MOLECULAR DYNAMICS
|
|
! ----------------------------------------------
|
|
! Car-Parrinello Parallel Program
|
|
! Carlo Cavazzoni - Gerardo Ballabio
|
|
! SISSA, Trieste, Italy - 1997-99
|
|
! Last modified: Sat Nov 13 11:04:11 MET 1999
|
|
! ----------------------------------------------
|
|
|
|
#include "f_defs.h"
|
|
|
|
! routines in this file:
|
|
! SUBROUTINE geninv(a,ld,n,mrank,cond,u,v,work,toleig,info,iopt)
|
|
! SUBROUTINE zgeninv(a,ld,n,mrank,cond,u,v,work,toleig,info,iopt)
|
|
! ----------------------------------------------
|
|
! ----------------------------------------------
|
|
! BEGIN manual
|
|
|
|
SUBROUTINE geninv( a, ld, n, mrank, cond, u, v, &
|
|
work, toleig, info, iopt)
|
|
|
|
! get a general inverse matrix
|
|
!
|
|
! iopt = 0: using calculated rank for pseudo inverse
|
|
! iopt = 1 : rank is assumed to be n-6 for pseudo inverse
|
|
! iopt > 10: scale matrix before decomposition
|
|
! ----------------------------------------------
|
|
! END manual
|
|
! end of declarations
|
|
! ----------------------------------------------
|
|
|
|
USE kinds
|
|
IMPLICIT NONE
|
|
|
|
INTEGER :: ld, n, mrank, info, iopt
|
|
REAL(DP) :: cond, toleig
|
|
|
|
REAL(DP) :: a(ld,n), u(ld,n), v(ld,n), work(4*n)
|
|
REAL(DP) :: zero=0.0d0
|
|
REAL(DP) :: one=1.0d0
|
|
|
|
INTEGER :: n3, i, j, k, m
|
|
|
|
! ... scale matrix before inversion
|
|
IF ( iopt >= 10 ) THEN
|
|
n3 = 3 * n
|
|
DO i = 1, n
|
|
work( n3 + i ) = one
|
|
IF ( ABS( a(i,i) ) >= 1.d-13 ) &
|
|
work( n3 + i ) = one / SQRT( ABS( a(i,i) ) )
|
|
END DO
|
|
DO i = 1, n
|
|
DO j = 1, n
|
|
a(i,j) = a(i,j) * work( n3 + i ) * work( n3 + j )
|
|
END DO
|
|
END DO
|
|
END IF
|
|
|
|
! ... get singular values
|
|
CALL dsvdc( a, ld, n, n, work, work(n+1), u, ld, v, ld, work(2*n+1), 11, info)
|
|
mrank = 0
|
|
DO i = 1, n
|
|
IF ( ABS( work(i) ) > toleig ) mrank = mrank + 1
|
|
END DO
|
|
|
|
m = mrank
|
|
IF ( iopt == 1 .OR. iopt == 11 ) m = n - 6
|
|
cond = work(1) / work(m)
|
|
DO i = 1, m
|
|
work(i) = one / work(i)
|
|
END DO
|
|
|
|
DO j = 1, n
|
|
DO i = 1, n
|
|
a(i,j) = zero
|
|
END DO
|
|
DO k = 1, m
|
|
DO i = 1, n
|
|
a(i,j) = a(i,j) + v(i,k) * work(k) * u(j,k)
|
|
END DO
|
|
END DO
|
|
END DO
|
|
|
|
! ... rescale matrix after inversion
|
|
IF ( iopt >= 10 ) THEN
|
|
DO i = 1, n
|
|
DO j = 1, n
|
|
a(i,j) = a(i,j) * work( n3 + i ) * work( n3 + j )
|
|
END DO
|
|
END DO
|
|
END IF
|
|
|
|
RETURN
|
|
END SUBROUTINE geninv
|
|
|
|
! ----------------------------------------------
|
|
! ----------------------------------------------
|
|
! BEGIN manual
|
|
SUBROUTINE zgeninv(a,ld,n,mrank,cond,u,v,work,toleig,info,iopt)
|
|
|
|
! get a general inverse matrix
|
|
!
|
|
! iopt = 0: using calculated rank for pseudo inverse
|
|
! iopt = 1 : rank is assumed to be n-6 for pseudo inverse
|
|
! iopt > 10: scale matrix before decomposition
|
|
! ----------------------------------------------
|
|
! END manual
|
|
! end of declarations
|
|
! ----------------------------------------------
|
|
|
|
USE kinds
|
|
IMPLICIT NONE
|
|
|
|
INTEGER :: ld, n, mrank, info, iopt
|
|
REAL(DP) :: cond, toleig
|
|
|
|
COMPLEX(DP) a(ld,n),u(ld,n),v(ld,n),work(4*n)
|
|
REAL(DP) :: zero=0.0d0
|
|
REAL(DP) :: one=1.0d0
|
|
|
|
INTEGER :: n3, i, j, k, m
|
|
|
|
! ... scale matrix before inversion
|
|
IF (iopt.GE.10) THEN
|
|
n3=3*n
|
|
DO i=1,n
|
|
work(n3+i)=one
|
|
IF (abs(a(i,i)).GE.1.d-13) &
|
|
work(n3+i)=one/dsqrt(abs(a(i,i)))
|
|
END DO
|
|
DO i=1,n
|
|
DO j=1,n
|
|
a(i,j)=a(i,j)*work(n3+i)*work(n3+j)
|
|
END DO
|
|
END DO
|
|
END IF
|
|
|
|
! ... get singular values
|
|
CALL dsvdc(a,ld,n,n,work,work(n+1),u,ld,v,ld,work(2*n+1),11,info)
|
|
mrank=0
|
|
DO i=1,n
|
|
IF (abs(work(i)).GT.toleig) mrank=mrank+1
|
|
END DO
|
|
|
|
m=mrank
|
|
IF (iopt.EQ.1.OR.iopt.EQ.11) m=n-6
|
|
cond=work(1)/work(m)
|
|
DO i=1,m
|
|
work(i)=one/work(i)
|
|
END DO
|
|
|
|
DO j=1,n
|
|
DO i=1,n
|
|
a(i,j)=zero
|
|
END DO
|
|
DO k=1,m
|
|
DO i=1,n
|
|
a(i,j)=a(i,j)+v(i,k)*work(k)*u(j,k)
|
|
END DO
|
|
END DO
|
|
END DO
|
|
|
|
! ... rescale matrix after inversion
|
|
IF (iopt.GE.10) THEN
|
|
DO i=1,n
|
|
DO j=1,n
|
|
a(i,j)=a(i,j)*work(n3+i)*work(n3+j)
|
|
END DO
|
|
END DO
|
|
END IF
|
|
|
|
RETURN
|
|
END SUBROUTINE zgeninv
|
|
|
|
! ----------------------------------------------
|
|
! ----------------------------------------------
|
|
SUBROUTINE dsvdc(x,ldx,n,p,s,e,u,ldu,v,ldv,work,job,info)
|
|
|
|
! (describe briefly what this routine does...)
|
|
! ----------------------------------------------
|
|
|
|
! DSVDC IS A SUBROUTINE TO REDUCE A DOUBLE PRECISION NXP MATRIX X
|
|
! BY ORTHOGONAL TRANSFORMATIONS U AND V TO DIAGONAL FORM. THE
|
|
! DIAGONAL ELEMENTS S(I) ARE THE SINGULAR VALUES OF X. THE
|
|
! COLUMNS OF U ARE THE CORRESPONDING LEFT SINGULAR VECTORS,
|
|
! AND THE COLUMNS OF V THE RIGHT SINGULAR VECTORS.
|
|
!
|
|
! ON ENTRY
|
|
!
|
|
! X DOUBLE PRECISION(LDX,P), WHERE LDX.GE.N.
|
|
! X CONTAINS THE MATRIX WHOSE SINGULAR VALUE
|
|
! DECOMPOSITION IS TO BE COMPUTED. X IS
|
|
! DESTROYED BY DSVDC.
|
|
!
|
|
! LDX INTEGER.
|
|
! LDX IS THE LEADING DIMENSION OF THE ARRAY X.
|
|
!
|
|
! N INTEGER.
|
|
! N IS THE NUMBER OF ROWS OF THE MATRIX X.
|
|
!
|
|
! P INTEGER.
|
|
! P IS THE NUMBER OF COLUMNS OF THE MATRIX X.
|
|
!
|
|
! LDU INTEGER.
|
|
! LDU IS THE LEADING DIMENSION OF THE ARRAY U.
|
|
! (SEE BELOW).
|
|
!
|
|
! LDV INTEGER.
|
|
! LDV IS THE LEADING DIMENSION OF THE ARRAY V.
|
|
! (SEE BELOW).
|
|
!
|
|
! WORK DOUBLE PRECISION(N).
|
|
! WORK IS A SCRATCH ARRAY.
|
|
!
|
|
! JOB INTEGER.
|
|
! JOB CONTROLS THE COMPUTATION OF THE SINGULAR
|
|
! VECTORS. IT HAS THE DECIMAL EXPANSION AB
|
|
! WITH THE FOLLOWING MEANING
|
|
!
|
|
! A.EQ.0 DO NOT COMPUTE THE LEFT SINGULAR
|
|
! VECTORS.
|
|
! A.EQ.1 RETURN THE N LEFT SINGULAR VECTORS
|
|
! IN U.
|
|
! A.GE.2 RETURN THE FIRST MIN(N,P) SINGULAR
|
|
! VECTORS IN U.
|
|
! B.EQ.0 DO NOT COMPUTE THE RIGHT SINGULAR
|
|
! VECTORS.
|
|
! B.EQ.1 RETURN THE RIGHT SINGULAR VECTORS
|
|
! IN V.
|
|
!
|
|
! ON RETURN
|
|
!
|
|
! S DOUBLE PRECISION(MM), WHERE MM=MIN(N+1,P).
|
|
! THE FIRST MIN(N,P) ENTRIES OF S CONTAIN THE
|
|
! SINGULAR VALUES OF X ARRANGED IN DESCENDING
|
|
! ORDER OF MAGNITUDE.
|
|
!
|
|
! E DOUBLE PRECISION(P),
|
|
! E ORDINARILY CONTAINS ZEROS. HOWEVER SEE THE
|
|
! DISCUSSION OF INFO FOR EXCEPTIONS.
|
|
!
|
|
! U DOUBLE PRECISION(LDU,K), WHERE LDU.GE.N. IF
|
|
! JOBA.EQ.1 THEN K.EQ.N, IF JOBA.GE.2
|
|
! THEN K.EQ.MIN(N,P).
|
|
! U CONTAINS THE MATRIX OF LEFT SINGULAR VECTORS.
|
|
! U IS NOT REFERENCED IF JOBA.EQ.0. IF N.LE.P
|
|
! OR IF JOBA.EQ.2, THEN U MAY BE IDENTIFIED WITH X
|
|
! IN THE SUBROUTINE CALL.
|
|
!
|
|
! V DOUBLE PRECISION(LDV,P), WHERE LDV.GE.P.
|
|
! V CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS.
|
|
! V IS NOT REFERENCED IF JOB.EQ.0. IF P.LE.N,
|
|
! THEN V MAY BE IDENTIFIED WITH X IN THE
|
|
! SUBROUTINE CALL.
|
|
!
|
|
! INFO INTEGER.
|
|
! THE SINGULAR VALUES (AND THEIR CORRESPONDING
|
|
! SINGULAR VECTORS) S(INFO+1),S(INFO+2),...,S(M)
|
|
! ARE CORRECT (HERE M=MIN(N,P)). THUS IF
|
|
! INFO.EQ.0, ALL THE SINGULAR VALUES AND THEIR
|
|
! VECTORS ARE CORRECT. IN ANY EVENT, THE MATRIX
|
|
! B = TRANS(U)*X*V IS THE BIDIAGONAL MATRIX
|
|
! WITH THE ELEMENTS OF S ON ITS DIAGONAL AND THE
|
|
! ELEMENTS OF E ON ITS SUPER-DIAGONAL (TRANS(U)
|
|
! IS THE TRANSPOSE OF U). THUS THE SINGULAR
|
|
! VALUES OF X AND B ARE THE SAME.
|
|
!
|
|
! LINPACK. THIS VERSION DATED 08/14/78 .
|
|
! CORRECTION MADE TO SHIFT 2/84.
|
|
! G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
|
|
!
|
|
! DSVDC USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS.
|
|
!
|
|
! EXTERNAL DROT
|
|
! BLAS DAXPY,DDOT,DSCAL,DSWAP,DNRM2,DROTG
|
|
! FORTRAN DABS,DMAX1,MAX0,MIN0,MOD,DSQRT
|
|
!
|
|
|
|
USE kinds
|
|
INTEGER LDX,N,P,LDU,LDV,JOB,INFO
|
|
REAL(DP) X(LDX,1),S(1),E(1),U(LDU,1),V(LDV,1),WORK(1)
|
|
|
|
! INTERNAL VARIABLES
|
|
!
|
|
INTEGER I,ITER,J,JOBU,K,KASE,KK,L,LL,LLS,LM1,LP1,LS,LU,M,MAXIT, &
|
|
& MM,MM1,MP1,NCT,NCTP1,NCU,NRT,NRTP1
|
|
REAL(DP) DDOT,T,R
|
|
REAL(DP) B,C,CS,EL,EMM1,F,G,DNRM2,SCALE,SHIFT,SL,SM,SN, &
|
|
& SMM1,T1,TEST,ZTEST
|
|
LOGICAL WANTU,WANTV
|
|
!
|
|
!
|
|
! SET THE MAXIMUM NUMBER OF ITERATIONS.
|
|
!
|
|
MAXIT = 30
|
|
!
|
|
! DETERMINE WHAT IS TO BE COMPUTED.
|
|
!
|
|
WANTU = .FALSE.
|
|
WANTV = .FALSE.
|
|
JOBU = MOD(JOB,100)/10
|
|
NCU = N
|
|
IF (JOBU .GT. 1) NCU = MIN0(N,P)
|
|
IF (JOBU .NE. 0) WANTU = .TRUE.
|
|
IF (MOD(JOB,10) .NE. 0) WANTV = .TRUE.
|
|
!
|
|
! REDUCE X TO BIDIAGONAL FORM, STORING THE DIAGONAL ELEMENTS
|
|
! IN S AND THE SUPER-DIAGONAL ELEMENTS IN E.
|
|
!
|
|
INFO = 0
|
|
NCT = MIN0(N-1,P)
|
|
NRT = MAX0(0,MIN0(P-2,N))
|
|
LU = MAX0(NCT,NRT)
|
|
IF (LU .LT. 1) GO TO 170
|
|
DO 160 L = 1, LU
|
|
LP1 = L + 1
|
|
IF (L .GT. NCT) GO TO 20
|
|
!
|
|
! COMPUTE THE TRANSFORMATION FOR THE L-TH COLUMN AND
|
|
! PLACE THE L-TH DIAGONAL IN S(L).
|
|
!
|
|
s(l) = DNRM2(n-l+1,x(l,l),1)
|
|
IF (S(L) .EQ. 0.0D0) GO TO 10
|
|
IF (X(L,L) .NE. 0.0D0) S(L) = SIGN(S(L),X(L,L))
|
|
call DSCAL(n-l+1,1.0d0/s(l),x(l,l),1)
|
|
X(L,L) = 1.0D0 + X(L,L)
|
|
10 CONTINUE
|
|
S(L) = -S(L)
|
|
20 CONTINUE
|
|
IF (P .LT. LP1) GO TO 50
|
|
DO 40 J = LP1, P
|
|
IF (L .GT. NCT) GO TO 30
|
|
IF (S(L) .EQ. 0.0D0) GO TO 30
|
|
!
|
|
! APPLY THE TRANSFORMATION.
|
|
!
|
|
t = - DDOT(n-l+1,x(l,l),1,x(l,j),1)/x(l,l)
|
|
call DAXPY(n-l+1,t,x(l,l),1,x(l,j),1)
|
|
30 CONTINUE
|
|
!
|
|
! PLACE THE L-TH ROW OF X INTO E FOR THE
|
|
! SUBSEQUENT CALCULATION OF THE ROW TRANSFORMATION.
|
|
!
|
|
E(J) = X(L,J)
|
|
40 CONTINUE
|
|
50 CONTINUE
|
|
IF (.NOT.WANTU .OR. L .GT. NCT) GO TO 70
|
|
!
|
|
! PLACE THE TRANSFORMATION IN U FOR SUBSEQUENT BACK
|
|
! MULTIPLICATION.
|
|
!
|
|
DO 60 I = L, N
|
|
U(I,L) = X(I,L)
|
|
60 CONTINUE
|
|
70 CONTINUE
|
|
IF (L .GT. NRT) GO TO 150
|
|
!
|
|
! COMPUTE THE L-TH ROW TRANSFORMATION AND PLACE THE
|
|
! L-TH SUPER-DIAGONAL IN E(L).
|
|
!
|
|
e(l) = DNRM2(p-l,e(lp1),1)
|
|
IF (E(L) .EQ. 0.0D0) GO TO 80
|
|
IF (E(LP1) .NE. 0.0D0) E(L) = SIGN(E(L),E(LP1))
|
|
call DSCAL(p-l,1.0d0/e(l),e(lp1),1)
|
|
E(LP1) = 1.0D0 + E(LP1)
|
|
80 CONTINUE
|
|
E(L) = -E(L)
|
|
IF (LP1 .GT. N .OR. E(L) .EQ. 0.0D0) GO TO 120
|
|
!
|
|
! APPLY THE TRANSFORMATION.
|
|
!
|
|
DO 90 I = LP1, N
|
|
WORK(I) = 0.0D0
|
|
90 CONTINUE
|
|
DO 100 J = LP1, P
|
|
call DAXPY(n-l,e(j),x(lp1,j),1,work(lp1),1)
|
|
100 CONTINUE
|
|
DO 110 J = LP1, P
|
|
call DAXPY(n-l,-e(j)/e(lp1),work(lp1),1,x(lp1,j),1)
|
|
110 CONTINUE
|
|
120 CONTINUE
|
|
IF (.NOT.WANTV) GO TO 140
|
|
!
|
|
! PLACE THE TRANSFORMATION IN V FOR SUBSEQUENT
|
|
! BACK MULTIPLICATION.
|
|
!
|
|
DO 130 I = LP1, P
|
|
V(I,L) = E(I)
|
|
130 CONTINUE
|
|
140 CONTINUE
|
|
150 CONTINUE
|
|
160 CONTINUE
|
|
170 CONTINUE
|
|
!
|
|
! SET UP THE FINAL BIDIAGONAL MATRIX OR ORDER M.
|
|
!
|
|
M = MIN0(P,N+1)
|
|
NCTP1 = NCT + 1
|
|
NRTP1 = NRT + 1
|
|
IF (NCT .LT. P) S(NCTP1) = X(NCTP1,NCTP1)
|
|
IF (N .LT. M) S(M) = 0.0D0
|
|
IF (NRTP1 .LT. M) E(NRTP1) = X(NRTP1,M)
|
|
E(M) = 0.0D0
|
|
!
|
|
! IF REQUIRED, GENERATE U.
|
|
!
|
|
IF (.NOT.WANTU) GO TO 300
|
|
IF (NCU .LT. NCTP1) GO TO 200
|
|
DO 190 J = NCTP1, NCU
|
|
DO 180 I = 1, N
|
|
U(I,J) = 0.0D0
|
|
180 CONTINUE
|
|
U(J,J) = 1.0D0
|
|
190 CONTINUE
|
|
200 CONTINUE
|
|
IF (NCT .LT. 1) GO TO 290
|
|
DO 280 LL = 1, NCT
|
|
L = NCT - LL + 1
|
|
IF (S(L) .EQ. 0.0D0) GO TO 250
|
|
LP1 = L + 1
|
|
IF (NCU .LT. LP1) GO TO 220
|
|
DO 210 J = LP1, NCU
|
|
t = - DDOT(n-l+1,u(l,l),1,u(l,j),1)/u(l,l)
|
|
call DAXPY(n-l+1,t,u(l,l),1,u(l,j),1)
|
|
210 CONTINUE
|
|
220 CONTINUE
|
|
call DSCAL(n-l+1,-1.0d0,u(l,l),1)
|
|
U(L,L) = 1.0D0 + U(L,L)
|
|
LM1 = L - 1
|
|
IF (LM1 .LT. 1) GO TO 240
|
|
DO 230 I = 1, LM1
|
|
U(I,L) = 0.0D0
|
|
230 CONTINUE
|
|
240 CONTINUE
|
|
GO TO 270
|
|
250 CONTINUE
|
|
DO 260 I = 1, N
|
|
U(I,L) = 0.0D0
|
|
260 CONTINUE
|
|
U(L,L) = 1.0D0
|
|
270 CONTINUE
|
|
280 CONTINUE
|
|
290 CONTINUE
|
|
300 CONTINUE
|
|
!
|
|
! IF IT IS REQUIRED, GENERATE V.
|
|
!
|
|
IF (.NOT.WANTV) GO TO 350
|
|
DO 340 LL = 1, P
|
|
L = P - LL + 1
|
|
LP1 = L + 1
|
|
IF (L .GT. NRT) GO TO 320
|
|
IF (E(L) .EQ. 0.0D0) GO TO 320
|
|
DO 310 J = LP1, P
|
|
t = - DDOT(p-l,v(lp1,l),1,v(lp1,j),1)/v(lp1,l)
|
|
call DAXPY(p-l,t,v(lp1,l),1,v(lp1,j),1)
|
|
310 CONTINUE
|
|
320 CONTINUE
|
|
DO 330 I = 1, P
|
|
V(I,L) = 0.0D0
|
|
330 CONTINUE
|
|
V(L,L) = 1.0D0
|
|
340 CONTINUE
|
|
350 CONTINUE
|
|
!
|
|
! MAIN ITERATION LOOP FOR THE SINGULAR VALUES.
|
|
!
|
|
MM = M
|
|
ITER = 0
|
|
360 CONTINUE
|
|
!
|
|
! QUIT IF ALL THE SINGULAR VALUES HAVE BEEN FOUND.
|
|
!
|
|
! ...EXIT
|
|
IF (M .EQ. 0) GO TO 620
|
|
!
|
|
! IF TOO MANY ITERATIONS HAVE BEEN PERFORMED, SET
|
|
! FLAG AND RETURN.
|
|
!
|
|
IF (ITER .LT. MAXIT) GO TO 370
|
|
INFO = M
|
|
! ......EXIT
|
|
GO TO 620
|
|
370 CONTINUE
|
|
!
|
|
! THIS SECTION OF THE PROGRAM INSPECTS FOR
|
|
! NEGLIGIBLE ELEMENTS IN THE S AND E ARRAYS. ON
|
|
! COMPLETION THE VARIABLES KASE AND L ARE SET AS FOLLOWS.
|
|
!
|
|
! KASE = 1 IF S(M) AND E(L-1) ARE NEGLIGIBLE AND L.LT.M
|
|
! KASE = 2 IF S(L) IS NEGLIGIBLE AND L.LT.M
|
|
! KASE = 3 IF E(L-1) IS NEGLIGIBLE, L.LT.M, AND
|
|
! S(L), ..., S(M) ARE NOT NEGLIGIBLE (QR STEP).
|
|
! KASE = 4 IF E(M-1) IS NEGLIGIBLE (CONVERGENCE).
|
|
!
|
|
DO 390 LL = 1, M
|
|
L = M - LL
|
|
! ...EXIT
|
|
IF (L .EQ. 0) GO TO 400
|
|
TEST = DABS(S(L)) + DABS(S(L+1))
|
|
ZTEST = TEST + DABS(E(L))
|
|
IF (ZTEST .NE. TEST) GO TO 380
|
|
E(L) = 0.0D0
|
|
! ......EXIT
|
|
GO TO 400
|
|
380 CONTINUE
|
|
390 CONTINUE
|
|
400 CONTINUE
|
|
IF (L .NE. M - 1) GO TO 410
|
|
KASE = 4
|
|
GO TO 480
|
|
410 CONTINUE
|
|
LP1 = L + 1
|
|
MP1 = M + 1
|
|
DO 430 LLS = LP1, MP1
|
|
LS = M - LLS + LP1
|
|
! ...EXIT
|
|
IF (LS .EQ. L) GO TO 440
|
|
TEST = 0.0D0
|
|
IF (LS .NE. M) TEST = TEST + DABS(E(LS))
|
|
IF (LS .NE. L + 1) TEST = TEST + DABS(E(LS-1))
|
|
ZTEST = TEST + DABS(S(LS))
|
|
IF (ZTEST .NE. TEST) GO TO 420
|
|
S(LS) = 0.0D0
|
|
! ......EXIT
|
|
GO TO 440
|
|
420 CONTINUE
|
|
430 CONTINUE
|
|
440 CONTINUE
|
|
IF (LS .NE. L) GO TO 450
|
|
KASE = 3
|
|
GO TO 470
|
|
450 CONTINUE
|
|
IF (LS .NE. M) GO TO 460
|
|
KASE = 1
|
|
GO TO 470
|
|
460 CONTINUE
|
|
KASE = 2
|
|
L = LS
|
|
470 CONTINUE
|
|
480 CONTINUE
|
|
L = L + 1
|
|
!
|
|
! PERFORM THE TASK INDICATED BY KASE.
|
|
!
|
|
GO TO (490,520,540,570), KASE
|
|
!
|
|
! DEFLATE NEGLIGIBLE S(M).
|
|
!
|
|
490 CONTINUE
|
|
MM1 = M - 1
|
|
F = E(M-1)
|
|
E(M-1) = 0.0D0
|
|
DO 510 KK = L, MM1
|
|
K = MM1 - KK + L
|
|
T1 = S(K)
|
|
CALL DROTG(T1,F,CS,SN)
|
|
S(K) = T1
|
|
IF (K .EQ. L) GO TO 500
|
|
F = -SN*E(K-1)
|
|
E(K-1) = CS*E(K-1)
|
|
500 CONTINUE
|
|
IF (WANTV) CALL DROT(P,V(1,K),1,V(1,M),1,CS,SN)
|
|
510 CONTINUE
|
|
GO TO 610
|
|
!
|
|
! SPLIT AT NEGLIGIBLE S(L).
|
|
!
|
|
520 CONTINUE
|
|
F = E(L-1)
|
|
E(L-1) = 0.0D0
|
|
DO 530 K = L, M
|
|
T1 = S(K)
|
|
CALL DROTG(T1,F,CS,SN)
|
|
S(K) = T1
|
|
F = -SN*E(K)
|
|
E(K) = CS*E(K)
|
|
IF (WANTU) CALL DROT(N,U(1,K),1,U(1,L-1),1,CS,SN)
|
|
530 CONTINUE
|
|
GO TO 610
|
|
!
|
|
! PERFORM ONE QR STEP.
|
|
!
|
|
540 CONTINUE
|
|
!
|
|
! CALCULATE THE SHIFT.
|
|
!
|
|
SCALE = DMAX1(DABS(S(M)),DABS(S(M-1)),DABS(E(M-1)), &
|
|
& DABS(S(L)),DABS(E(L)))
|
|
SM = S(M)/SCALE
|
|
SMM1 = S(M-1)/SCALE
|
|
EMM1 = E(M-1)/SCALE
|
|
SL = S(L)/SCALE
|
|
EL = E(L)/SCALE
|
|
B = ((SMM1 + SM)*(SMM1 - SM) + EMM1**2)/2.0D0
|
|
C = (SM*EMM1)**2
|
|
SHIFT = 0.0D0
|
|
IF (B .EQ. 0.0D0 .AND. C .EQ. 0.0D0) GO TO 550
|
|
SHIFT = DSQRT(B**2+C)
|
|
IF (B .LT. 0.0D0) SHIFT = -SHIFT
|
|
SHIFT = C/(B + SHIFT)
|
|
550 CONTINUE
|
|
F = (SL + SM)*(SL - SM) + SHIFT
|
|
G = SL*EL
|
|
!
|
|
! CHASE ZEROS.
|
|
!
|
|
MM1 = M - 1
|
|
DO 560 K = L, MM1
|
|
CALL DROTG(F,G,CS,SN)
|
|
IF (K .NE. L) E(K-1) = F
|
|
F = CS*S(K) + SN*E(K)
|
|
E(K) = CS*E(K) - SN*S(K)
|
|
G = SN*S(K+1)
|
|
S(K+1) = CS*S(K+1)
|
|
IF (WANTV) CALL DROT(P,V(1,K),1,V(1,K+1),1,CS,SN)
|
|
CALL DROTG(F,G,CS,SN)
|
|
S(K) = F
|
|
F = CS*E(K) + SN*S(K+1)
|
|
S(K+1) = -SN*E(K) + CS*S(K+1)
|
|
G = SN*E(K+1)
|
|
E(K+1) = CS*E(K+1)
|
|
IF (WANTU .AND. K .LT. N) &
|
|
& CALL DROT(N,U(1,K),1,U(1,K+1),1,CS,SN)
|
|
560 CONTINUE
|
|
E(M-1) = F
|
|
ITER = ITER + 1
|
|
GO TO 610
|
|
!
|
|
! CONVERGENCE.
|
|
!
|
|
570 CONTINUE
|
|
!
|
|
! MAKE THE SINGULAR VALUE POSITIVE.
|
|
!
|
|
IF (S(L) .GE. 0.0D0) GO TO 580
|
|
S(L) = -S(L)
|
|
if(wantv)call DSCAL(p,-1.0d0,v(1,l),1)
|
|
580 CONTINUE
|
|
!
|
|
! ORDER THE SINGULAR VALUE.
|
|
!
|
|
590 IF (L .EQ. MM) GO TO 600
|
|
! ...EXIT
|
|
IF (S(L) .GE. S(L+1)) GO TO 600
|
|
T = S(L)
|
|
S(L) = S(L+1)
|
|
S(L+1) = T
|
|
if(wantv.and.l.lt.p)call DSWAP(p,v(1,l),1,v(1,l+1),1)
|
|
if(wantu.and.l.lt.n)call DSWAP(n,u(1,l),1,u(1,l+1),1)
|
|
L = L + 1
|
|
GO TO 590
|
|
600 CONTINUE
|
|
ITER = 0
|
|
M = M - 1
|
|
610 CONTINUE
|
|
GO TO 360
|
|
620 CONTINUE
|
|
RETURN
|
|
END SUBROUTINE dsvdc
|