diff options
Diffstat (limited to 'Utilib/src/ALSVDF.f')
| -rw-r--r-- | Utilib/src/ALSVDF.f | 302 |
1 files changed, 302 insertions, 0 deletions
diff --git a/Utilib/src/ALSVDF.f b/Utilib/src/ALSVDF.f new file mode 100644 index 0000000..351424e --- /dev/null +++ b/Utilib/src/ALSVDF.f @@ -0,0 +1,302 @@ +*DECK ALSVDF + SUBROUTINE ALSVDF(A,M,N,MP,NP,W,V) +* +*----------------------------------------------------------------------- +* +*Purpose: +* singular value decomposition. +* +*Copyright: +* Copyright (C) 1993 Ecole Polytechnique de Montreal +* This library is free software; you can redistribute it and/or +* modify it under the terms of the GNU Lesser General Public +* License as published by the Free Software Foundation; either +* version 2.1 of the License, or (at your option) any later version +* +*Author(s): A. Hebert +* +*Parameters: input +* A matrix to decompose. +* M,N first/second mathematical dimension of matrix A +* MP,NP first/second physical dimension of matrix A +* +*Parameters: output +* A first decomposed matrix. +* W singular values. +* V second decomposed matrix. +* +*----------------------------------------------------------------------- +* + IMPLICIT NONE + INTEGER M,MP,N,NP + DOUBLE PRECISION A(MP,NP),V(NP,NP),W(NP) + INTEGER I,ITS,J,JJ,K,L,NM + DOUBLE PRECISION ANORM,C,F,G,H,S,SCALE,X,Y,Z,RV0 + DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: RV1,RV2 +* + ALLOCATE(RV1(NP)) + NM=0 + G=0.0D0 + SCALE=0.0D0 + ANORM=0.0D0 + DO 25 I=1,N + L=I+1 + RV1(I)=SCALE*G + G=0.0D0 + S=0.0D0 + SCALE=0.0D0 + IF(I.LE.M)THEN + DO 11 K=I,M + SCALE=SCALE+ABS(A(K,I)) +11 CONTINUE + IF(SCALE.NE.0.0D0)THEN + DO 12 K=I,M + A(K,I)=A(K,I)/SCALE + S=S+A(K,I)*A(K,I) +12 CONTINUE + F=A(I,I) + G=-SIGN(SQRT(S),F) + H=F*G-S + A(I,I)=F-G + DO 15 J=L,N + S=0.0D0 + DO 13 K=I,M + S=S+A(K,I)*A(K,J) +13 CONTINUE + F=S/H + DO 14 K=I,M + A(K,J)=A(K,J)+F*A(K,I) +14 CONTINUE +15 CONTINUE + DO 16 K=I,M + A(K,I)=SCALE*A(K,I) +16 CONTINUE + ENDIF + ENDIF + W(I)=SCALE*G + G=0.0D0 + S=0.0D0 + SCALE=0.0D0 + IF((I.LE.M).AND.(I.NE.N))THEN + DO 17 K=L,N + SCALE=SCALE+ABS(A(I,K)) +17 CONTINUE + IF(SCALE.NE.0.0D0)THEN + DO 18 K=L,N + A(I,K)=A(I,K)/SCALE + S=S+A(I,K)*A(I,K) +18 CONTINUE + F=A(I,L) + G=-SIGN(SQRT(S),F) + H=F*G-S + A(I,L)=F-G + DO 19 K=L,N + RV1(K)=A(I,K)/H +19 CONTINUE + DO 23 J=L,M + S=0.0D0 + DO 21 K=L,N + S=S+A(J,K)*A(I,K) +21 CONTINUE + DO 22 K=L,N + A(J,K)=A(J,K)+S*RV1(K) +22 CONTINUE +23 CONTINUE + DO 24 K=L,N + A(I,K)=SCALE*A(I,K) +24 CONTINUE + ENDIF + ENDIF + ANORM=MAX(ANORM,(ABS(W(I))+ABS(RV1(I)))) +25 CONTINUE + DO 32 I=N,1,-1 + IF(I.LT.N)THEN + IF(G.NE.0.0D0)THEN + DO 26 J=L,N + V(J,I)=(A(I,J)/A(I,L))/G +26 CONTINUE + DO 29 J=L,N + S=0.0D0 + DO 27 K=L,N + S=S+A(I,K)*V(K,J) +27 CONTINUE + DO 28 K=L,N + V(K,J)=V(K,J)+S*V(K,I) +28 CONTINUE +29 CONTINUE + ENDIF + DO 31 J=L,N + V(I,J)=0.0D0 + V(J,I)=0.0D0 +31 CONTINUE + ENDIF + V(I,I)=1.0D0 + G=RV1(I) + L=I +32 CONTINUE + DO 39 I=MIN(M,N),1,-1 + L=I+1 + G=W(I) + DO 33 J=L,N + A(I,J)=0.0D0 +33 CONTINUE + IF(G.NE.0.0D0)THEN + G=1.0D0/G + DO 36 J=L,N + S=0.0D0 + DO 34 K=L,M + S=S+A(K,I)*A(K,J) +34 CONTINUE + F=(S/A(I,I))*G + DO 35 K=I,M + A(K,J)=A(K,J)+F*A(K,I) +35 CONTINUE +36 CONTINUE + DO 37 J=I,M + A(J,I)=A(J,I)*G +37 CONTINUE + ELSE + DO 38 J= I,M + A(J,I)=0.0D0 +38 CONTINUE + ENDIF + A(I,I)=A(I,I)+1.0D0 +39 CONTINUE + DO 49 K=N,1,-1 + DO 48 ITS=1,30 + DO 41 L=K,1,-1 + NM=L-1 + IF((ABS(RV1(L))+ANORM).EQ.ANORM) GOTO 2 + IF((ABS(W(NM))+ANORM).EQ.ANORM) GOTO 1 +41 CONTINUE +1 C=0.0D0 + S=1.0D0 + DO 43 I=L,K + F=S*RV1(I) + RV1(I)=C*RV1(I) + IF((ABS(F)+ANORM).EQ.ANORM) GOTO 2 + G=W(I) + IF(ABS(F).GT.ABS(G))THEN + W(I)=ABS(F)*SQRT(1.0D0+(ABS(G)/ABS(F))**2) + ELSE + IF(ABS(G).EQ.0.0D0)THEN + W(I)=0.0D0 + ELSE + W(I)=ABS(G)*SQRT(1.0D0+(ABS(F)/ABS(G))**2) + ENDIF + ENDIF + H=1.0D0/W(I) + C= (G*H) + S=-(F*H) + DO 42 J=1,M + Y=A(J,NM) + Z=A(J,I) + A(J,NM)=(Y*C)+(Z*S) + A(J,I)=-(Y*S)+(Z*C) +42 CONTINUE +43 CONTINUE +2 Z=W(K) + IF(L.EQ.K)THEN + IF(Z.LT.0.0D0)THEN + W(K)=-Z + DO 44 J=1,N + V(J,K)=-V(J,K) +44 CONTINUE + ENDIF + GOTO 3 + ENDIF + IF(ITS.EQ.30) CALL XABORT('ALSVDF: NO CONVERGENCE.') + X=W(L) + NM=K-1 + Y=W(NM) + G=RV1(NM) + H=RV1(K) + F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.0D0*H*Y) + IF(ABS(F).GT.1.0D0)THEN + G=ABS(F)*SQRT(1.0D0+(1.0D0/ABS(F))**2) + ELSE + G=SQRT(1.0D0+(ABS(F))**2) + ENDIF + F=((X-Z)*(X+Z)+H*((Y/(F+SIGN(G,F)))-H))/X + C=1.0D0 + S=1.0D0 + DO 47 J=L,NM + I=J+1 + G=RV1(I) + Y=W(I) + H=S*G + G=C*G + IF(ABS(F).GT.ABS(H))THEN + Z=ABS(F)*SQRT(1.0D0+(ABS(H)/ABS(F))**2) + ELSE + IF(ABS(H).EQ.0.D0)THEN + Z=0.0D0 + ELSE + Z=ABS(H)*SQRT(1.0D0+(ABS(F)/ABS(H))**2) + ENDIF + ENDIF + RV1(J)=Z + C=F/Z + S=H/Z + F= (X*C)+(G*S) + G=-(X*S)+(G*C) + H=Y*S + Y=Y*C + DO 45 JJ=1,N + X=V(JJ,J) + Z=V(JJ,I) + V(JJ,J)= (X*C)+(Z*S) + V(JJ,I)=-(X*S)+(Z*C) +45 CONTINUE + IF(ABS(F).GT.ABS(H))THEN + Z=ABS(F)*SQRT(1.0D0+(ABS(H)/ABS(F))**2) + ELSE + IF(ABS(H).EQ.0.D0)THEN + Z=0.0D0 + ELSE + Z=ABS(H)*SQRT(1.0D0+(ABS(F)/ABS(H))**2) + ENDIF + ENDIF + W(J)=Z + IF(Z.NE.0.0D0)THEN + Z=1.0D0/Z + C=F*Z + S=H*Z + ENDIF + F= (C*G)+(S*Y) + X=-(S*G)+(C*Y) + DO 46 JJ=1,M + Y=A(JJ,J) + Z=A(JJ,I) + A(JJ,J)= (Y*C)+(Z*S) + A(JJ,I)=-(Y*S)+(Z*C) +46 CONTINUE +47 CONTINUE + RV1(L)=0.0D0 + RV1(K)=F + W(K)=X +48 CONTINUE +3 CONTINUE +49 CONTINUE + DEALLOCATE(RV1) +* +* Sort the data from highest to lowest singular value + ALLOCATE(RV1(M),RV2(N)) + DO I=1,NP + DO J=1,NP-I + IF(W(J).LE.W(J+1)) THEN + RV0=W(J) + RV1(:M)=A(:M,J) + RV2(:N)=V(:N,J) + W(J)=W(J+1) + A(:M,J)=A(:M,J+1) + V(:N,J)=V(:N,J+1) + W(J+1)=RV0 + A(:M,J+1)=RV1(:M) + V(:N,J+1)=RV2(:N) + ENDIF + ENDDO + ENDDO + DEALLOCATE(RV2,RV1) + RETURN + END |