|
|
|
|
|
SUBROUTINE BNDACC (G,MDG,NB,IP,IR,MT,JT)
|
|
|
c
|
|
|
C SEQUENTIAL ALGORITHM FOR BANDED LEAST SQUARES PROBLEM..
|
|
|
C ACCUMULATION PHASE. FOR SOLUTION PHASE USE BNDSOL.
|
|
|
c
|
|
|
c The original version of this code was developed by
|
|
|
c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
|
|
|
c 1973 JUN 12, and published in the book
|
|
|
c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
|
|
|
c Revised FEB 1995 to accompany reprinting of the book by SIAM.
|
|
|
C
|
|
|
C THE CALLING PROGRAM MUST SET IR=1 AND IP=1 BEFORE THE FIRST CALL
|
|
|
C TO BNDACC FOR A NEW CASE.
|
|
|
C
|
|
|
C THE SECOND SUBSCRIPT OF G( ) MUST BE DIMENSIONED AT LEAST
|
|
|
C NB+1 IN THE CALLING PROGRAM.
|
|
|
c ------------------------------------------------------------------
|
|
|
integer I, J, IE, IG, IG1, IG2, IP, IR, JG, JT, K, KH, L, LP1
|
|
|
integer MDG, MH, MT, MU, NB, NBP1
|
|
|
c double precision G(MDG,NB+1)
|
|
|
double precision G(MDG,*)
|
|
|
double precision RHO, ZERO
|
|
|
parameter(ZERO = 0.0d0)
|
|
|
c ------------------------------------------------------------------
|
|
|
C
|
|
|
C ALG. STEPS 1-4 ARE PERFORMED EXTERNAL TO THIS SUBROUTINE.
|
|
|
C
|
|
|
NBP1=NB+1
|
|
|
IF (MT.LE.0) RETURN
|
|
|
C ALG. STEP 5
|
|
|
IF (JT.EQ.IP) GO TO 70
|
|
|
C ALG. STEPS 6-7
|
|
|
IF (JT.LE.IR) GO TO 30
|
|
|
C ALG. STEPS 8-9
|
|
|
DO 10 I=1,MT
|
|
|
IG1=JT+MT-I
|
|
|
IG2=IR+MT-I
|
|
|
DO 10 J=1,NBP1
|
|
|
10 G(IG1,J)=G(IG2,J)
|
|
|
C ALG. STEP 10
|
|
|
IE=JT-IR
|
|
|
DO 20 I=1,IE
|
|
|
IG=IR+I-1
|
|
|
DO 20 J=1,NBP1
|
|
|
20 G(IG,J)=ZERO
|
|
|
C ALG. STEP 11
|
|
|
IR=JT
|
|
|
C ALG. STEP 12
|
|
|
30 MU=min(NB-1,IR-IP-1)
|
|
|
IF (MU.EQ.0) GO TO 60
|
|
|
C ALG. STEP 13
|
|
|
DO 50 L=1,MU
|
|
|
C ALG. STEP 14
|
|
|
K=min(L,JT-IP)
|
|
|
C ALG. STEP 15
|
|
|
LP1=L+1
|
|
|
IG=IP+L
|
|
|
DO 40 I=LP1,NB
|
|
|
JG=I-K
|
|
|
40 G(IG,JG)=G(IG,I)
|
|
|
C ALG. STEP 16
|
|
|
DO 50 I=1,K
|
|
|
JG=NBP1-I
|
|
|
50 G(IG,JG)=ZERO
|
|
|
C ALG. STEP 17
|
|
|
60 IP=JT
|
|
|
C ALG. STEPS 18-19
|
|
|
70 MH=IR+MT-IP
|
|
|
KH=min(NBP1,MH)
|
|
|
C ALG. STEP 20
|
|
|
DO 80 I=1,KH
|
|
|
CALL H12 (1,I,max(I+1,IR-IP+1),MH,G(IP,I),1,RHO,
|
|
|
* G(IP,I+1),1,MDG,NBP1-I)
|
|
|
80 continue
|
|
|
C ALG. STEP 21
|
|
|
IR=IP+KH
|
|
|
C ALG. STEP 22
|
|
|
IF (KH.LT.NBP1) GO TO 100
|
|
|
C ALG. STEP 23
|
|
|
DO 90 I=1,NB
|
|
|
90 G(IR-1,I)=ZERO
|
|
|
C ALG. STEP 24
|
|
|
100 CONTINUE
|
|
|
C ALG. STEP 25
|
|
|
RETURN
|
|
|
END
|
|
|
|
|
|
|
|
|
SUBROUTINE BNDSOL (MODE,G,MDG,NB,IP,IR,X,N,RNORM)
|
|
|
c
|
|
|
C SEQUENTIAL SOLUTION OF A BANDED LEAST SQUARES PROBLEM..
|
|
|
C SOLUTION PHASE. FOR THE ACCUMULATION PHASE USE BNDACC.
|
|
|
c
|
|
|
c The original version of this code was developed by
|
|
|
c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
|
|
|
c 1973 JUN 12, and published in the book
|
|
|
c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
|
|
|
c Revised FEB 1995 to accompany reprinting of the book by SIAM.
|
|
|
C
|
|
|
C MODE = 1 SOLVE R*X=Y WHERE R AND Y ARE IN THE G( ) ARRAY
|
|
|
C AND X WILL BE STORED IN THE X( ) ARRAY.
|
|
|
C 2 SOLVE (R**T)*X=Y WHERE R IS IN G( ),
|
|
|
C Y IS INITIALLY IN X( ), AND X REPLACES Y IN X( ),
|
|
|
C 3 SOLVE R*X=Y WHERE R IS IN G( ).
|
|
|
C Y IS INITIALLY IN X( ), AND X REPLACES Y IN X( ).
|
|
|
C
|
|
|
C THE SECOND SUBSCRIPT OF G( ) MUST BE DIMENSIONED AT LEAST
|
|
|
C NB+1 IN THE CALLING PROGRAM.
|
|
|
integer I, I1, I2, IE, II, IP, IR, IX, J, JG, L
|
|
|
integer MDG, MODE, N, NB, NP1, IRM1
|
|
|
double precision G(MDG,*), RNORM, RSQ, S, X(N), ZERO
|
|
|
parameter(ZERO = 0.0d0)
|
|
|
C
|
|
|
RNORM=ZERO
|
|
|
GO TO (10,90,50), MODE
|
|
|
C ********************* MODE = 1
|
|
|
C ALG. STEP 26
|
|
|
10 DO 20 J=1,N
|
|
|
20 X(J)=G(J,NB+1)
|
|
|
RSQ=ZERO
|
|
|
NP1=N+1
|
|
|
IRM1=IR-1
|
|
|
IF (NP1.GT.IRM1) GO TO 40
|
|
|
DO 30 J=NP1,IRM1
|
|
|
30 RSQ=RSQ+G(J,NB+1)**2
|
|
|
RNORM=SQRT(RSQ)
|
|
|
40 CONTINUE
|
|
|
C ********************* MODE = 3
|
|
|
C ALG. STEP 27
|
|
|
50 DO 80 II=1,N
|
|
|
I=N+1-II
|
|
|
C ALG. STEP 28
|
|
|
S=ZERO
|
|
|
L=max(0,I-IP)
|
|
|
C ALG. STEP 29
|
|
|
IF (I.EQ.N) GO TO 70
|
|
|
C ALG. STEP 30
|
|
|
IE=min(N+1-I,NB)
|
|
|
DO 60 J=2,IE
|
|
|
JG=J+L
|
|
|
IX=I-1+J
|
|
|
60 S=S+G(I,JG)*X(IX)
|
|
|
C ALG. STEP 31
|
|
|
70 continue
|
|
|
IF (G(I,L+1) .eq. ZERO) go to 130
|
|
|
80 X(I)=(X(I)-S)/G(I,L+1)
|
|
|
C ALG. STEP 32
|
|
|
RETURN
|
|
|
C ********************* MODE = 2
|
|
|
90 DO 120 J=1,N
|
|
|
S=ZERO
|
|
|
IF (J.EQ.1) GO TO 110
|
|
|
I1=max(1,J-NB+1)
|
|
|
I2=J-1
|
|
|
DO 100 I=I1,I2
|
|
|
L=J-I+1+max(0,I-IP)
|
|
|
100 S=S+X(I)*G(I,L)
|
|
|
110 L=max(0,J-IP)
|
|
|
IF (G(J,L+1) .eq. ZERO) go to 130
|
|
|
120 X(J)=(X(J)-S)/G(J,L+1)
|
|
|
RETURN
|
|
|
C
|
|
|
130 write (*,'(/a/a,4i6)')' ZERO DIAGONAL TERM IN BNDSOL.',
|
|
|
* ' MODE,I,J,L = ',MODE,I,J,L
|
|
|
STOP
|
|
|
END
|
|
|
|
|
|
double precision FUNCTION DIFF(X,Y)
|
|
|
c
|
|
|
c Function used in tests that depend on machine precision.
|
|
|
c
|
|
|
c The original version of this code was developed by
|
|
|
c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
|
|
|
c 1973 JUN 7, and published in the book
|
|
|
c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
|
|
|
c Revised FEB 1995 to accompany reprinting of the book by SIAM.
|
|
|
C
|
|
|
double precision X, Y
|
|
|
DIFF=X-Y
|
|
|
RETURN
|
|
|
END
|
|
|
SUBROUTINE G1 (A,B,CTERM,STERM,SIG)
|
|
|
c
|
|
|
C COMPUTE ORTHOGONAL ROTATION MATRIX..
|
|
|
c
|
|
|
c The original version of this code was developed by
|
|
|
c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
|
|
|
c 1973 JUN 12, and published in the book
|
|
|
c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
|
|
|
c Revised FEB 1995 to accompany reprinting of the book by SIAM.
|
|
|
C
|
|
|
C COMPUTE.. MATRIX (C, S) SO THAT (C, S)(A) = (SQRT(A**2+B**2))
|
|
|
C (-S,C) (-S,C)(B) ( 0 )
|
|
|
C COMPUTE SIG = SQRT(A**2+B**2)
|
|
|
C SIG IS COMPUTED LAST TO ALLOW FOR THE POSSIBILITY THAT
|
|
|
C SIG MAY BE IN THE SAME LOCATION AS A OR B .
|
|
|
C ------------------------------------------------------------------
|
|
|
double precision A, B, CTERM, ONE, SIG, STERM, XR, YR, ZERO
|
|
|
parameter(ONE = 1.0d0, ZERO = 0.0d0)
|
|
|
C ------------------------------------------------------------------
|
|
|
if (abs(A) .gt. abs(B)) then
|
|
|
XR=B/A
|
|
|
YR=sqrt(ONE+XR**2)
|
|
|
CTERM=sign(ONE/YR,A)
|
|
|
STERM=CTERM*XR
|
|
|
SIG=abs(A)*YR
|
|
|
RETURN
|
|
|
endif
|
|
|
|
|
|
if (B .ne. ZERO) then
|
|
|
XR=A/B
|
|
|
YR=sqrt(ONE+XR**2)
|
|
|
STERM=sign(ONE/YR,B)
|
|
|
CTERM=STERM*XR
|
|
|
SIG=abs(B)*YR
|
|
|
RETURN
|
|
|
endif
|
|
|
|
|
|
SIG=ZERO
|
|
|
CTERM=ZERO
|
|
|
STERM=ONE
|
|
|
RETURN
|
|
|
END
|
|
|
|
|
|
|
|
|
SUBROUTINE G2 (CTERM,STERM,X,Y)
|
|
|
c
|
|
|
C APPLY THE ROTATION COMPUTED BY G1 TO (X,Y).
|
|
|
c
|
|
|
c The original version of this code was developed by
|
|
|
c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
|
|
|
c 1972 DEC 15, and published in the book
|
|
|
c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
|
|
|
c Revised FEB 1995 to accompany reprinting of the book by SIAM.
|
|
|
c ------------------------------------------------------------------
|
|
|
double precision CTERM, STERM, X, XR, Y
|
|
|
c ------------------------------------------------------------------
|
|
|
XR=CTERM*X+STERM*Y
|
|
|
Y=-STERM*X+CTERM*Y
|
|
|
X=XR
|
|
|
RETURN
|
|
|
END
|
|
|
|
|
|
double precision FUNCTION GEN2(ANOISE)
|
|
|
c
|
|
|
C GENERATE NUMBERS FOR CONSTRUCTION OF TEST CASES.
|
|
|
c
|
|
|
c The original version of this code was developed by
|
|
|
c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
|
|
|
c 1972 DEC 15, and published in the book
|
|
|
c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
|
|
|
c Revised FEB 1995 to accompany reprinting of the book by SIAM.
|
|
|
c ------------------------------------------------------------------
|
|
|
integer I, J, MI, MJ
|
|
|
double precision AI, AJ, ANOISE, ZERO
|
|
|
parameter(ZERO = 0.0d0)
|
|
|
SAVE
|
|
|
c ------------------------------------------------------------------
|
|
|
IF (ANOISE) 10,30,20
|
|
|
10 MI=891
|
|
|
MJ=457
|
|
|
I=5
|
|
|
J=7
|
|
|
AJ= ZERO
|
|
|
GEN2= ZERO
|
|
|
RETURN
|
|
|
C
|
|
|
C THE SEQUENCE OF VALUES OF J IS BOUNDED BETWEEN 1 AND 996
|
|
|
C IF INITIAL J = 1,2,3,4,5,6,7,8, OR 9, THE PERIOD IS 332
|
|
|
20 J=J*MJ
|
|
|
J=J-997*(J/997)
|
|
|
AJ=J-498
|
|
|
C THE SEQUENCE OF VALUES OF I IS BOUNDED BETWEEN 1 AND 999
|
|
|
C IF INITIAL I = 1,2,3,6,7, OR 9, THE PERIOD WILL BE 50
|
|
|
C IF INITIAL I = 4 OR 8 THE PERIOD WILL BE 25
|
|
|
C IF INITIAL I = 5 THE PERIOD WILL BE 10
|
|
|
30 I=I*MI
|
|
|
I=I-1000*(I/1000)
|
|
|
AI=I-500
|
|
|
GEN2=AI+AJ*ANOISE
|
|
|
RETURN
|
|
|
END
|
|
|
|
|
|
|
|
|
C SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV)
|
|
|
C
|
|
|
C CONSTRUCTION AND/OR APPLICATION OF A SINGLE
|
|
|
C HOUSEHOLDER TRANSFORMATION.. Q = I + U*(U**T)/B
|
|
|
C
|
|
|
c The original version of this code was developed by
|
|
|
c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
|
|
|
c 1973 JUN 12, and published in the book
|
|
|
c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
|
|
|
c Revised FEB 1995 to accompany reprinting of the book by SIAM.
|
|
|
C ------------------------------------------------------------------
|
|
|
c Subroutine Arguments
|
|
|
c
|
|
|
C MODE = 1 OR 2 Selects Algorithm H1 to construct and apply a
|
|
|
c Householder transformation, or Algorithm H2 to apply a
|
|
|
c previously constructed transformation.
|
|
|
C LPIVOT IS THE INDEX OF THE PIVOT ELEMENT.
|
|
|
C L1,M IF L1 .LE. M THE TRANSFORMATION WILL BE CONSTRUCTED TO
|
|
|
C ZERO ELEMENTS INDEXED FROM L1 THROUGH M. IF L1 GT. M
|
|
|
C THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION.
|
|
|
C U(),IUE,UP On entry with MODE = 1, U() contains the pivot
|
|
|
c vector. IUE is the storage increment between elements.
|
|
|
c On exit when MODE = 1, U() and UP contain quantities
|
|
|
c defining the vector U of the Householder transformation.
|
|
|
c on entry with MODE = 2, U() and UP should contain
|
|
|
c quantities previously computed with MODE = 1. These will
|
|
|
c not be modified during the entry with MODE = 2.
|
|
|
C C() ON ENTRY with MODE = 1 or 2, C() CONTAINS A MATRIX WHICH
|
|
|
c WILL BE REGARDED AS A SET OF VECTORS TO WHICH THE
|
|
|
c HOUSEHOLDER TRANSFORMATION IS TO BE APPLIED.
|
|
|
c ON EXIT C() CONTAINS THE SET OF TRANSFORMED VECTORS.
|
|
|
C ICE STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C().
|
|
|
C ICV STORAGE INCREMENT BETWEEN VECTORS IN C().
|
|
|
C NCV NUMBER OF VECTORS IN C() TO BE TRANSFORMED. IF NCV .LE. 0
|
|
|
C NO OPERATIONS WILL BE DONE ON C().
|
|
|
C ------------------------------------------------------------------
|
|
|
SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV)
|
|
|
C ------------------------------------------------------------------
|
|
|
integer I, I2, I3, I4, ICE, ICV, INCR, IUE, J
|
|
|
integer L1, LPIVOT, M, MODE, NCV
|
|
|
double precision B, C(*), CL, CLINV, ONE, SM
|
|
|
c double precision U(IUE,M)
|
|
|
double precision U(IUE,*)
|
|
|
double precision UP
|
|
|
parameter(ONE = 1.0d0)
|
|
|
C ------------------------------------------------------------------
|
|
|
IF (0.GE.LPIVOT.OR.LPIVOT.GE.L1.OR.L1.GT.M) RETURN
|
|
|
CL=abs(U(1,LPIVOT))
|
|
|
IF (MODE.EQ.2) GO TO 60
|
|
|
C ****** CONSTRUCT THE TRANSFORMATION. ******
|
|
|
DO 10 J=L1,M
|
|
|
10 CL=MAX(abs(U(1,J)),CL)
|
|
|
IF (CL) 130,130,20
|
|
|
20 CLINV=ONE/CL
|
|
|
SM=(U(1,LPIVOT)*CLINV)**2
|
|
|
DO 30 J=L1,M
|
|
|
30 SM=SM+(U(1,J)*CLINV)**2
|
|
|
CL=CL*SQRT(SM)
|
|
|
IF (U(1,LPIVOT)) 50,50,40
|
|
|
40 CL=-CL
|
|
|
50 UP=U(1,LPIVOT)-CL
|
|
|
U(1,LPIVOT)=CL
|
|
|
GO TO 70
|
|
|
C ****** APPLY THE TRANSFORMATION I+U*(U**T)/B TO C. ******
|
|
|
C
|
|
|
60 IF (CL) 130,130,70
|
|
|
70 IF (NCV.LE.0) RETURN
|
|
|
B= UP*U(1,LPIVOT)
|
|
|
C B MUST BE NONPOSITIVE HERE. IF B = 0., RETURN.
|
|
|
C
|
|
|
IF (B) 80,130,130
|
|
|
80 B=ONE/B
|
|
|
I2=1-ICV+ICE*(LPIVOT-1)
|
|
|
INCR=ICE*(L1-LPIVOT)
|
|
|
DO 120 J=1,NCV
|
|
|
I2=I2+ICV
|
|
|
I3=I2+INCR
|
|
|
I4=I3
|
|
|
SM=C(I2)*UP
|
|
|
DO 90 I=L1,M
|
|
|
SM=SM+C(I3)*U(1,I)
|
|
|
90 I3=I3+ICE
|
|
|
IF (SM) 100,120,100
|
|
|
100 SM=SM*B
|
|
|
C(I2)=C(I2)+SM*UP
|
|
|
DO 110 I=L1,M
|
|
|
C(I4)=C(I4)+SM*U(1,I)
|
|
|
110 I4=I4+ICE
|
|
|
120 CONTINUE
|
|
|
130 RETURN
|
|
|
END
|
|
|
|
|
|
SUBROUTINE HFTI (A,MDA,M,N,B,MDB,NB,TAU,KRANK,RNORM,H,G,IP)
|
|
|
c
|
|
|
C SOLVE LEAST SQUARES PROBLEM USING ALGORITHM, HFTI.
|
|
|
c Householder Forward Triangulation with column Interchanges.
|
|
|
c
|
|
|
c The original version of this code was developed by
|
|
|
c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
|
|
|
c 1973 JUN 12, and published in the book
|
|
|
c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
|
|
|
c Revised FEB 1995 to accompany reprinting of the book by SIAM.
|
|
|
C ------------------------------------------------------------------
|
|
|
integer I, II, IP1, J, JB, JJ, K, KP1, KRANK
|
|
|
integer L, LDIAG, LMAX, M, MDA, MDB, N, NB
|
|
|
c integer IP(N)
|
|
|
c double precision A(MDA,N),B(MDB,NB),H(N),G(N),RNORM(NB)
|
|
|
integer IP(*)
|
|
|
double precision A(MDA,*),B(MDB, *),H(*),G(*),RNORM( *)
|
|
|
double precision DIFF, FACTOR, HMAX, SM, TAU, TMP, ZERO
|
|
|
parameter(FACTOR = 0.001d0, ZERO = 0.0d0)
|
|
|
C ------------------------------------------------------------------
|
|
|
C
|
|
|
K=0
|
|
|
LDIAG=min(M,N)
|
|
|
IF (LDIAG.LE.0) GO TO 270
|
|
|
DO 80 J=1,LDIAG
|
|
|
IF (J.EQ.1) GO TO 20
|
|
|
C
|
|
|
C UPDATE SQUARED COLUMN LENGTHS AND FIND LMAX
|
|
|
C ..
|
|
|
LMAX=J
|
|
|
DO 10 L=J,N
|
|
|
H(L)=H(L)-A(J-1,L)**2
|
|
|
IF (H(L).GT.H(LMAX)) LMAX=L
|
|
|
10 CONTINUE
|
|
|
IF(DIFF(HMAX+FACTOR*H(LMAX),HMAX)) 20,20,50
|
|
|
C
|
|
|
C COMPUTE SQUARED COLUMN LENGTHS AND FIND LMAX
|
|
|
C ..
|
|
|
20 LMAX=J
|
|
|
DO 40 L=J,N
|
|
|
H(L)=0.
|
|
|
DO 30 I=J,M
|
|
|
30 H(L)=H(L)+A(I,L)**2
|
|
|
IF (H(L).GT.H(LMAX)) LMAX=L
|
|
|
40 CONTINUE
|
|
|
HMAX=H(LMAX)
|
|
|
C ..
|
|
|
C LMAX HAS BEEN DETERMINED
|
|
|
C
|
|
|
C DO COLUMN INTERCHANGES IF NEEDED.
|
|
|
C ..
|
|
|
50 CONTINUE
|
|
|
IP(J)=LMAX
|
|
|
IF (IP(J).EQ.J) GO TO 70
|
|
|
DO 60 I=1,M
|
|
|
TMP=A(I,J)
|
|
|
A(I,J)=A(I,LMAX)
|
|
|
60 A(I,LMAX)=TMP
|
|
|
H(LMAX)=H(J)
|
|
|
C
|
|
|
C COMPUTE THE J-TH TRANSFORMATION AND APPLY IT TO A AND B.
|
|
|
C ..
|
|
|
70 CALL H12 (1,J,J+1,M,A(1,J),1,H(J),A(1,J+1),1,MDA,N-J)
|
|
|
80 CALL H12 (2,J,J+1,M,A(1,J),1,H(J),B,1,MDB,NB)
|
|
|
C
|
|
|
C DETERMINE THE PSEUDORANK, K, USING THE TOLERANCE, TAU.
|
|
|
C ..
|
|
|
DO 90 J=1,LDIAG
|
|
|
IF (ABS(A(J,J)).LE.TAU) GO TO 100
|
|
|
90 CONTINUE
|
|
|
K=LDIAG
|
|
|
GO TO 110
|
|
|
100 K=J-1
|
|
|
110 KP1=K+1
|
|
|
C
|
|
|
C COMPUTE THE NORMS OF THE RESIDUAL VECTORS.
|
|
|
C
|
|
|
IF (NB.LE.0) GO TO 140
|
|
|
DO 130 JB=1,NB
|
|
|
TMP=ZERO
|
|
|
IF (KP1.GT.M) GO TO 130
|
|
|
DO 120 I=KP1,M
|
|
|
120 TMP=TMP+B(I,JB)**2
|
|
|
130 RNORM(JB)=SQRT(TMP)
|
|
|
140 CONTINUE
|
|
|
C SPECIAL FOR PSEUDORANK = 0
|
|
|
IF (K.GT.0) GO TO 160
|
|
|
IF (NB.LE.0) GO TO 270
|
|
|
DO 150 JB=1,NB
|
|
|
DO 150 I=1,N
|
|
|
150 B(I,JB)=ZERO
|
|
|
GO TO 270
|
|
|
C
|
|
|
C IF THE PSEUDORANK IS LESS THAN N COMPUTE HOUSEHOLDER
|
|
|
C DECOMPOSITION OF FIRST K ROWS.
|
|
|
C ..
|
|
|
160 IF (K.EQ.N) GO TO 180
|
|
|
DO 170 II=1,K
|
|
|
I=KP1-II
|
|
|
170 CALL H12 (1,I,KP1,N,A(I,1),MDA,G(I),A,MDA,1,I-1)
|
|
|
180 CONTINUE
|
|
|
C
|
|
|
C
|
|
|
IF (NB.LE.0) GO TO 270
|
|
|
DO 260 JB=1,NB
|
|
|
C
|
|
|
C SOLVE THE K BY K TRIANGULAR SYSTEM.
|
|
|
C ..
|
|
|
DO 210 L=1,K
|
|
|
SM=ZERO
|
|
|
I=KP1-L
|
|
|
IF (I.EQ.K) GO TO 200
|
|
|
IP1=I+1
|
|
|
DO 190 J=IP1,K
|
|
|
190 SM=SM+A(I,J)*B(J,JB)
|
|
|
200 continue
|
|
|
210 B(I,JB)=(B(I,JB)-SM)/A(I,I)
|
|
|
C
|
|
|
C COMPLETE COMPUTATION OF SOLUTION VECTOR.
|
|
|
C ..
|
|
|
IF (K.EQ.N) GO TO 240
|
|
|
DO 220 J=KP1,N
|
|
|
220 B(J,JB)=ZERO
|
|
|
DO 230 I=1,K
|
|
|
230 CALL H12 (2,I,KP1,N,A(I,1),MDA,G(I),B(1,JB),1,MDB,1)
|
|
|
C
|
|
|
C RE-ORDER THE SOLUTION VECTOR TO COMPENSATE FOR THE
|
|
|
C COLUMN INTERCHANGES.
|
|
|
C ..
|
|
|
240 DO 250 JJ=1,LDIAG
|
|
|
J=LDIAG+1-JJ
|
|
|
IF (IP(J).EQ.J) GO TO 250
|
|
|
L=IP(J)
|
|
|
TMP=B(L,JB)
|
|
|
B(L,JB)=B(J,JB)
|
|
|
B(J,JB)=TMP
|
|
|
250 CONTINUE
|
|
|
260 CONTINUE
|
|
|
C ..
|
|
|
C THE SOLUTION VECTORS, X, ARE NOW
|
|
|
C IN THE FIRST N ROWS OF THE ARRAY B(,).
|
|
|
C
|
|
|
270 KRANK=K
|
|
|
RETURN
|
|
|
END
|
|
|
|
|
|
SUBROUTINE LDP (G,MDG,M,N,H,X,XNORM,W,INDEX,MODE)
|
|
|
c
|
|
|
C Algorithm LDP: LEAST DISTANCE PROGRAMMING
|
|
|
c
|
|
|
c The original version of this code was developed by
|
|
|
c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
|
|
|
c 1974 MAR 1, and published in the book
|
|
|
c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
|
|
|
c Revised FEB 1995 to accompany reprinting of the book by SIAM.
|
|
|
C ------------------------------------------------------------------
|
|
|
integer I, IW, IWDUAL, IY, IZ, J, JF, M, MDG, MODE, N, NP1
|
|
|
c integer INDEX(M)
|
|
|
c double precision G(MDG,N), H(M), X(N), W(*)
|
|
|
integer INDEX(*)
|
|
|
double precision G(MDG,*), H(*), X(*), W(*)
|
|
|
double precision DIFF, FAC, ONE, RNORM, XNORM, ZERO
|
|
|
parameter(ONE = 1.0d0, ZERO = 0.0d0)
|
|
|
C ------------------------------------------------------------------
|
|
|
IF (N.LE.0) GO TO 120
|
|
|
DO 10 J=1,N
|
|
|
10 X(J)=ZERO
|
|
|
XNORM=ZERO
|
|
|
IF (M.LE.0) GO TO 110
|
|
|
C
|
|
|
C THE DECLARED DIMENSION OF W() MUST BE AT LEAST (N+1)*(M+2)+2*M.
|
|
|
C
|
|
|
C FIRST (N+1)*M LOCS OF W() = MATRIX E FOR PROBLEM NNLS.
|
|
|
C NEXT N+1 LOCS OF W() = VECTOR F FOR PROBLEM NNLS.
|
|
|
C NEXT N+1 LOCS OF W() = VECTOR Z FOR PROBLEM NNLS.
|
|
|
C NEXT M LOCS OF W() = VECTOR Y FOR PROBLEM NNLS.
|
|
|
C NEXT M LOCS OF W() = VECTOR WDUAL FOR PROBLEM NNLS.
|
|
|
C COPY G**T INTO FIRST N ROWS AND M COLUMNS OF E.
|
|
|
C COPY H**T INTO ROW N+1 OF E.
|
|
|
C
|
|
|
IW=0
|
|
|
DO 30 J=1,M
|
|
|
DO 20 I=1,N
|
|
|
IW=IW+1
|
|
|
20 W(IW)=G(J,I)
|
|
|
IW=IW+1
|
|
|
30 W(IW)=H(J)
|
|
|
JF=IW+1
|
|
|
C STORE N ZEROS FOLLOWED BY A ONE INTO F.
|
|
|
DO 40 I=1,N
|
|
|
IW=IW+1
|
|
|
40 W(IW)=ZERO
|
|
|
W(IW+1)=ONE
|
|
|
C
|
|
|
NP1=N+1
|
|
|
IZ=IW+2
|
|
|
IY=IZ+NP1
|
|
|
IWDUAL=IY+M
|
|
|
C
|
|
|
CALL NNLS (W,NP1,NP1,M,W(JF),W(IY),RNORM,W(IWDUAL),W(IZ),INDEX,
|
|
|
* MODE)
|
|
|
C USE THE FOLLOWING RETURN IF UNSUCCESSFUL IN NNLS.
|
|
|
IF (MODE.NE.1) RETURN
|
|
|
IF (RNORM) 130,130,50
|
|
|
50 FAC=ONE
|
|
|
IW=IY-1
|
|
|
DO 60 I=1,M
|
|
|
IW=IW+1
|
|
|
C HERE WE ARE USING THE SOLUTION VECTOR Y.
|
|
|
60 FAC=FAC-H(I)*W(IW)
|
|
|
C
|
|
|
IF (DIFF(ONE+FAC,ONE)) 130,130,70
|
|
|
70 FAC=ONE/FAC
|
|
|
DO 90 J=1,N
|
|
|
IW=IY-1
|
|
|
DO 80 I=1,M
|
|
|
IW=IW+1
|
|
|
C HERE WE ARE USING THE SOLUTION VECTOR Y.
|
|
|
80 X(J)=X(J)+G(I,J)*W(IW)
|
|
|
90 X(J)=X(J)*FAC
|
|
|
DO 100 J=1,N
|
|
|
100 XNORM=XNORM+X(J)**2
|
|
|
XNORM=sqrt(XNORM)
|
|
|
C SUCCESSFUL RETURN.
|
|
|
110 MODE=1
|
|
|
RETURN
|
|
|
C ERROR RETURN. N .LE. 0.
|
|
|
120 MODE=2
|
|
|
RETURN
|
|
|
C RETURNING WITH CONSTRAINTS NOT COMPATIBLE.
|
|
|
130 MODE=4
|
|
|
RETURN
|
|
|
END
|
|
|
|
|
|
subroutine MFEOUT (A, MDA, M, N, NAMES, MODE, UNIT, WIDTH)
|
|
|
c
|
|
|
C LABELED MATRIX OUTPUT FOR USE WITH SINGULAR VALUE ANALYSIS.
|
|
|
C
|
|
|
c The original version of this code was developed by
|
|
|
c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
|
|
|
c 1973 JUN 12, and published in the book
|
|
|
c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
|
|
|
c Revised FEB 1995 to accompany reprinting of the book by SIAM.
|
|
|
C ------------------------------------------------------------------
|
|
|
c This 1995 version has additional arguments, UNIT and WIDTH,
|
|
|
c to support user options regarding the output unit and the width of
|
|
|
c print lines. Also allows user to choose length of names in NAMES().
|
|
|
c ------------------------------------------------------------------
|
|
|
c Subroutine Arguments
|
|
|
C All are input arguments. None are modified by this subroutine.
|
|
|
c
|
|
|
C A(,) Array containing matrix to be output.
|
|
|
C MDA First dimension of the array, A(,).
|
|
|
C M, N No. of rows and columns, respectively in the matrix
|
|
|
c contained in A(,).
|
|
|
C NAMES() [character array] Array of names.
|
|
|
c If NAMES(1) contains only blanks, the rest of the NAMES()
|
|
|
c array will be ignored.
|
|
|
C MODE =1 Write header for V matrix and use an F format.
|
|
|
C =2 Write header for for candidate solutions and use
|
|
|
c P format.
|
|
|
c UNIT [integer] Selects output unit. If UNIT .ge. 0 then UNIT
|
|
|
c is the output unit number. If UNIT = -1, output to
|
|
|
c the '*' unit.
|
|
|
c WIDTH [integer] Selects width of output lines.
|
|
|
c Each output line from this subroutine will have at most
|
|
|
c max(26,min(124,WIDTH)) characters plus one additional
|
|
|
c leading character for Fortran "carriage control". The
|
|
|
c carriage control character will always be a blank.
|
|
|
c ------------------------------------------------------------------
|
|
|
integer I, J, J1, J2, KBLOCK, L, LENNAM
|
|
|
integer M, MAXCOL, MDA, MODE, N, NAMSIZ, NBLOCK, UNIT, WIDTH
|
|
|
double precision A(MDA,N)
|
|
|
character NAMES(M)*(*)
|
|
|
character*4 HEAD (2)
|
|
|
character*26 FMT1(2)
|
|
|
character*26 FMT2(2)
|
|
|
logical BLKNAM, STAR
|
|
|
C
|
|
|
data HEAD(1)/' COL'/
|
|
|
data HEAD(2)/'SOLN'/
|
|
|
data FMT1 / '(/7x,00x,8(5x,a4,i4,1x)/)',
|
|
|
* '(/7x,00x,8(2x,a4,i4,4x)/)'/
|
|
|
data FMT2 / '(1x,i4,1x,a00,1x,4p8f14.0)',
|
|
|
* '(1x,i4,1x,a00,1x,8g14.6 )'/
|
|
|
c ------------------------------------------------------------------
|
|
|
if (M .le. 0 .or. N .le. 0) return
|
|
|
STAR = UNIT .lt. 0
|
|
|
C
|
|
|
c The LEN function returns the char length of a single element of
|
|
|
c the NAMES() array.
|
|
|
c
|
|
|
LENNAM = len(NAMES(1))
|
|
|
BLKNAM = NAMES(1) .eq. ' '
|
|
|
NAMSIZ = 1
|
|
|
if(.not. BLKNAM) then
|
|
|
do 30 I = 1,M
|
|
|
do 10 L = LENNAM, NAMSIZ+1, -1
|
|
|
if(NAMES(I)(L:L) .ne. ' ') then
|
|
|
NAMSIZ = L
|
|
|
go to 20
|
|
|
endif
|
|
|
10 continue
|
|
|
20 continue
|
|
|
30 continue
|
|
|
endif
|
|
|
c
|
|
|
write(FMT1(MODE)(6:7),'(i2.2)') NAMSIZ
|
|
|
write(FMT2(MODE)(12:13),'(i2.2)') NAMSIZ
|
|
|
|
|
|
70 format(/' V-Matrix of the Singular Value Decomposition of A*D.'/
|
|
|
* ' (Elements of V scaled up by a factor of 10**4)')
|
|
|
80 format(/' Sequence of candidate solutions, X')
|
|
|
|
|
|
if(STAR) then
|
|
|
if (MODE .eq. 1) then
|
|
|
write (*,70)
|
|
|
else
|
|
|
write (*,80)
|
|
|
endif
|
|
|
else
|
|
|
if (MODE .eq. 1) then
|
|
|
write (UNIT,70)
|
|
|
else
|
|
|
write (UNIT,80)
|
|
|
endif
|
|
|
endif
|
|
|
c
|
|
|
c With NAMSIZ characters allowed for the "name" and MAXCOL
|
|
|
c columns of numbers, the total line width, exclusive of a
|
|
|
c carriage control character, will be 6 + LENNAM + 14 * MAXCOL.
|
|
|
c
|
|
|
MAXCOL = max(1,min(8,(WIDTH - 6 - NAMSIZ)/14))
|
|
|
C
|
|
|
NBLOCK = (N + MAXCOL -1) / MAXCOL
|
|
|
J2 = 0
|
|
|
C
|
|
|
do 50 KBLOCK = 1, NBLOCK
|
|
|
J1 = J2 + 1
|
|
|
J2 = min(N, J2 + MAXCOL)
|
|
|
if(STAR) then
|
|
|
write (*,FMT1(MODE)) (HEAD(MODE),J,J=J1,J2)
|
|
|
else
|
|
|
write (UNIT,FMT1(MODE)) (HEAD(MODE),J,J=J1,J2)
|
|
|
endif
|
|
|
C
|
|
|
do 40 I=1,M
|
|
|
if(STAR) then
|
|
|
if(BLKNAM) then
|
|
|
write (*,FMT2(MODE)) I,' ',(A(I,J),J=J1,J2)
|
|
|
else
|
|
|
write (*,FMT2(MODE)) I,NAMES(I),(A(I,J),J=J1,J2)
|
|
|
endif
|
|
|
else
|
|
|
if(BLKNAM) then
|
|
|
write (UNIT,FMT2(MODE)) I,' ',(A(I,J),J=J1,J2)
|
|
|
else
|
|
|
write (UNIT,FMT2(MODE)) I,NAMES(I),(A(I,J),J=J1,J2)
|
|
|
endif
|
|
|
endif
|
|
|
40 continue
|
|
|
50 continue
|
|
|
end
|
|
|
|
|
|
C SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE)
|
|
|
C
|
|
|
C Algorithm NNLS: NONNEGATIVE LEAST SQUARES
|
|
|
C
|
|
|
c The original version of this code was developed by
|
|
|
c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
|
|
|
c 1973 JUN 15, and published in the book
|
|
|
c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
|
|
|
c Revised FEB 1995 to accompany reprinting of the book by SIAM.
|
|
|
c
|
|
|
C GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B, COMPUTE AN
|
|
|
C N-VECTOR, X, THAT SOLVES THE LEAST SQUARES PROBLEM
|
|
|
C
|
|
|
C A * X = B SUBJECT TO X .GE. 0
|
|
|
C ------------------------------------------------------------------
|
|
|
c Subroutine Arguments
|
|
|
c
|
|
|
C A(),MDA,M,N MDA IS THE FIRST DIMENSIONING PARAMETER FOR THE
|
|
|
C ARRAY, A(). ON ENTRY A() CONTAINS THE M BY N
|
|
|
C MATRIX, A. ON EXIT A() CONTAINS
|
|
|
C THE PRODUCT MATRIX, Q*A , WHERE Q IS AN
|
|
|
C M BY M ORTHOGONAL MATRIX GENERATED IMPLICITLY BY
|
|
|
C THIS SUBROUTINE.
|
|
|
C B() ON ENTRY B() CONTAINS THE M-VECTOR, B. ON EXIT B() CON-
|
|
|
C TAINS Q*B.
|
|
|
C X() ON ENTRY X() NEED NOT BE INITIALIZED. ON EXIT X() WILL
|
|
|
C CONTAIN THE SOLUTION VECTOR.
|
|
|
C RNORM ON EXIT RNORM CONTAINS THE EUCLIDEAN NORM OF THE
|
|
|
C RESIDUAL VECTOR.
|
|
|
C W() AN N-ARRAY OF WORKING SPACE. ON EXIT W() WILL CONTAIN
|
|
|
C THE DUAL SOLUTION VECTOR. W WILL SATISFY W(I) = 0.
|
|
|
C FOR ALL I IN SET P AND W(I) .LE. 0. FOR ALL I IN SET Z
|
|
|
C ZZ() AN M-ARRAY OF WORKING SPACE.
|
|
|
C INDEX() AN INTEGER WORKING ARRAY OF LENGTH AT LEAST N.
|
|
|
C ON EXIT THE CONTENTS OF THIS ARRAY DEFINE THE SETS
|
|
|
C P AND Z AS FOLLOWS..
|
|
|
C
|
|
|
C INDEX(1) THRU INDEX(NSETP) = SET P.
|
|
|
C INDEX(IZ1) THRU INDEX(IZ2) = SET Z.
|
|
|
C IZ1 = NSETP + 1 = NPP1
|
|
|
C IZ2 = N
|
|
|
C MODE THIS IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING
|
|
|
C MEANINGS.
|
|
|
C 1 THE SOLUTION HAS BEEN COMPUTED SUCCESSFULLY.
|
|
|
C 2 THE DIMENSIONS OF THE PROBLEM ARE BAD.
|
|
|
C EITHER M .LE. 0 OR N .LE. 0.
|
|
|
C 3 ITERATION COUNT EXCEEDED. MORE THAN 3*N ITERATIONS.
|
|
|
C
|
|
|
C ------------------------------------------------------------------
|
|
|
SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE)
|
|
|
C ------------------------------------------------------------------
|
|
|
integer I, II, IP, ITER, ITMAX, IZ, IZ1, IZ2, IZMAX, J, JJ, JZ, L
|
|
|
integer M, MDA, MODE,N, NPP1, NSETP, RTNKEY
|
|
|
c integer INDEX(N)
|
|
|
c double precision A(MDA,N), B(M), W(N), X(N), ZZ(M)
|
|
|
integer INDEX(*)
|
|
|
double precision A(MDA,*), B(*), W(*), X(*), ZZ(*)
|
|
|
double precision ALPHA, ASAVE, CC, DIFF, DUMMY, FACTOR, RNORM
|
|
|
double precision SM, SS, T, TEMP, TWO, UNORM, UP, WMAX
|
|
|
double precision ZERO, ZTEST
|
|
|
parameter(FACTOR = 0.01d0)
|
|
|
parameter(TWO = 2.0d0, ZERO = 0.0d0)
|
|
|
C ------------------------------------------------------------------
|
|
|
MODE=1
|
|
|
IF (M .le. 0 .or. N .le. 0) then
|
|
|
MODE=2
|
|
|
RETURN
|
|
|
endif
|
|
|
ITER=0
|
|
|
ITMAX=3*N
|
|
|
C
|
|
|
C INITIALIZE THE ARRAYS INDEX() AND X().
|
|
|
C
|
|
|
DO 20 I=1,N
|
|
|
X(I)=ZERO
|
|
|
20 INDEX(I)=I
|
|
|
C
|
|
|
IZ2=N
|
|
|
IZ1=1
|
|
|
NSETP=0
|
|
|
NPP1=1
|
|
|
C ****** MAIN LOOP BEGINS HERE ******
|
|
|
30 CONTINUE
|
|
|
C QUIT IF ALL COEFFICIENTS ARE ALREADY IN THE SOLUTION.
|
|
|
C OR IF M COLS OF A HAVE BEEN TRIANGULARIZED.
|
|
|
C
|
|
|
IF (IZ1 .GT.IZ2.OR.NSETP.GE.M) GO TO 350
|
|
|
C
|
|
|
C COMPUTE COMPONENTS OF THE DUAL (NEGATIVE GRADIENT) VECTOR W().
|
|
|
C
|
|
|
DO 50 IZ=IZ1,IZ2
|
|
|
J=INDEX(IZ)
|
|
|
SM=ZERO
|
|
|
DO 40 L=NPP1,M
|
|
|
40 SM=SM+A(L,J)*B(L)
|
|
|
W(J)=SM
|
|
|
50 continue
|
|
|
C FIND LARGEST POSITIVE W(J).
|
|
|
60 continue
|
|
|
WMAX=ZERO
|
|
|
DO 70 IZ=IZ1,IZ2
|
|
|
J=INDEX(IZ)
|
|
|
IF (W(J) .gt. WMAX) then
|
|
|
WMAX=W(J)
|
|
|
IZMAX=IZ
|
|
|
endif
|
|
|
70 CONTINUE
|
|
|
C
|
|
|
C IF WMAX .LE. 0. GO TO TERMINATION.
|
|
|
C THIS INDICATES SATISFACTION OF THE KUHN-TUCKER CONDITIONS.
|
|
|
C
|
|
|
IF (WMAX .le. ZERO) go to 350
|
|
|
IZ=IZMAX
|
|
|
J=INDEX(IZ)
|
|
|
C
|
|
|
C THE SIGN OF W(J) IS OK FOR J TO BE MOVED TO SET P.
|
|
|
C BEGIN THE TRANSFORMATION AND CHECK NEW DIAGONAL ELEMENT TO AVOID
|
|
|
C NEAR LINEAR DEPENDENCE.
|
|
|
C
|
|
|
ASAVE=A(NPP1,J)
|
|
|
CALL H12 (1,NPP1,NPP1+1,M,A(1,J),1,UP,DUMMY,1,1,0)
|
|
|
UNORM=ZERO
|
|
|
IF (NSETP .ne. 0) then
|
|
|
DO 90 L=1,NSETP
|
|
|
90 UNORM=UNORM+A(L,J)**2
|
|
|
endif
|
|
|
UNORM=sqrt(UNORM)
|
|
|
IF (DIFF(UNORM+ABS(A(NPP1,J))*FACTOR,UNORM) .gt. ZERO) then
|
|
|
C
|
|
|
C COL J IS SUFFICIENTLY INDEPENDENT. COPY B INTO ZZ, UPDATE ZZ
|
|
|
C AND SOLVE FOR ZTEST ( = PROPOSED NEW VALUE FOR X(J) ).
|
|
|
C
|
|
|
DO 120 L=1,M
|
|
|
120 ZZ(L)=B(L)
|
|
|
CALL H12 (2,NPP1,NPP1+1,M,A(1,J),1,UP,ZZ,1,1,1)
|
|
|
ZTEST=ZZ(NPP1)/A(NPP1,J)
|
|
|
C
|
|
|
C SEE IF ZTEST IS POSITIVE
|
|
|
C
|
|
|
IF (ZTEST .gt. ZERO) go to 140
|
|
|
endif
|
|
|
C
|
|
|
C REJECT J AS A CANDIDATE TO BE MOVED FROM SET Z TO SET P.
|
|
|
C RESTORE A(NPP1,J), SET W(J)=0., AND LOOP BACK TO TEST DUAL
|
|
|
C COEFFS AGAIN.
|
|
|
C
|
|
|
A(NPP1,J)=ASAVE
|
|
|
W(J)=ZERO
|
|
|
GO TO 60
|
|
|
C
|
|
|
C THE INDEX J=INDEX(IZ) HAS BEEN SELECTED TO BE MOVED FROM
|
|
|
C SET Z TO SET P. UPDATE B, UPDATE INDICES, APPLY HOUSEHOLDER
|
|
|
C TRANSFORMATIONS TO COLS IN NEW SET Z, ZERO SUBDIAGONAL ELTS IN
|
|
|
C COL J, SET W(J)=0.
|
|
|
C
|
|
|
140 continue
|
|
|
DO 150 L=1,M
|
|
|
150 B(L)=ZZ(L)
|
|
|
C
|
|
|
INDEX(IZ)=INDEX(IZ1)
|
|
|
INDEX(IZ1)=J
|
|
|
IZ1=IZ1+1
|
|
|
NSETP=NPP1
|
|
|
NPP1=NPP1+1
|
|
|
C
|
|
|
IF (IZ1 .le. IZ2) then
|
|
|
DO 160 JZ=IZ1,IZ2
|
|
|
JJ=INDEX(JZ)
|
|
|
CALL H12 (2,NSETP,NPP1,M,A(1,J),1,UP,A(1,JJ),1,MDA,1)
|
|
|
160 continue
|
|
|
endif
|
|
|
C
|
|
|
IF (NSETP .ne. M) then
|
|
|
DO 180 L=NPP1,M
|
|
|
180 A(L,J)=ZERO
|
|
|
endif
|
|
|
C
|
|
|
W(J)=ZERO
|
|
|
C SOLVE THE TRIANGULAR SYSTEM.
|
|
|
C STORE THE SOLUTION TEMPORARILY IN ZZ().
|
|
|
RTNKEY = 1
|
|
|
GO TO 400
|
|
|
200 CONTINUE
|
|
|
C
|
|
|
C ****** SECONDARY LOOP BEGINS HERE ******
|
|
|
C
|
|
|
C ITERATION COUNTER.
|
|
|
C
|
|
|
210 continue
|
|
|
ITER=ITER+1
|
|
|
IF (ITER .gt. ITMAX) then
|
|
|
MODE=3
|
|
|
write (*,'(/a)') ' NNLS quitting on iteration count.'
|
|
|
GO TO 350
|
|
|
endif
|
|
|
C
|
|
|
C SEE IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE.
|
|
|
C IF NOT COMPUTE ALPHA.
|
|
|
C
|
|
|
ALPHA=TWO
|
|
|
DO 240 IP=1,NSETP
|
|
|
L=INDEX(IP)
|
|
|
IF (ZZ(IP) .le. ZERO) then
|
|
|
T=-X(L)/(ZZ(IP)-X(L))
|
|
|
IF (ALPHA .gt. T) then
|
|
|
ALPHA=T
|
|
|
JJ=IP
|
|
|
endif
|
|
|
endif
|
|
|
240 CONTINUE
|
|
|
C
|
|
|
C IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE THEN ALPHA WILL
|
|
|
C STILL = 2. IF SO EXIT FROM SECONDARY LOOP TO MAIN LOOP.
|
|
|
C
|
|
|
IF (ALPHA.EQ.TWO) GO TO 330
|
|
|
C
|
|
|
C OTHERWISE USE ALPHA WHICH WILL BE BETWEEN 0. AND 1. TO
|
|
|
C INTERPOLATE BETWEEN THE OLD X AND THE NEW ZZ.
|
|
|
C
|
|
|
DO 250 IP=1,NSETP
|
|
|
L=INDEX(IP)
|
|
|
X(L)=X(L)+ALPHA*(ZZ(IP)-X(L))
|
|
|
250 continue
|
|
|
C
|
|
|
C MODIFY A AND B AND THE INDEX ARRAYS TO MOVE COEFFICIENT I
|
|
|
C FROM SET P TO SET Z.
|
|
|
C
|
|
|
I=INDEX(JJ)
|
|
|
260 continue
|
|
|
X(I)=ZERO
|
|
|
C
|
|
|
IF (JJ .ne. NSETP) then
|
|
|
JJ=JJ+1
|
|
|
DO 280 J=JJ,NSETP
|
|
|
II=INDEX(J)
|
|
|
INDEX(J-1)=II
|
|
|
CALL G1 (A(J-1,II),A(J,II),CC,SS,A(J-1,II))
|
|
|
A(J,II)=ZERO
|
|
|
DO 270 L=1,N
|
|
|
IF (L.NE.II) then
|
|
|
c
|
|
|
c Apply procedure G2 (CC,SS,A(J-1,L),A(J,L))
|
|
|
c
|
|
|
TEMP = A(J-1,L)
|
|
|
A(J-1,L) = CC*TEMP + SS*A(J,L)
|
|
|
A(J,L) =-SS*TEMP + CC*A(J,L)
|
|
|
endif
|
|
|
270 CONTINUE
|
|
|
c
|
|
|
c Apply procedure G2 (CC,SS,B(J-1),B(J))
|
|
|
c
|
|
|
TEMP = B(J-1)
|
|
|
B(J-1) = CC*TEMP + SS*B(J)
|
|
|
B(J) =-SS*TEMP + CC*B(J)
|
|
|
280 continue
|
|
|
endif
|
|
|
c
|
|
|
NPP1=NSETP
|
|
|
NSETP=NSETP-1
|
|
|
IZ1=IZ1-1
|
|
|
INDEX(IZ1)=I
|
|
|
C
|
|
|
C SEE IF THE REMAINING COEFFS IN SET P ARE FEASIBLE. THEY SHOULD
|
|
|
C BE BECAUSE OF THE WAY ALPHA WAS DETERMINED.
|
|
|
C IF ANY ARE INFEASIBLE IT IS DUE TO ROUND-OFF ERROR. ANY
|
|
|
C THAT ARE NONPOSITIVE WILL BE SET TO ZERO
|
|
|
C AND MOVED FROM SET P TO SET Z.
|
|
|
C
|
|
|
DO 300 JJ=1,NSETP
|
|
|
I=INDEX(JJ)
|
|
|
IF (X(I) .le. ZERO) go to 260
|
|
|
300 CONTINUE
|
|
|
C
|
|
|
C COPY B( ) INTO ZZ( ). THEN SOLVE AGAIN AND LOOP BACK.
|
|
|
C
|
|
|
DO 310 I=1,M
|
|
|
310 ZZ(I)=B(I)
|
|
|
RTNKEY = 2
|
|
|
GO TO 400
|
|
|
320 CONTINUE
|
|
|
GO TO 210
|
|
|
C ****** END OF SECONDARY LOOP ******
|
|
|
C
|
|
|
330 continue
|
|
|
DO 340 IP=1,NSETP
|
|
|
I=INDEX(IP)
|
|
|
340 X(I)=ZZ(IP)
|
|
|
C ALL NEW COEFFS ARE POSITIVE. LOOP BACK TO BEGINNING.
|
|
|
GO TO 30
|
|
|
C
|
|
|
C ****** END OF MAIN LOOP ******
|
|
|
C
|
|
|
C COME TO HERE FOR TERMINATION.
|
|
|
C COMPUTE THE NORM OF THE FINAL RESIDUAL VECTOR.
|
|
|
C
|
|
|
350 continue
|
|
|
SM=ZERO
|
|
|
IF (NPP1 .le. M) then
|
|
|
DO 360 I=NPP1,M
|
|
|
360 SM=SM+B(I)**2
|
|
|
else
|
|
|
DO 380 J=1,N
|
|
|
380 W(J)=ZERO
|
|
|
endif
|
|
|
RNORM=sqrt(SM)
|
|
|
RETURN
|
|
|
C
|
|
|
C THE FOLLOWING BLOCK OF CODE IS USED AS AN INTERNAL SUBROUTINE
|
|
|
C TO SOLVE THE TRIANGULAR SYSTEM, PUTTING THE SOLUTION IN ZZ().
|
|
|
C
|
|
|
400 continue
|
|
|
DO 430 L=1,NSETP
|
|
|
IP=NSETP+1-L
|
|
|
IF (L .ne. 1) then
|
|
|
DO 410 II=1,IP
|
|
|
ZZ(II)=ZZ(II)-A(II,JJ)*ZZ(IP+1)
|
|
|
410 continue
|
|
|
endif
|
|
|
JJ=INDEX(IP)
|
|
|
ZZ(IP)=ZZ(IP)/A(IP,JJ)
|
|
|
430 continue
|
|
|
go to (200, 320), RTNKEY
|
|
|
END
|
|
|
|
|
|
C SUBROUTINE QRBD (IPASS,Q,E,NN,V,MDV,NRV,C,MDC,NCC)
|
|
|
c
|
|
|
C QR ALGORITHM FOR SINGULAR VALUES OF A BIDIAGONAL MATRIX.
|
|
|
c
|
|
|
c The original version of this code was developed by
|
|
|
c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
|
|
|
c 1973 JUN 12, and published in the book
|
|
|
c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
|
|
|
c Revised FEB 1995 to accompany reprinting of the book by SIAM.
|
|
|
C ------------------------------------------------------------------
|
|
|
C THE BIDIAGONAL MATRIX
|
|
|
C
|
|
|
C (Q1,E2,0... )
|
|
|
C ( Q2,E3,0... )
|
|
|
C D= ( . )
|
|
|
C ( . 0)
|
|
|
C ( .EN)
|
|
|
C ( 0,QN)
|
|
|
C
|
|
|
C IS PRE AND POST MULTIPLIED BY
|
|
|
C ELEMENTARY ROTATION MATRICES
|
|
|
C RI AND PI SO THAT
|
|
|
C
|
|
|
C RK...R1*D*P1**(T)...PK**(T) = DIAG(S1,...,SN)
|
|
|
C
|
|
|
C TO WITHIN WORKING ACCURACY.
|
|
|
C
|
|
|
C 1. EI AND QI OCCUPY E(I) AND Q(I) AS INPUT.
|
|
|
C
|
|
|
C 2. RM...R1*C REPLACES 'C' IN STORAGE AS OUTPUT.
|
|
|
C
|
|
|
C 3. V*P1**(T)...PM**(T) REPLACES 'V' IN STORAGE AS OUTPUT.
|
|
|
C
|
|
|
C 4. SI OCCUPIES Q(I) AS OUTPUT.
|
|
|
C
|
|
|
C 5. THE SI'S ARE NONINCREASING AND NONNEGATIVE.
|
|
|
C
|
|
|
C THIS CODE IS BASED ON THE PAPER AND 'ALGOL' CODE..
|
|
|
C REF..
|
|
|
C 1. REINSCH,C.H. AND GOLUB,G.H. 'SINGULAR VALUE DECOMPOSITION
|
|
|
C AND LEAST SQUARES SOLUTIONS' (NUMER. MATH.), VOL. 14,(1970).
|
|
|
C
|
|
|
C ------------------------------------------------------------------
|
|
|
SUBROUTINE QRBD (IPASS,Q,E,NN,V,MDV,NRV,C,MDC,NCC)
|
|
|
C ------------------------------------------------------------------
|
|
|
integer MDC, MDV, NCC, NN, NRV
|
|
|
c double precision C(MDC,NCC), E(NN), Q(NN),V(MDV,NN)
|
|
|
double precision C(MDC,* ), E(* ), Q(* ),V(MDV,* )
|
|
|
integer I, II, IPASS, J, K, KK, L, LL, LP1, N, N10, NQRS
|
|
|
double precision CS, DIFF, DNORM, F, G, H, SMALL
|
|
|
double precision ONE, SN, T, TEMP, TWO, X, Y, Z, ZERO
|
|
|
|
|
|
logical WNTV ,HAVERS,FAIL
|
|
|
parameter(ONE = 1.0d0, TWO = 2.0d0, ZERO = 0.0d0)
|
|
|
C ------------------------------------------------------------------
|
|
|
N=NN
|
|
|
IPASS=1
|
|
|
IF (N.LE.0) RETURN
|
|
|
N10=10*N
|
|
|
WNTV=NRV.GT.0
|
|
|
HAVERS=NCC.GT.0
|
|
|
FAIL=.FALSE.
|
|
|
NQRS=0
|
|
|
E(1)=ZERO
|
|
|
DNORM=ZERO
|
|
|
DO 10 J=1,N
|
|
|
10 DNORM=max(abs(Q(J))+abs(E(J)),DNORM)
|
|
|
DO 200 KK=1,N
|
|
|
K=N+1-KK
|
|
|
C
|
|
|
C TEST FOR SPLITTING OR RANK DEFICIENCIES..
|
|
|
C FIRST MAKE TEST FOR LAST DIAGONAL TERM, Q(K), BEING SMALL.
|
|
|
20 IF(K.EQ.1) GO TO 50
|
|
|
IF(DIFF(DNORM+Q(K),DNORM) .ne. ZERO) go to 50
|
|
|
C
|
|
|
C SINCE Q(K) IS SMALL WE WILL MAKE A SPECIAL PASS TO
|
|
|
C TRANSFORM E(K) TO ZERO.
|
|
|
C
|
|
|
CS=ZERO
|
|
|
SN=-ONE
|
|
|
DO 40 II=2,K
|
|
|
I=K+1-II
|
|
|
F=-SN*E(I+1)
|
|
|
E(I+1)=CS*E(I+1)
|
|
|
CALL G1 (Q(I),F,CS,SN,Q(I))
|
|
|
C TRANSFORMATION CONSTRUCTED TO ZERO POSITION (I,K).
|
|
|
C
|
|
|
IF (.NOT.WNTV) GO TO 40
|
|
|
DO 30 J=1,NRV
|
|
|
c
|
|
|
c Apply procedure G2 (CS,SN,V(J,I),V(J,K))
|
|
|
c
|
|
|
TEMP = V(J,I)
|
|
|
V(J,I) = CS*TEMP + SN*V(J,K)
|
|
|
V(J,K) =-SN*TEMP + CS*V(J,K)
|
|
|
30 continue
|
|
|
C ACCUMULATE RT. TRANSFORMATIONS IN V.
|
|
|
C
|
|
|
40 CONTINUE
|
|
|
C
|
|
|
C THE MATRIX IS NOW BIDIAGONAL, AND OF LOWER ORDER
|
|
|
C SINCE E(K) .EQ. ZERO..
|
|
|
C
|
|
|
50 DO 60 LL=1,K
|
|
|
L=K+1-LL
|
|
|
IF(DIFF(DNORM+E(L),DNORM) .eq. ZERO) go to 100
|
|
|
IF(DIFF(DNORM+Q(L-1),DNORM) .eq. ZERO) go to 70
|
|
|
60 CONTINUE
|
|
|
C THIS LOOP CAN'T COMPLETE SINCE E(1) = ZERO.
|
|
|
C
|
|
|
GO TO 100
|
|
|
C
|
|
|
C CANCELLATION OF E(L), L.GT.1.
|
|
|
70 CS=ZERO
|
|
|
SN=-ONE
|
|
|
DO 90 I=L,K
|
|
|
F=-SN*E(I)
|
|
|
E(I)=CS*E(I)
|
|
|
IF(DIFF(DNORM+F,DNORM) .eq. ZERO) go to 100
|
|
|
CALL G1 (Q(I),F,CS,SN,Q(I))
|
|
|
IF (HAVERS) then
|
|
|
DO 80 J=1,NCC
|
|
|
c
|
|
|
c Apply procedure G2 ( CS, SN, C(I,J), C(L-1,J)
|
|
|
c
|
|
|
TEMP = C(I,J)
|
|
|
C(I,J) = CS*TEMP + SN*C(L-1,J)
|
|
|
C(L-1,J) =-SN*TEMP + CS*C(L-1,J)
|
|
|
80 continue
|
|
|
endif
|
|
|
90 CONTINUE
|
|
|
C
|
|
|
C TEST FOR CONVERGENCE..
|
|
|
100 Z=Q(K)
|
|
|
IF (L.EQ.K) GO TO 170
|
|
|
C
|
|
|
C SHIFT FROM BOTTOM 2 BY 2 MINOR OF B**(T)*B.
|
|
|
X=Q(L)
|
|
|
Y=Q(K-1)
|
|
|
G=E(K-1)
|
|
|
H=E(K)
|
|
|
F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(TWO*H*Y)
|
|
|
G=sqrt(ONE+F**2)
|
|
|
IF (F .ge. ZERO) then
|
|
|
T=F+G
|
|
|
else
|
|
|
T=F-G
|
|
|
endif
|
|
|
F=((X-Z)*(X+Z)+H*(Y/T-H))/X
|
|
|
C
|
|
|
C NEXT QR SWEEP..
|
|
|
CS=ONE
|
|
|
SN=ONE
|
|
|
LP1=L+1
|
|
|
DO 160 I=LP1,K
|
|
|
G=E(I)
|
|
|
Y=Q(I)
|
|
|
H=SN*G
|
|
|
G=CS*G
|
|
|
CALL G1 (F,H,CS,SN,E(I-1))
|
|
|
F=X*CS+G*SN
|
|
|
G=-X*SN+G*CS
|
|
|
H=Y*SN
|
|
|
Y=Y*CS
|
|
|
IF (WNTV) then
|
|
|
C
|
|
|
C ACCUMULATE ROTATIONS (FROM THE RIGHT) IN 'V'
|
|
|
c
|
|
|
DO 130 J=1,NRV
|
|
|
c
|
|
|
c Apply procedure G2 (CS,SN,V(J,I-1),V(J,I))
|
|
|
c
|
|
|
TEMP = V(J,I-1)
|
|
|
V(J,I-1) = CS*TEMP + SN*V(J,I)
|
|
|
V(J,I) =-SN*TEMP + CS*V(J,I)
|
|
|
130 continue
|
|
|
endif
|
|
|
CALL G1 (F,H,CS,SN,Q(I-1))
|
|
|
F=CS*G+SN*Y
|
|
|
X=-SN*G+CS*Y
|
|
|
IF (HAVERS) then
|
|
|
DO 150 J=1,NCC
|
|
|
c
|
|
|
c Apply procedure G2 (CS,SN,C(I-1,J),C(I,J))
|
|
|
c
|
|
|
TEMP = C(I-1,J)
|
|
|
C(I-1,J) = CS*TEMP + SN*C(I,J)
|
|
|
C(I,J) =-SN*TEMP + CS*C(I,J)
|
|
|
150 continue
|
|
|
endif
|
|
|
c
|
|
|
C APPLY ROTATIONS FROM THE LEFT TO
|
|
|
C RIGHT HAND SIDES IN 'C'..
|
|
|
C
|
|
|
160 CONTINUE
|
|
|
E(L)=ZERO
|
|
|
E(K)=F
|
|
|
Q(K)=X
|
|
|
NQRS=NQRS+1
|
|
|
IF (NQRS.LE.N10) GO TO 20
|
|
|
C RETURN TO 'TEST FOR SPLITTING'.
|
|
|
C
|
|
|
SMALL=ABS(E(K))
|
|
|
I=K
|
|
|
C IF FAILURE TO CONVERGE SET SMALLEST MAGNITUDE
|
|
|
C TERM IN OFF-DIAGONAL TO ZERO. CONTINUE ON.
|
|
|
C ..
|
|
|
DO 165 J=L,K
|
|
|
TEMP=ABS(E(J))
|
|
|
IF(TEMP .EQ. ZERO) GO TO 165
|
|
|
IF(TEMP .LT. SMALL) THEN
|
|
|
SMALL=TEMP
|
|
|
I=J
|
|
|
end if
|
|
|
165 CONTINUE
|
|
|
E(I)=ZERO
|
|
|
NQRS=0
|
|
|
FAIL=.TRUE.
|
|
|
GO TO 20
|
|
|
C ..
|
|
|
C CUTOFF FOR CONVERGENCE FAILURE. 'NQRS' WILL BE 2*N USUALLY.
|
|
|
170 IF (Z.GE.ZERO) GO TO 190
|
|
|
Q(K)=-Z
|
|
|
IF (WNTV) then
|
|
|
DO 180 J=1,NRV
|
|
|
180 V(J,K)=-V(J,K)
|
|
|
endif
|
|
|
190 CONTINUE
|
|
|
C CONVERGENCE. Q(K) IS MADE NONNEGATIVE..
|
|
|
C
|
|
|
200 CONTINUE
|
|
|
IF (N.EQ.1) RETURN
|
|
|
DO 210 I=2,N
|
|
|
IF (Q(I).GT.Q(I-1)) GO TO 220
|
|
|
210 CONTINUE
|
|
|
IF (FAIL) IPASS=2
|
|
|
RETURN
|
|
|
C ..
|
|
|
C EVERY SINGULAR VALUE IS IN ORDER..
|
|
|
220 DO 270 I=2,N
|
|
|
T=Q(I-1)
|
|
|
K=I-1
|
|
|
DO 230 J=I,N
|
|
|
IF (T.GE.Q(J)) GO TO 230
|
|
|
T=Q(J)
|
|
|
K=J
|
|
|
230 CONTINUE
|
|
|
IF (K.EQ.I-1) GO TO 270
|
|
|
Q(K)=Q(I-1)
|
|
|
Q(I-1)=T
|
|
|
IF (HAVERS) then
|
|
|
DO 240 J=1,NCC
|
|
|
T=C(I-1,J)
|
|
|
C(I-1,J)=C(K,J)
|
|
|
240 C(K,J)=T
|
|
|
endif
|
|
|
|
|
|
250 IF (WNTV) then
|
|
|
DO 260 J=1,NRV
|
|
|
T=V(J,I-1)
|
|
|
V(J,I-1)=V(J,K)
|
|
|
260 V(J,K)=T
|
|
|
endif
|
|
|
270 CONTINUE
|
|
|
C END OF ORDERING ALGORITHM.
|
|
|
C
|
|
|
IF (FAIL) IPASS=2
|
|
|
RETURN
|
|
|
END
|
|
|
subroutine SVA(A,MDA,M,N,MDATA,B,SING,KPVEC,NAMES,ISCALE,D,WORK)
|
|
|
c
|
|
|
C SINGULAR VALUE ANALYSIS. COMPUTES THE SINGULAR VALUE
|
|
|
C DECOMPOSITION OF THE MATRIX OF A LEAST SQUARES PROBLEM, AND
|
|
|
C PRODUCES A PRINTED REPORT.
|
|
|
c
|
|
|
c The original version of this code was developed by
|
|
|
c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
|
|
|
c 1973 JUN 12, and published in the book
|
|
|
c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
|
|
|
c Revised FEB 1995 to accompany reprinting of the book by SIAM.
|
|
|
C ------------------------------------------------------------------
|
|
|
c This 1995 version differs from the original 1973 version by the
|
|
|
c addition of the arguments KPVEC() and WORK(), and by allowing user to
|
|
|
c choose the length of names in NAMES().
|
|
|
c KPVEC() allows the user to exercise options regarding printing.
|
|
|
c WORK() provides 2*N locations of work space. Originally SING() was
|
|
|
c required to have 3*N elements, of which the last 2*N were used for
|
|
|
c work space. Now SING() only needs N elements.
|
|
|
c ------------------------------------------------------------------
|
|
|
c Subroutine Arguments
|
|
|
c A(,) [inout] On entry, contains the M x N matrix of the least
|
|
|
c squares problem to be analyzed. This could be a matrix
|
|
|
c obtained by preliminary orthogonal transformations
|
|
|
c applied to the actual problem matrix which may have had
|
|
|
c more rows (See MDATA below.)
|
|
|
c
|
|
|
c MDA [in] First dimensioning parameter for A(,). Require
|
|
|
c MDA .ge. max(M, N).
|
|
|
c
|
|
|
c M,N [in] No. of rows and columns, respectively, in the
|
|
|
c matrix, A. Either M > N or M .le. N is permitted.
|
|
|
c Require M > 0 and N > 0.
|
|
|
c
|
|
|
c MDATA [in] No. of rows in actual least squares problem.
|
|
|
c Generally MDATA .ge. M. MDATA is used only in computing
|
|
|
c statistics for the report and is not used as a loop
|
|
|
c count or array dimension.
|
|
|
c
|
|
|
c B() [inout] On entry, contains the right-side vector, b, of the
|
|
|
c least squares problem. This vector is of length, M.
|
|
|
c On return, contains the vector, g = (U**t)*b, where U
|
|
|
c comes from the singular value decomposition of A. The
|
|
|
c vector , g, is also of length M.
|
|
|
c
|
|
|
c SING() [out] On return, contains the singular values of A, in
|
|
|
c descending order, in locations indexed 1 thru min(M,N).
|
|
|
c If M < N, locations indexed from M+1 through N will be
|
|
|
c set to zero.
|
|
|
c
|
|
|
c KPVEC() [integer array, in] Array of integers to select print
|
|
|
c options. KPVEC(1) determines whether the rest of
|
|
|
c the array is to be used or ignored.
|
|
|
c If KPVEC(1) = 1, the contents of (KPVEC(I), I=2,4)
|
|
|
c will be used to set internal variables as follows:
|
|
|
c PRBLK = KPVEC(2)
|
|
|
c UNIT = KPVEC(3)
|
|
|
c WIDTH = KPVEC(4)
|
|
|
c If KPVEC(1) = 0 default settings will be used. The user
|
|
|
c need not dimension KPVEC() greater than 1. The subr will
|
|
|
c set PRBLK = 111111, UNIT = -1, and WIDTH = 69.
|
|
|
c
|
|
|
c The internal variables PRBLK, UNIT, and WIDTH are
|
|
|
c interpreted as follows:
|
|
|
c
|
|
|
c PRBLK The decimal representation of PRBLK must be
|
|
|
c representable as at most 6 digits, each being 0 or 1.
|
|
|
c The decimal digits will be interpreted as independant
|
|
|
c on/off flags for the 6 possible blocks of printed output.
|
|
|
c Examples: 111111 selects all blocks, 0 suppresses all
|
|
|
c printing, 101010 selects the 1st, 3rd, and 5th blocks,
|
|
|
c etc.
|
|
|
c The six blocks are:
|
|
|
c 1. Header, with size and scaling option parameters.
|
|
|
c 2. V-matrix. Amount of output depends on M and N.
|
|
|
c 3. Singular values and related quantities. Amount of
|
|
|
c output depends on N.
|
|
|
c 4. Listing of YNORM and RNORM and their logarithms.
|
|
|
c Amount of output depends on N.
|
|
|
c 5. Levenberg-Marquart analysis.
|
|
|
c 6. Candidate solutions. Amount of output depends on
|
|
|
c M and N.
|
|
|
c
|
|
|
c UNIT Selects the output unit. If UNIT .ge. 0,
|
|
|
c UNIT will be used as the output unit number.
|
|
|
c If UNIT = -1, output will be written to the "*" output
|
|
|
c unit, i.e., the standard system output unit.
|
|
|
c The calling program unit is responsible for opening
|
|
|
c and/or closing the selected output unit if the host
|
|
|
c system requires these actions.
|
|
|
c
|
|
|
c WIDTH Default value is 79. Determines the width of
|
|
|
c blocks 2, 3, and 6 of the output report.
|
|
|
c Block 3 will use 95(+1) cols if WIDTH .ge. 95, and otherwise
|
|
|
c 69(+1) cols.
|
|
|
c Blocks 2 and 6 are printed by subroutine MFEOUT. These blocks
|
|
|
c generally use at most WIDTH(+1) cols, but will use more if
|
|
|
c the names are so long that more space is needed to print one
|
|
|
c name and one numeric column. The (+1)'s above are reminders
|
|
|
c that in all cases there is one extra initial column for Fortran
|
|
|
c "carriage control". The carriage control character will always
|
|
|
c be a blank.
|
|
|
c Blocks 1, 4, and 5 have fixed widths of 63(+1), 66(+1) and
|
|
|
c 66(+1), respectively.
|
|
|
c
|
|
|
c NAMES() [in] NAMES(j), for j = 1, ..., N, may contain a
|
|
|
c name for the jth component of the solution
|
|
|
c vector. The declared length of the elements of the
|
|
|
c NAMES() array is not specifically limited, but a
|
|
|
c greater length reduces the space available for columns
|
|
|
c of the matris to be printed.
|
|
|
c If NAMES(1) contains only blank characters,
|
|
|
c it will be assumed that no names have been provided,
|
|
|
c and this subr will not access the NAMES() array beyond
|
|
|
c the first element.
|
|
|
C
|
|
|
C ISCALE [in] Set by the user to 1, 2, or 3 to select the column
|
|
|
c scaling option.
|
|
|
C 1 SUBR WILL USE IDENTITY SCALING AND IGNORE THE D()
|
|
|
C ARRAY.
|
|
|
C 2 SUBR WILL SCALE NONZERO COLS TO HAVE UNIT EUCLID-
|
|
|
C EAN LENGTH AND WILL STORE RECIPROCAL LENGTHS OF
|
|
|
C ORIGINAL NONZERO COLS IN D().
|
|
|
C 3 USER SUPPLIES COL SCALE FACTORS IN D(). SUBR
|
|
|
C WILL MULT COL J BY D(J) AND REMOVE THE SCALING
|
|
|
C FROM THE SOLN AT THE END.
|
|
|
c
|
|
|
c D() [ignored or out or in] Usage of D() depends on ISCALE as
|
|
|
c described above. When used, its length must be
|
|
|
c at least N.
|
|
|
c
|
|
|
c WORK() [scratch] Work space of length at least 2*N. Used
|
|
|
c directly in this subr and also in _SVDRS.
|
|
|
c ------------------------------------------------------------------
|
|
|
integer I, IE, IPASS, ISCALE, J, K, KPVEC(4), M, MDA, MDATA
|
|
|
integer MINMN, MINMN1, MPASS, N, NSOL
|
|
|
integer PRBLK, UNIT, WIDTH
|
|
|
double precision A(MDA,N), A1, A2, A3, A4, ALAMB, ALN10, B(M)
|
|
|
double precision D(N), DEL, EL, EL2
|
|
|
double precision ONE, PCOEF, RL, RNORM, RS
|
|
|
double precision SB, SING(N), SL, TEN, THOU, TWENTY
|
|
|
double precision WORK(2*N), YL, YNORM, YS, YSQ, ZERO
|
|
|
character*(*) NAMES(N)
|
|
|
logical BLK(6), NARROW, STAR
|
|
|
parameter( ZERO = 0.0d0, ONE = 1.0d0)
|
|
|
parameter( TEN = 10.0d0, TWENTY = 20.0d0, THOU = 1000.0d0)
|
|
|
c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
220 format (1X/' INDEX SING. VAL. P COEF ',
|
|
|
* ' RECIPROCAL G COEF G**2 ',
|
|
|
* ' CUMULATIVE SCALED SQRT'/
|
|
|
* 31x,' SING. VAL.',26x,
|
|
|
* ' SUM of SQRS of CUM.S.S.')
|
|
|
221 format (1X/' INDEX SING. VAL. P COEF ',
|
|
|
* ' RECIPROCAL G COEF SCALED SQRT'/
|
|
|
* 31x,' SING. VAL.',13x,' of CUM.S.S.')
|
|
|
222 format (1X/' INDEX SING. VAL. G COEF G**2 ',
|
|
|
* ' CUMULATIVE SCALED SQRT'/
|
|
|
* 44x,' SUM of SQRS of CUM.S.S.')
|
|
|
|
|
|
230 format (' ',4X,'0',64X,2g13.4)
|
|
|
231 format (' ',4X,'0',51X, g13.4)
|
|
|
232 format (' ',4X,'0',38X,2g13.4)
|
|
|
|
|
|
240 format (' ',i5,g12.4,6g13.4)
|
|
|
|
|
|
260 format (1X,' M = ',I6,', N =',I4,', MDATA =',I8)
|
|
|
270 format (1X/' Singular Value Analysis of the least squares',
|
|
|
* ' problem, A*X = B,'/
|
|
|
* ' scaled as (A*D)*Y = B.')
|
|
|
280 format (1X/' Scaling option No.',I2,'. D is a diagonal',
|
|
|
* ' matrix with the following diagonal elements..'/(5X,10E12.4))
|
|
|
290 format (1X/' Scaling option No. 1. D is the identity matrix.'/
|
|
|
* 1X)
|
|
|
300 format (1X/' INDEX',12X,'YNORM RNORM',11X,
|
|
|
* ' LOG10 LOG10'/
|
|
|
* 45X,' YNORM RNORM'/1X)
|
|
|
310 format (' ',I5,6X,2E11.3,11X,2F11.3)
|
|
|
320 format (1X/
|
|
|
*' Norms of solution and residual vectors for a range of values'/
|
|
|
*' of the Levenberg-Marquardt parameter, LAMBDA.'//
|
|
|
* ' LAMBDA YNORM RNORM',
|
|
|
* ' LOG10 LOG10 LOG10'/
|
|
|
* 34X,' LAMBDA YNORM RNORM')
|
|
|
330 format (1X, 3E11.3, 3F11.3)
|
|
|
c ------------------------------------------------------------------
|
|
|
IF (M.LE.0 .OR. N.LE.0) RETURN
|
|
|
MINMN = min(M,N)
|
|
|
MINMN1 = MINMN + 1
|
|
|
if(KPVEC(1) .eq. 0) then
|
|
|
PRBLK = 111111
|
|
|
UNIT = -1
|
|
|
WIDTH = 79
|
|
|
else
|
|
|
PRBLK = KPVEC(2)
|
|
|
UNIT = KPVEC(3)
|
|
|
WIDTH = KPVEC(4)
|
|
|
endif
|
|
|
STAR = UNIT .lt. 0
|
|
|
c Build logical array BLK() by testing
|
|
|
c decimal digits of PRBLK.
|
|
|
do 20 I=6, 1, -1
|
|
|
J = PRBLK/10
|
|
|
BLK(I) = (PRBLK - 10*J) .gt. 0
|
|
|
PRBLK = J
|
|
|
20 continue
|
|
|
c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
c Optionally print header and M, N, MDATA
|
|
|
if(BLK(1)) then
|
|
|
if(STAR) then
|
|
|
write (*,270)
|
|
|
write (*,260) M,N,MDATA
|
|
|
else
|
|
|
write (UNIT,270)
|
|
|
write (UNIT,260) M,N,MDATA
|
|
|
endif
|
|
|
endif
|
|
|
c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
c Handle scaling as selected by ISCALE.
|
|
|
if( ISCALE .eq. 1) then
|
|
|
if(BLK(1)) then
|
|
|
if(STAR) then
|
|
|
write (*,290)
|
|
|
else
|
|
|
write (UNIT,290)
|
|
|
endif
|
|
|
endif
|
|
|
else
|
|
|
C
|
|
|
C Apply column scaling to A.
|
|
|
C
|
|
|
DO 52 J = 1,N
|
|
|
A1 = D(J)
|
|
|
if( ISCALE .le. 2) then
|
|
|
SB = ZERO
|
|
|
DO 30 I = 1,M
|
|
|
30 SB = SB + A(I,J)**2
|
|
|
A1 = sqrt(SB)
|
|
|
IF (A1.EQ.ZERO) A1 = ONE
|
|
|
A1 = ONE/A1
|
|
|
D(J) = A1
|
|
|
endif
|
|
|
DO 50 I = 1,M
|
|
|
A(I,J) = A(I,J)*A1
|
|
|
50 continue
|
|
|
52 continue
|
|
|
if(BLK(1)) then
|
|
|
if(STAR) then
|
|
|
write (*,280) ISCALE,(D(J),J = 1,N)
|
|
|
else
|
|
|
write (UNIT,280) ISCALE,(D(J),J = 1,N)
|
|
|
endif
|
|
|
endif
|
|
|
endif
|
|
|
c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
C Compute the Singular Value Decomposition of the scaled matrix.
|
|
|
C
|
|
|
call SVDRS (A,MDA,M,N,B,M,1,SING,WORK)
|
|
|
c
|
|
|
c Determine NSOL.
|
|
|
NSOL = MINMN
|
|
|
do 60 J = 1,MINMN
|
|
|
if(SING(J) .eq. ZERO) then
|
|
|
NSOL = J-1
|
|
|
go to 65
|
|
|
endif
|
|
|
60 continue
|
|
|
65 continue
|
|
|
C
|
|
|
c The array B() contains the vector G.
|
|
|
C Compute cumulative sums of squares of components of
|
|
|
C G and store them in WORK(I), I = 1,...,MINMN+1
|
|
|
C
|
|
|
SB = ZERO
|
|
|
DO 70 I = MINMN1,M
|
|
|
SB = SB + B(I)**2
|
|
|
70 CONTINUE
|
|
|
WORK(MINMN+1) = SB
|
|
|
DO 75 J = MINMN, 1, -1
|
|
|
SB = SB + B(J)**2
|
|
|
WORK(J) = SB
|
|
|
75 CONTINUE
|
|
|
c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
C PRINT THE V MATRIX.
|
|
|
C
|
|
|
if(BLK(2)) CALL MFEOUT (A,MDA,N,N,NAMES,1, UNIT, WIDTH)
|
|
|
c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
C REPLACE V BY D*V IN THE ARRAY A()
|
|
|
if (ISCALE .gt.1) then
|
|
|
do 82 I = 1,N
|
|
|
do 80 J = 1,N
|
|
|
A(I,J) = D(I)*A(I,J)
|
|
|
80 continue
|
|
|
82 continue
|
|
|
endif
|
|
|
c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
if(BLK(3)) then
|
|
|
c Print singular values and other summary results.
|
|
|
c
|
|
|
c Output will be done using one of two layouts. The narrow
|
|
|
c layout uses 69 cols + 1 for carriage control, and makes two passes
|
|
|
c through the computation.
|
|
|
c The wide layout uses 95 cols + 1 for carriage control, and makes
|
|
|
c only one pass through the computation.
|
|
|
c
|
|
|
C G NOW IN B() ARRAY. V NOW IN A(,) ARRAY.
|
|
|
C
|
|
|
NARROW = WIDTH .lt. 95
|
|
|
MPASS = 1
|
|
|
if(NARROW) MPASS = 2
|
|
|
do 170 IPASS = 1, MPASS
|
|
|
if(STAR) then
|
|
|
if(NARROW) then
|
|
|
if(IPASS .eq. 1) then
|
|
|
write(*,221)
|
|
|
else
|
|
|
write(*,222)
|
|
|
endif
|
|
|
else
|
|
|
write (*,220)
|
|
|
endif
|
|
|
else
|
|
|
if(NARROW) then
|
|
|
if(IPASS .eq. 1) then
|
|
|
write(UNIT,221)
|
|
|
else
|
|
|
write(UNIT,222)
|
|
|
endif
|
|
|
else
|
|
|
write (UNIT,220)
|
|
|
endif
|
|
|
endif
|
|
|
c The following stmt converts from
|
|
|
c integer to floating-point.
|
|
|
A3 = WORK(1)
|
|
|
A4 = sqrt(A3/ max(1,MDATA))
|
|
|
if(STAR) then
|
|
|
if(NARROW) then
|
|
|
if(IPASS .eq. 1) then
|
|
|
write(*,231) A4
|
|
|
else
|
|
|
write(*,232) A3, A4
|
|
|
endif
|
|
|
else
|
|
|
write (*,230) A3,A4
|
|
|
endif
|
|
|
else
|
|
|
if(NARROW) then
|
|
|
if(IPASS .eq. 1) then
|
|
|
write(UNIT,231) A4
|
|
|
else
|
|
|
write(UNIT,232) A3, A4
|
|
|
endif
|
|
|
else
|
|
|
write (UNIT,230) A3,A4
|
|
|
endif
|
|
|
endif
|
|
|
C
|
|
|
DO 160 K = 1,MINMN
|
|
|
if (SING(K).EQ.ZERO) then
|
|
|
PCOEF = ZERO
|
|
|
if(STAR) then
|
|
|
write (*,240) K,SING(K)
|
|
|
else
|
|
|
write (UNIT,240) K,SING(K)
|
|
|
endif
|
|
|
else
|
|
|
PCOEF = B(K) / SING(K)
|
|
|
A1 = ONE / SING(K)
|
|
|
A2 = B(K)**2
|
|
|
A3 = WORK(K+1)
|
|
|
A4 = sqrt(A3/max(1,MDATA-K))
|
|
|
if(STAR) then
|
|
|
if(NARROW) then
|
|
|
if(IPASS .eq. 1) then
|
|
|
write(*,240) K,SING(K),PCOEF,A1,B(K), A4
|
|
|
else
|
|
|
write(*,240) K,SING(K), B(K),A2,A3,A4
|
|
|
endif
|
|
|
else
|
|
|
write (*,240) K,SING(K),PCOEF,A1,B(K),A2,A3,A4
|
|
|
endif
|
|
|
else
|
|
|
if(NARROW) then
|
|
|
if(IPASS .eq. 1) then
|
|
|
write(UNIT,240) K,SING(K),PCOEF,A1,B(K), A4
|
|
|
else
|
|
|
write(UNIT,240) K,SING(K), B(K),A2,A3,A4
|
|
|
endif
|
|
|
else
|
|
|
write (UNIT,240) K,SING(K),PCOEF,A1,B(K),A2,A3,A4
|
|
|
endif
|
|
|
|
|
|
endif
|
|
|
endif
|
|
|
160 continue
|
|
|
170 continue
|
|
|
endif
|
|
|
c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
if( BLK(4) ) then
|
|
|
C
|
|
|
C Compute and print values of YNORM, RNORM and their logarithms.
|
|
|
C
|
|
|
if(STAR) then
|
|
|
write (*,300)
|
|
|
else
|
|
|
write (UNIT,300)
|
|
|
endif
|
|
|
YSQ = ZERO
|
|
|
do 180 J = 0, NSOL
|
|
|
if(J .ne. 0) YSQ = YSQ + (B(J) / SING(J))**2
|
|
|
YNORM = sqrt(YSQ)
|
|
|
RNORM = sqrt(WORK(J+1))
|
|
|
YL = -THOU
|
|
|
IF (YNORM .GT. ZERO) YL = log10(YNORM)
|
|
|
RL = -THOU
|
|
|
IF (RNORM .GT. ZERO) RL = log10(RNORM)
|
|
|
if(STAR) then
|
|
|
write (*,310) J,YNORM,RNORM,YL,RL
|
|
|
else
|
|
|
write (UNIT,310) J,YNORM,RNORM,YL,RL
|
|
|
endif
|
|
|
180 continue
|
|
|
endif
|
|
|
c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
if( BLK(5) .and. SING(1) .ne. ZERO ) then
|
|
|
C
|
|
|
C COMPUTE VALUES OF XNORM AND RNORM FOR A SEQUENCE OF VALUES OF
|
|
|
C THE LEVENBERG-MARQUARDT PARAMETER.
|
|
|
C
|
|
|
EL = log10(SING(1)) + ONE
|
|
|
EL2 = log10(SING(NSOL)) - ONE
|
|
|
DEL = (EL2-EL) / TWENTY
|
|
|
ALN10 = log(TEN)
|
|
|
if(STAR) then
|
|
|
write (*,320)
|
|
|
else
|
|
|
write (UNIT,320)
|
|
|
endif
|
|
|
DO 200 IE = 1,21
|
|
|
C COMPUTE ALAMB = 10.0**EL
|
|
|
ALAMB = EXP(ALN10*EL)
|
|
|
YS = ZERO
|
|
|
RS = WORK(NSOL+1)
|
|
|
DO 190 I = 1,MINMN
|
|
|
SL = SING(I)**2 + ALAMB**2
|
|
|
YS = YS + (B(I)*SING(I)/SL)**2
|
|
|
RS = RS + (B(I)*(ALAMB**2)/SL)**2
|
|
|
190 CONTINUE
|
|
|
YNORM = sqrt(YS)
|
|
|
RNORM = sqrt(RS)
|
|
|
RL = -THOU
|
|
|
IF (RNORM.GT.ZERO) RL = log10(RNORM)
|
|
|
YL = -THOU
|
|
|
IF (YNORM.GT.ZERO) YL = log10(YNORM)
|
|
|
if(STAR) then
|
|
|
write (*,330) ALAMB,YNORM,RNORM,EL,YL,RL
|
|
|
else
|
|
|
write (UNIT,330) ALAMB,YNORM,RNORM,EL,YL,RL
|
|
|
endif
|
|
|
EL = EL + DEL
|
|
|
200 CONTINUE
|
|
|
endif
|
|
|
c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
C Compute and optionally print candidate solutions.
|
|
|
C
|
|
|
do 215 K = 1,NSOL
|
|
|
PCOEF = B(K) / SING(K)
|
|
|
DO 210 I = 1,N
|
|
|
A(I,K) = A(I,K) * PCOEF
|
|
|
210 IF (K.GT.1) A(I,K) = A(I,K) + A(I,K-1)
|
|
|
215 continue
|
|
|
if (BLK(6) .and. NSOL.GE.1)
|
|
|
* CALL MFEOUT (A,MDA,N,NSOL,NAMES,2,UNIT,WIDTH)
|
|
|
return
|
|
|
END
|
|
|
C SUBROUTINE SVDRS (A, MDA, M1, N1, B, MDB, NB, S, WORK)
|
|
|
c
|
|
|
C SINGULAR VALUE DECOMPOSITION ALSO TREATING RIGHT SIDE VECTOR.
|
|
|
c
|
|
|
c The original version of this code was developed by
|
|
|
c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
|
|
|
c 1974 SEP 25, and published in the book
|
|
|
c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
|
|
|
c Revised FEB 1995 to accompany reprinting of the book by SIAM.
|
|
|
C ------------------------------------------------------------------
|
|
|
c This 1995 version differs from the original 1974 version by adding
|
|
|
c the argument WORK().
|
|
|
c WORK() provides 2*N1 locations of work space. Originally S() was
|
|
|
c required to have 3*N1 elements, of which the last 2*N1 were used for
|
|
|
c work space. Now S() only needs N1 elements.
|
|
|
C ------------------------------------------------------------------
|
|
|
c This subroutine computes the singular value decomposition of the
|
|
|
c given M1 x N1 matrix, A, and optionally applys the transformations
|
|
|
c from the left to the NB column vectors of the M1 x NB matrix B.
|
|
|
c Either M1 .ge. N1 or M1 .lt. N1 is permitted.
|
|
|
c
|
|
|
c The singular value decomposition of A is of the form
|
|
|
c
|
|
|
c A = U * S * V**t
|
|
|
c
|
|
|
c where U is M1 x M1 orthogonal, S is M1 x N1 diagonal with the
|
|
|
c diagonal terms nonnegative and ordered from large to small, and
|
|
|
c V is N1 x N1 orthogonal. Note that these matrices also satisfy
|
|
|
c
|
|
|
c S = (U**t) * A * V
|
|
|
c
|
|
|
c The matrix V is returned in the leading N1 rows and
|
|
|
c columns of the array A(,).
|
|
|
c
|
|
|
c The singular values, i.e. the diagonal terms of the matrix S,
|
|
|
c are returned in the array S(). If M1 .lt. N1, positions M1+1
|
|
|
c through N1 of S() will be set to zero.
|
|
|
c
|
|
|
c The product matrix G = U**t * B replaces the given matrix B
|
|
|
c in the array B(,).
|
|
|
c
|
|
|
c If the user wishes to obtain a minimum length least squares
|
|
|
c solution of the linear system
|
|
|
c
|
|
|
c A * X ~=~ B
|
|
|
c
|
|
|
c the solution X can be constructed, following use of this subroutine,
|
|
|
c by computing the sum for i = 1, ..., R of the outer products
|
|
|
c
|
|
|
c (Col i of V) * (1/S(i)) * (Row i of G)
|
|
|
c
|
|
|
c Here R denotes the pseudorank of A which the user may choose
|
|
|
c in the range 0 through Min(M1, N1) based on the sizes of the
|
|
|
c singular values.
|
|
|
C ------------------------------------------------------------------
|
|
|
C Subroutine Arguments
|
|
|
c
|
|
|
c A(,) (In/Out) On input contains the M1 x N1 matrix A.
|
|
|
c On output contains the N1 x N1 matrix V.
|
|
|
c
|
|
|
c LDA (In) First dimensioning parameter for A(,).
|
|
|
c Require LDA .ge. Max(M1, N1).
|
|
|
c
|
|
|
c M1 (In) No. of rows of matrices A, B, and G.
|
|
|
c Require M1 > 0.
|
|
|
c
|
|
|
c N1 (In) No. of cols of matrix A, No. of rows and cols of
|
|
|
c matrix V. Permit M1 .ge. N1 or M1 .lt. N1.
|
|
|
c Require N1 > 0.
|
|
|
c
|
|
|
c B(,) (In/Out) If NB .gt. 0 this array must contain an
|
|
|
c M1 x NB matrix on input and will contain the
|
|
|
c M1 x NB product matrix, G = (U**t) * B on output.
|
|
|
c
|
|
|
c LDB (In) First dimensioning parameter for B(,).
|
|
|
c Require LDB .ge. M1.
|
|
|
c
|
|
|
c NB (In) No. of cols in the matrices B and G.
|
|
|
c Require NB .ge. 0.
|
|
|
c
|
|
|
c S() (Out) Must be dimensioned at least N1. On return will
|
|
|
c contain the singular values of A, with the ordering
|
|
|
c S(1) .ge. S(2) .ge. ... .ge. S(N1) .ge. 0.
|
|
|
c If M1 .lt. N1 the singular values indexed from M1+1
|
|
|
c through N1 will be zero.
|
|
|
c If the given integer arguments are not consistent, this
|
|
|
c subroutine will return immediately, setting S(1) = -1.0.
|
|
|
c
|
|
|
c WORK() (Scratch) Work space of total size at least 2*N1.
|
|
|
c Locations 1 thru N1 will hold the off-diagonal terms of
|
|
|
c the bidiagonal matrix for subroutine QRBD. Locations N1+1
|
|
|
c thru 2*N1 will save info from one call to the next of
|
|
|
c H12.
|
|
|
c ------------------------------------------------------------------
|
|
|
C This code gives special treatment to rows and columns that are
|
|
|
c entirely zero. This causes certain zero sing. vals. to appear as
|
|
|
c exact zeros rather than as about MACHEPS times the largest sing. val.
|
|
|
c It similarly cleans up the associated columns of U and V.
|
|
|
c
|
|
|
c METHOD..
|
|
|
c 1. EXCHANGE COLS OF A TO PACK NONZERO COLS TO THE LEFT.
|
|
|
c SET N = NO. OF NONZERO COLS.
|
|
|
c USE LOCATIONS A(1,N1),A(1,N1-1),...,A(1,N+1) TO RECORD THE
|
|
|
c COL PERMUTATIONS.
|
|
|
c 2. EXCHANGE ROWS OF A TO PACK NONZERO ROWS TO THE TOP.
|
|
|
c QUIT PACKING IF FIND N NONZERO ROWS. MAKE SAME ROW EXCHANGES
|
|
|
c IN B. SET M SO THAT ALL NONZERO ROWS OF THE PERMUTED A
|
|
|
c ARE IN FIRST M ROWS. IF M .LE. N THEN ALL M ROWS ARE
|
|
|
c NONZERO. IF M .GT. N THEN THE FIRST N ROWS ARE KNOWN
|
|
|
c TO BE NONZERO,AND ROWS N+1 THRU M MAY BE ZERO OR NONZERO.
|
|
|
c 3. APPLY ORIGINAL ALGORITHM TO THE M BY N PROBLEM.
|
|
|
c 4. MOVE PERMUTATION RECORD FROM A(,) TO S(I),I=N+1,...,N1.
|
|
|
c 5. BUILD V UP FROM N BY N TO N1 BY N1 BY PLACING ONES ON
|
|
|
c THE DIAGONAL AND ZEROS ELSEWHERE. THIS IS ONLY PARTLY DONE
|
|
|
c EXPLICITLY. IT IS COMPLETED DURING STEP 6.
|
|
|
c 6. EXCHANGE ROWS OF V TO COMPENSATE FOR COL EXCHANGES OF STEP 2.
|
|
|
c 7. PLACE ZEROS IN S(I),I=N+1,N1 TO REPRESENT ZERO SING VALS.
|
|
|
c ------------------------------------------------------------------
|
|
|
subroutine SVDRS (A, MDA, M1, N1, B, MDB, NB, S, WORK)
|
|
|
integer I, IPASS, J, K, L, M, MDA, MDB, M1
|
|
|
integer N, NB, N1, NP1, NS, NSP1
|
|
|
c double precision A(MDA,N1),B(MDB,NB), S(N1)
|
|
|
double precision A(MDA, *),B(MDB, *), S( *)
|
|
|
double precision ONE, T, WORK(N1,2), ZERO
|
|
|
parameter(ONE = 1.0d0, ZERO = 0.0d0)
|
|
|
c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
C BEGIN.. SPECIAL FOR ZERO ROWS AND COLS.
|
|
|
C
|
|
|
C PACK THE NONZERO COLS TO THE LEFT
|
|
|
C
|
|
|
N=N1
|
|
|
IF (N.LE.0.OR.M1.LE.0) RETURN
|
|
|
J=N
|
|
|
10 CONTINUE
|
|
|
DO 20 I=1,M1
|
|
|
IF (A(I,J) .ne. ZERO) go to 50
|
|
|
20 CONTINUE
|
|
|
C
|
|
|
C COL J IS ZERO. EXCHANGE IT WITH COL N.
|
|
|
C
|
|
|
IF (J .ne. N) then
|
|
|
DO 30 I=1,M1
|
|
|
30 A(I,J)=A(I,N)
|
|
|
endif
|
|
|
A(1,N)=J
|
|
|
N=N-1
|
|
|
50 CONTINUE
|
|
|
J=J-1
|
|
|
IF (J.GE.1) GO TO 10
|
|
|
C IF N=0 THEN A IS ENTIRELY ZERO AND SVD
|
|
|
C COMPUTATION CAN BE SKIPPED
|
|
|
NS=0
|
|
|
IF (N.EQ.0) GO TO 240
|
|
|
C PACK NONZERO ROWS TO THE TOP
|
|
|
C QUIT PACKING IF FIND N NONZERO ROWS
|
|
|
I=1
|
|
|
M=M1
|
|
|
60 IF (I.GT.N.OR.I.GE.M) GO TO 150
|
|
|
IF (A(I,I)) 90,70,90
|
|
|
70 DO 80 J=1,N
|
|
|
IF (A(I,J)) 90,80,90
|
|
|
80 CONTINUE
|
|
|
GO TO 100
|
|
|
90 I=I+1
|
|
|
GO TO 60
|
|
|
C ROW I IS ZERO
|
|
|
C EXCHANGE ROWS I AND M
|
|
|
100 IF(NB.LE.0) GO TO 115
|
|
|
DO 110 J=1,NB
|
|
|
T=B(I,J)
|
|
|
B(I,J)=B(M,J)
|
|
|
110 B(M,J)=T
|
|
|
115 DO 120 J=1,N
|
|
|
120 A(I,J)=A(M,J)
|
|
|
IF (M.GT.N) GO TO 140
|
|
|
DO 130 J=1,N
|
|
|
130 A(M,J)=ZERO
|
|
|
140 CONTINUE
|
|
|
C EXCHANGE IS FINISHED
|
|
|
M=M-1
|
|
|
GO TO 60
|
|
|
C
|
|
|
150 CONTINUE
|
|
|
C END.. SPECIAL FOR ZERO ROWS AND COLUMNS
|
|
|
C BEGIN.. SVD ALGORITHM
|
|
|
C METHOD..
|
|
|
C (1) REDUCE THE MATRIX TO UPPER BIDIAGONAL FORM WITH
|
|
|
C HOUSEHOLDER TRANSFORMATIONS.
|
|
|
C H(N)...H(1)AQ(1)...Q(N-2) = (D**T,0)**T
|
|
|
C WHERE D IS UPPER BIDIAGONAL.
|
|
|
C
|
|
|
C (2) APPLY H(N)...H(1) TO B. HERE H(N)...H(1)*B REPLACES B
|
|
|
C IN STORAGE.
|
|
|
C
|
|
|
C (3) THE MATRIX PRODUCT W= Q(1)...Q(N-2) OVERWRITES THE FIRST
|
|
|
C N ROWS OF A IN STORAGE.
|
|
|
C
|
|
|
C (4) AN SVD FOR D IS COMPUTED. HERE K ROTATIONS RI AND PI ARE
|
|
|
C COMPUTED SO THAT
|
|
|
C RK...R1*D*P1**(T)...PK**(T) = DIAG(S1,...,SM)
|
|
|
C TO WORKING ACCURACY. THE SI ARE NONNEGATIVE AND NONINCREASING.
|
|
|
C HERE RK...R1*B OVERWRITES B IN STORAGE WHILE
|
|
|
C A*P1**(T)...PK**(T) OVERWRITES A IN STORAGE.
|
|
|
C
|
|
|
C (5) IT FOLLOWS THAT,WITH THE PROPER DEFINITIONS,
|
|
|
C U**(T)*B OVERWRITES B, WHILE V OVERWRITES THE FIRST N ROW AND
|
|
|
C COLUMNS OF A.
|
|
|
C
|
|
|
L=min(M,N)
|
|
|
C THE FOLLOWING LOOP REDUCES A TO UPPER BIDIAGONAL AND
|
|
|
C ALSO APPLIES THE PREMULTIPLYING TRANSFORMATIONS TO B.
|
|
|
C
|
|
|
DO 170 J=1,L
|
|
|
IF (J.GE.M) GO TO 160
|
|
|
CALL H12 (1,J,J+1,M,A(1,J),1,T,A(1,J+1),1,MDA,N-J)
|
|
|
CALL H12 (2,J,J+1,M,A(1,J),1,T,B,1,MDB,NB)
|
|
|
160 IF (J.GE.N-1) GO TO 170
|
|
|
CALL H12 (1,J+1,J+2,N,A(J,1),MDA,work(J,2),A(J+1,1),MDA,1,M-J)
|
|
|
170 CONTINUE
|
|
|
C
|
|
|
C COPY THE BIDIAGONAL MATRIX INTO S() and WORK() FOR QRBD.
|
|
|
C 1986 Jan 8. C. L. Lawson. Changed N to L in following 2 statements.
|
|
|
IF (L.EQ.1) GO TO 190
|
|
|
DO 180 J=2,L
|
|
|
S(J)=A(J,J)
|
|
|
180 WORK(J,1)=A(J-1,J)
|
|
|
190 S(1)=A(1,1)
|
|
|
C
|
|
|
NS=N
|
|
|
IF (M.GE.N) GO TO 200
|
|
|
NS=M+1
|
|
|
S(NS)=ZERO
|
|
|
WORK(NS,1)=A(M,M+1)
|
|
|
200 CONTINUE
|
|
|
C
|
|
|
C CONSTRUCT THE EXPLICIT N BY N PRODUCT MATRIX, W=Q1*Q2*...*QL*I
|
|
|
C IN THE ARRAY A().
|
|
|
C
|
|
|
DO 230 K=1,N
|
|
|
I=N+1-K
|
|
|
IF (I .GT. min(M,N-2)) GO TO 210
|
|
|
CALL H12 (2,I+1,I+2,N,A(I,1),MDA,WORK(I,2),A(1,I+1),1,MDA,N-I)
|
|
|
210 DO 220 J=1,N
|
|
|
220 A(I,J)=ZERO
|
|
|
230 A(I,I)=ONE
|
|
|
C
|
|
|
C COMPUTE THE SVD OF THE BIDIAGONAL MATRIX
|
|
|
C
|
|
|
CALL QRBD (IPASS,S(1),WORK(1,1),NS,A,MDA,N,B,MDB,NB)
|
|
|
C
|
|
|
if(IPASS .eq. 2) then
|
|
|
write (*,'(/a)')
|
|
|
* ' FULL ACCURACY NOT ATTAINED IN BIDIAGONAL SVD'
|
|
|
endif
|
|
|
|
|
|
240 CONTINUE
|
|
|
IF (NS.GE.N) GO TO 260
|
|
|
NSP1=NS+1
|
|
|
DO 250 J=NSP1,N
|
|
|
250 S(J)=ZERO
|
|
|
260 CONTINUE
|
|
|
IF (N.EQ.N1) RETURN
|
|
|
NP1=N+1
|
|
|
C MOVE RECORD OF PERMUTATIONS
|
|
|
C AND STORE ZEROS
|
|
|
DO 280 J=NP1,N1
|
|
|
S(J)=A(1,J)
|
|
|
DO 270 I=1,N
|
|
|
270 A(I,J)=ZERO
|
|
|
280 CONTINUE
|
|
|
C PERMUTE ROWS AND SET ZERO SINGULAR VALUES.
|
|
|
DO 300 K=NP1,N1
|
|
|
I=S(K)
|
|
|
S(K)=ZERO
|
|
|
DO 290 J=1,N1
|
|
|
A(K,J)=A(I,J)
|
|
|
290 A(I,J)=ZERO
|
|
|
A(I,K)=ONE
|
|
|
300 CONTINUE
|
|
|
C END.. SPECIAL FOR ZERO ROWS AND COLUMNS
|
|
|
RETURN
|
|
|
END
|
|
|
|
|
|
|
|
|
|
