This package of Fortran 77 and Fortran 90 code accompanies the SIAM Publications printing of "Solving Least Squares Problems," by C. Lawson and R. Hanson. The routines in the package are filed as shown in the list below. The units listed 1.-14. may be compiled and placed in a library. Routine G2 is no longer used but is included for completeness. The main program units 15.-18., 20.-21. require this library for all external references. The item 19. is a data file read by program unit PROG4.FOR. The main program 23., PROG7.F90, requires the subprogram 22. or BVLS.F90. 1. BNDACC.FOR 2. BNDSOL.FOR 3. DIFF.FOR 4. G1.FOR 5. G2.FOR 6. GEN.FOR 7. H12.FOR 8. HFTI.FOR 9. LDP.FOR 10. MFEOUT.FOR 11. NNLS.FOR 12. QRBD.FOR 13. SVA.FOR 14. SVDRS.FOR 15. PROG1.FOR 16. PROG2.FOR 17. PROG3.FOR 18. PROG4.FOR 19. DATA4.DAT 20. PROG5.FOR 21. PROG6.FOR 22. BVLS.F90 23. PROG7.F90 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 GEN(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 GEN= 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 GEN=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 program PROG1 c C DEMONSTRATE ALGORITHMS HFT AND HS1 FOR SOLVING LEAST SQUARES C PROBLEMS AND ALGORITHM COV FOR COMPUTING THE ASSOCIATED COVARIANCE C MATRICES. 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 MDA parameter(MDA = 8) integer I, IP1, J, JM1, K, L integer M, MMN, MN1, MN2, N, NM1, NOISE, NP1, NPJ double precision A(MDA,MDA), ANOISE, B(MDA), DUMMY, GEN, H(MDA) double precision ONE, SM, SMALL, SRSMSQ, ZERO parameter(ONE = 1.0d0, SMALL = 1.0d-4, ZERO = 0.0d0) c ------------------------------------------------------------------ 300 format (/1x,8X,'ESTIMATED PARAMETERS, X=A**(+)*B,', * ' COMPUTED BY ''HFT,HS1'''// * (9X,I6,E16.8,I6,E16.8,I6,E16.8,I6,E16.8,I6,E16.8)) 305 format (/1x,8X,'RESIDUAL LENGTH =',E12.4) 310 format (/1x,8X,'COVARIANCE MATRIX (UNSCALED) OF', * ' ESTIMATED PARAMETERS'/ * 9x,'COMPUTED BY ''COV''.'/1X) 320 format (9X,2I3,E16.8,2I3,E16.8,2I3,E16.8,2I3,E16.8,2I3,E16.8) 330 format (/' PROG1. THIS PROGRAM DEMONSTRATES THE ALGORITHMS', * ' HFT, HS1, AND COV.'// * ' CAUTION.. ''PROG1'' DOES NO CHECKING FOR', * ' NEAR RANK DEFICIENT MATRICES.'/ * ' RESULTS IN SUCH CASES MAY BE MEANINGLESS.'/ * ' SUCH CASES ARE HANDLED BY ''PROG2'' OR ''PROG3''') 340 format (/ * ' THE RELATIVE NOISE LEVEL OF THE GENERATED DATA WILL BE ',g11.3) 350 format (/////' M N'/1X,2I4) 360 format (/2x, * '****** TERMINATING THIS CASE DUE TO', * ' A DIVISOR BEING EXACTLY ZERO. ******') c ------------------------------------------------------------------ DO 230 NOISE=1,2 if(NOISE .eq. 1) then ANOISE= ZERO else ANOISE= SMALL endif write (*,330) write (*,340) ANOISE C INITIALIZE THE DATA GENERATION FUNCTION C .. DUMMY=GEN(-ONE) DO 220 MN1=1,6,5 MN2=MN1+2 DO 210 M=MN1,MN2 DO 200 N=MN1, M NP1=N+1 write (*,350) M,N C GENERATE DATA C .. DO 10 I=1,M DO 10 J=1,N 10 A(I,J)=GEN(ANOISE) DO 20 I=1,M 20 B(I)=GEN(ANOISE) IF(M .LT. N) GO TO 180 C C ****** BEGIN ALGORITHM HFT ****** C .. DO 30 J=1,N CALL H12 (1,J,J+1,M,A(1,J),1,H(J), * A(1,min(J+1,N)),1,MDA,N-J) 30 continue C .. C THE ALGORITHM 'HFT' IS COMPLETED. C C ****** BEGIN ALGORITHM HS1 ****** C APPLY THE TRANSFORMATIONS Q(N)...Q(1)=Q TO B C REPLACING THE PREVIOUS CONTENTS OF THE ARRAY, B . C .. DO 40 J=1,N 40 CALL H12 (2,J,J+1,M,A(1,J),1,H(J),B,1,1,1) C SOLVE THE TRIANGULAR SYSTEM FOR THE SOLUTION X. C STORE X IN THE ARRAY, B . C .. DO 80 K=1,N I=NP1-K SM= ZERO IF (I.EQ.N) GO TO 60 IP1=I+1 DO 50 J=IP1,N 50 SM=SM+A(I,J)*B(J) 60 IF (A(I,I) .eq. ZERO) then write (*,360) GO TO 180 endif B(I)=(B(I)-SM)/A(I,I) 80 continue C COMPUTE LENGTH OF RESIDUAL VECTOR. C .. SRSMSQ= ZERO IF (N.EQ.M) GO TO 100 MMN=M-N DO 90 J=1,MMN NPJ=N+J 90 SRSMSQ=SRSMSQ+B(NPJ)**2 SRSMSQ=SQRT(SRSMSQ) C ****** BEGIN ALGORITHM COV ****** C COMPUTE UNSCALED COVARIANCE MATRIX ((A**T)*A)**(-1) C .. 100 DO 110 J=1,N 110 A(J,J)= ONE/A(J,J) IF (N.EQ.1) GO TO 140 NM1=N-1 DO 130 I=1,NM1 IP1=I+1 DO 130 J=IP1,N JM1=J-1 SM= ZERO DO 120 L=I,JM1 120 SM=SM+A(I,L)*DBLE(A(L,J)) 130 A(I,J)=-SM*A(J,J) C .. C THE UPPER TRIANGLE OF A HAS BEEN INVERTED C UPON ITSELF. 140 DO 160 I=1,N DO 160 J=I,N SM= ZERO DO 150 L=J,N 150 SM=SM+A(I,L)*DBLE(A(J,L)) 160 A(I,J)=SM C .. C THE UPPER TRIANGULAR PART OF THE C SYMMETRIC MATRIX (A**T*A)**(-1) HAS C REPLACED THE UPPER TRIANGULAR PART OF C THE A ARRAY. write (*,300) (I,B(I),I=1,N) write (*,305) SRSMSQ write (*,310) DO 170 I=1,N 170 write (*,320) (I,J,A(I,J),J=I,N) 180 CONTINUE 200 continue 210 continue 220 continue 230 continue STOP END program PROG2 c C DEMONSTRATE ALGORITHM HFTI FOR SOLVING LEAST SQUARES PROBLEMS C AND ALGORITHM COV FOR COMPUTING THE ASSOCIATED UNSCALED C COVARIANCE 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 ------------------------------------------------------------------ integer MDA parameter(MDA = 8) integer I, II, IP(MDA), IP1, J, JM1, K, KM1, KP1, KRANK, L integer M, MN1, MN2, N, NM1, NOISE double precision A(MDA,MDA), ANOISE, ANORM, B(MDA) double precision DUMMY, FIVE00, G(MDA), GEN double precision H(MDA), HALF double precision ONE, SM, SMALL, SRSMSQ, TAU, TMP, ZERO parameter(FIVE00 = 500.0d0, HALF = 0.5d0) parameter(ONE = 1.0d0, SMALL = 1.0d-4, ZERO = 0.0d0) C ------------------------------------------------------------------ 190 format (/1x,8X,'RESIDUAL LENGTH = ',E12.4) 200 format (/1x,8X, * 'ESTIMATED PARAMETERS, X=A**(+)*B, COMPUTED BY ''HFTI'''// * (9X,I6,E16.8,I6,E16.8,I6,E16.8,I6,E16.8,I6,E16.8)) 210 format (/1x,8X, * 'COVARIANCE MATRIX (UNSCALED) OF ESTIMATED PARAMETERS'/ * 9x,'COMPUTED BY ''COV''.'/1X) 220 format (9X,2I3,E16.8,2I3,E16.8,2I3,E16.8,2I3,E16.8,2I3,E16.8) 230 format (/' PROG2. THIS PROGRAM DEMONSTATES THE ALGORITHMS'/ * 9x,'HFTI AND COV.') 240 format (/ * ' THE RELATIVE NOISE LEVEL OF THE GENERATED DATA WILL BE ', * g11.3/ * ' THE MATRIX NORM IS APPROXIMATELY ',E12.4/ * ' THE ABSOLUTE PSEUDORANK TOLERANCE, TAU, IS ',E12.4) 250 format (/////' M N'/1X,2I4) 260 format (/1x,8X,'PSEUDORANK = ',I4) C ------------------------------------------------------------------ DO 180 NOISE=1,2 ANORM= FIVE00 ANOISE= ZERO TAU= HALF IF (NOISE.EQ.1) GO TO 10 ANOISE= SMALL TAU=ANORM*ANOISE*10. 10 CONTINUE C INITIALIZE THE DATA GENERATION FUNCTION C .. DUMMY=GEN(-ONE) write (*,230) write (*,240) ANOISE,ANORM,TAU C DO 180 MN1=1,6,5 MN2=MN1+2 DO 180 M=MN1,MN2 DO 180 N=MN1,MN2 write (*,250) M,N C GENERATE DATA C .. DO 20 I=1,M DO 20 J=1,N 20 A(I,J)=GEN(ANOISE) DO 30 I=1,M 30 B(I)=GEN(ANOISE) C C ****** CALL HFTI ****** C CALL HFTI(A,MDA,M,N,B,1,1,TAU,KRANK,SRSMSQ,H,G,IP) C write (*,260) KRANK write (*,200) (I,B(I),I=1,N) write (*,190) SRSMSQ IF (KRANK.LT.N) GO TO 180 C ****** ALGORITHM COV BEGINS HERE ****** C .. DO 40 J=1,N 40 A(J,J)= ONE/A(J,J) IF (N.EQ.1) GO TO 70 NM1=N-1 DO 60 I=1,NM1 IP1=I+1 DO 60 J=IP1,N JM1=J-1 SM= ZERO DO 50 L=I,JM1 50 SM=SM+A(I,L)*A(L,J) 60 A(I,J)=-SM*A(J,J) C .. C THE UPPER TRIANGLE OF A HAS BEEN INVERTED C UPON ITSELF. 70 DO 90 I=1,N DO 90 J=I,N SM= ZERO DO 80 L=J,N 80 SM=SM+A(I,L)*DBLE(A(J,L)) 90 A(I,J)=SM IF (N.LT.2) GO TO 160 DO 150 II=2,N I=N+1-II IF (IP(I).EQ.I) GO TO 150 K=IP(I) TMP=A(I,I) A(I,I)=A(K,K) A(K,K)=TMP IF (I.EQ.1) GO TO 110 DO 100 L=2,I TMP=A(L-1,I) A(L-1,I)=A(L-1,K) 100 A(L-1,K)=TMP 110 IP1=I+1 KM1=K-1 IF (IP1.GT.KM1) GO TO 130 DO 120 L=IP1,KM1 TMP=A(I,L) A(I,L)=A(L,K) 120 A(L,K)=TMP 130 IF (K.EQ.N) GO TO 150 KP1=K+1 DO 140 L=KP1,N TMP=A(I,L) A(I,L)=A(K,L) 140 A(K,L)=TMP 150 CONTINUE 160 CONTINUE C .. C COVARIANCE HAS BEEN COMPUTED AND REPERMUTED. C THE UPPER TRIANGULAR PART OF THE C SYMMETRIC MATRIX (A**T*A)**(-1) HAS C REPLACED THE UPPER TRIANGULAR PART OF C THE A ARRAY. write (*,210) DO 170 I=1,N 170 write (*,220) (I,J,A(I,J),J=I,N) 180 CONTINUE STOP END program PROG3 c C DEMONSTRATE THE USE OF SUBROUTINE SVDRS TO COMPUTE THE C SINGULAR VALUE DECOMPOSITION OF A MATRIX, A, AND SOLVE A LEAST C SQUARES PROBLEM, A*X=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 C THE S.V.D. A= U*(S,0)**T*V**T IS C COMPUTED SO THAT.. C (1) U**T*B REPLACES THE M+1 ST. COL. OF B. C C (2) U**T REPLACES THE M BY M IDENTITY IN C THE FIRST M COLS. OF B. C C (3) V REPLACES THE FIRST N ROWS AND COLS. C OF A. C C (4) THE DIAGONAL ENTRIES OF THE S MATRIX C REPLACE THE FIRST N ENTRIES OF THE ARRAY S. C C THE ARRAY S( ) MUST BE DIMENSIONED AT LEAST 3*N . C ------------------------------------------------------------------ integer MDA, MDB, MX parameter(MDA = 8, MDB = 8, MX = 8) integer I, J, KP1, KRANK, L, M, MINMN, MN1, MN2, N, NOISE double precision A(MDA, MX), AA(MDA, MX), ANOISE, B(MDB, MDB+1) double precision DN, DUMMY, FLOATN, GEN, ONE, RHO double precision S(MX), SM, SMALL3, SMALL4, SRSMSQ double precision T, TAU, TEN, WORK(2*MX), X(MX), ZERO parameter(ONE = 1.0d0, SMALL3 = 1.0d-3, SMALL4 = 1.0d-4) parameter(TEN = 10.0d0) parameter(ZERO = 0.0d0) c ------------------------------------------------------------------ 170 format (/////' M N'/1X,2I4) 180 format (/' SINGULAR VALUES OF MATRIX') 190 format (/' TRANSFORMED RIGHT SIDE, U**T*B') 200 format (/ * ' ESTIMATED PARAMETERS, X=A**(+)*B, COMPUTED BY ''SVDRS''') 210 format (/' RESIDUAL VECTOR LENGTH = ',E12.4) 220 format (/ * ' FROBENIUS NORM(A-U*(S,0)**T*V**T)'/ * ' --------------------------------- = ',g12.4/ * ' SQRT(N) * SPECTRAL NORM OF A') 230 format (2X,I5,E16.8,I5,E16.8,I5,E16.8) 240 format (/ * ' PROG3. THIS PROGRAM DEMONSTRATES THE ALGORITHM SVDRS.'/) 250 format (/ * ' THE RELATIVE NOISE LEVEL OF THE GENERATED DATA WILL BE', * g12.4/ * ' THE RELATIVE TOLERANCE FOR PSEUDORANK DETERMINATION IS RHO =', * g12.4) 260 format (/' ABSOLUTE PSEUDORANK TOLERANCE, TAU = ',g12.4/ * ' PSEUDORANK = ',I5) c ------------------------------------------------------------------ DO 166 NOISE=1,2 if(NOISE .eq. 1) then ANOISE= ZERO RHO= SMALL3 else ANOISE= SMALL4 RHO= TEN * ANOISE endif write (*,240) write (*,250) ANOISE,RHO C INITIALIZE DATA GENERATION FUNCTION C .. DUMMY=GEN(-ONE) C DO 164 MN1=1,6,5 MN2=MN1+2 DO 162 M=MN1,MN2 DO 160 N=MN1,MN2 write (*,170) M,N DO 35 I=1,M DO 20 J=1,M 20 B(I,J)= ZERO B(I,I)= ONE DO 30 J=1,N A(I,J)=GEN(ANOISE) AA(I,J)=A(I,J) 30 continue 35 continue DO 40 I=1,M 40 B(I,M+1)=GEN(ANOISE) C C The arrays are now filled in.. C Compute the singular value decomposition. C ****************************************************** CALL SVDRS (A,MDA,M,N,B,MDB,M+1,S, WORK) C ****************************************************** C write (*,180) write (*,230) (I,S(I),I=1,N) write (*,190) write (*,230) (I,B(I,M+1),I=1,M) C C TEST FOR DISPARITY OF RATIO OF SINGULAR VALUES. C .. KRANK=N TAU=RHO*S(1) DO 50 I=1,N IF (S(I).LE.TAU) GO TO 60 50 CONTINUE GO TO 70 60 KRANK=I-1 70 continue write (*,260) TAU,KRANK C COMPUTE SOLUTION VECTOR ASSUMING PSEUDORANK IS KRANK C .. DO 80 I=1,KRANK 80 B(I,M+1)=B(I,M+1)/S(I) DO 100 I=1,N SM= ZERO DO 90 J=1,KRANK 90 SM=SM+A(I,J)*B(J,M+1) 100 X(I)=SM C COMPUTE PREDICTED NORM OF RESIDUAL VECTOR. C .. SRSMSQ= ZERO IF (KRANK .ne. M) then KP1=KRANK+1 DO 110 I=KP1,M 110 SRSMSQ=SRSMSQ+B(I,M+1)**2 SRSMSQ=sqrt(SRSMSQ) endif write (*,200) write (*,230) (I,X(I),I=1,N) write (*,210) SRSMSQ C COMPUTE THE FROBENIUS NORM OF A**T- V*(S,0)*U**T. C C COMPUTE V*S FIRST. C MINMN=min(M,N) DO 130 J=1,MINMN DO 130 I=1,N A(I,J)=A(I,J)*S(J) 130 continue 135 continue DN= ZERO DO 150 J=1,M DO 145 I=1,N SM= ZERO DO 140 L=1,MINMN 140 SM=SM+A(I,L)*B(L,J) C COMPUTED DIFFERENCE OF (I,J) TH ENTRY C OF A**T-V*(S,0)*U**T. C .. T=AA(J,I)-SM DN=DN+T**2 145 continue 150 continue FLOATN = N DN=sqrt(DN)/(sqrt(FLOATN)*S(1)) write (*,220) DN 160 continue 162 continue 164 continue 166 continue stop end program PROG4 c C DEMONSTRATE 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 ------------------------------------------------------------------ integer MDA, MX parameter(MDA = 15, MX = 5) integer I, J, KPVEC(4), M, N double precision A(MDA,MX), B(MDA), D(MX), SING(MX), WORK(2*MX) character NAMES(MX)*8 data KPVEC / 1, 111111, -1, 76 / data NAMES/' Fire',' Water',' Earth',' Air','Cosmos'/ C ------------------------------------------------------------------ M = MDA N = MX write (*,'(a)') * ' PROG4. DEMONSTRATE SINGULAR VALUE ANALYSIS', * ' PROG4 will read data from file: data4.dat' open(unit= 10, file= 'data4.dat', status= 'old') read (10,'(6F12.0)') ((A(I,J),J=1,N),B(I),I=1,M) write (*,'(/a)') * ' LISTING OF INPUT MATRIX, A, AND VECTOR, B, FOLLOWS..' write (*,'(/(1x,f11.7,4F12.7,F13.4))') ((A(I,J),J=1,N),B(I),I=1,M) write (*,'(/)') C call SVA (A,MDA,M,N,MDA,B,SING,KPVEC, NAMES,1,D, WORK) C end -.13405547 -.20162827 -.16930778 -.18971990 -.17387234 -.4361 -.10379475 -.15766336 -.13346256 -.14848550 -.13597690 -.3437 -.08779597 -.12883867 -.10683007 -.12011796 -.10932972 -.2657 .02058554 .00335331 -.01641270 .00078606 .00271659 -.0392 -.03248093 -.01876799 .00410639 -.01405894 -.01384391 .0193 .05967662 .06667714 .04352153 .05740438 .05024962 .0747 .06712457 .07352437 .04489770 .06471862 .05876455 .0935 .08687186 .09368296 .05672327 .08141043 .07302320 .1079 .02149662 .06222662 .07213486 .06200069 .05570931 .1930 .06687407 .10344506 .09153849 .09508223 .08393667 .2058 .15879069 .18088339 .11540692 .16160727 .14796479 .2606 .17642887 .20361830 .13057860 .18385729 .17005549 .3142 .11414080 .17259611 .14816471 .16007466 .14374096 .3529 .07846038 .14669563 .14365800 .14003842 .12571177 .3615 .10803175 .16994623 .14971519 .15885312 .14301547 .3647 program PROG5 c C EXAMPLE OF THE USE OF SUBROUTINES BNDACC AND BNDSOL TO SOLVE C SEQUENTIALLY THE BANDED LEAST SQUARES PROBLEM THAT ARISES IN C SPLINE CURVE FITTING. 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 MDG, MXY parameter(MDG = 12, MXY = 12) integer I, IC, IG, IP, IR, J, JT, L, M, MT, NBAND, NBP, NC double precision B(10), C(MXY), COV(MXY), G(MDG,5), H double precision ONE, P1, P2, Q(4), R, RDUMMY, RNORM double precision SIGFAC, SIGSQ, U, X(MXY), Y(MXY), YFIT, ZERO parameter(ONE = 1.0d0, ZERO = 0.0d0) data Y /2.2d0, 4.0d0, 5.0d0, 4.6d0, 2.8d0, 2.7d0, * 3.8d0, 5.1d0, 6.1d0, 6.3d0, 5.0d0, 2.0d0/ C ------------------------------------------------------------------ 160 format (/' I',8X,'X',10X,'Y',6X,'YFIT',4X,'R=Y-YFIT/1X') 170 format (1X,I3,4X,F6.0,4X,F6.2,4X,F6.2,4X,F8.4) 180 format (/' C =',6F10.5/(4X,6F10.5)) 190 format (3(2X,2I4,g15.7)) 200 format (/' COVARIANCE MATRIX OF THE SPLINE COEFFICIENTS.') 210 format (/' RNORM =',g15.8) 220 format (/' SIGFAC =',g15.8) 230 format (/' PROG5. EXECUTE A SEQUENCE OF CUBIC SPLINE FITS'/ * 9x,'TO A DISCRETE SET OF DATA.') 240 format (////' NEW CASE..'/ * ' NUMBER OF BREAKPOINTS, INCLUDING ENDPOINTS, IS ',I5/ * ' BREAKPOINTS =',6f10.5,/(14X,6f10.5)) C ------------------------------------------------------------------ write (*,230) NBAND=4 M=MDG C SET ABCISSAS OF DATA DO 10 I=1,M 10 X(I)=2*I C C BEGIN LOOP THRU CASES USING INCREASING NOS OF BREAKPOINTS. C DO 150 NBP=5,10 NC=NBP+2 C SET BREAKPOINTS B(1)=X(1) B(NBP)=X(M) H=(B(NBP)-B(1))/(NBP-1) IF (NBP.LE.2) GO TO 30 DO 20 I=3,NBP 20 B(I-1)=B(I-2)+H 30 CONTINUE write (*,240) NBP,(B(I),I=1,NBP) C C INITIALIZE IR AND IP BEFORE FIRST CALL TO BNDACC. C IR=1 IP=1 I=1 JT=1 40 MT=0 50 CONTINUE IF (X(I).GT.B(JT+1)) GO TO 60 C C SET ROW FOR ITH DATA POINT C U=(X(I)-B(JT))/H IG=IR+MT G(IG,1)=P1(ONE - U) G(IG,2)=P2(ONE - U) G(IG,3)=P2(U) G(IG,4)=P1(U) G(IG,5)=Y(I) MT=MT+1 IF (I.EQ.M) GO TO 60 I=I+1 GO TO 50 C C SEND BLOCK OF DATA TO PROCESSOR C 60 CONTINUE CALL BNDACC (G,MDG,NBAND,IP,IR,MT,JT) IF (I.EQ.M) GO TO 70 JT=JT+1 GO TO 40 C COMPUTE SOLUTION C() 70 CONTINUE CALL BNDSOL (1,G,MDG,NBAND,IP,IR,C,NC,RNORM) C WRITE SOLUTION COEFFICIENTS write (*,180) (C(L),L=1,NC) write (*,210) RNORM C C COMPUTE AND PRINT X,Y,YFIT,R=Y-YFIT C write (*,160) JT=1 DO 110 I=1,M 80 IF (X(I).LE.B(JT+1)) GO TO 90 JT=JT+1 GO TO 80 C 90 U=(X(I)-B(JT))/H Q(1)=P1(ONE - U) Q(2)=P2(ONE - U) Q(3)=P2(U) Q(4)=P1(U) YFIT=ZERO DO 100 L=1,4 IC=JT-1+L 100 YFIT=YFIT+C(IC)*Q(L) R=Y(I)-YFIT write (*,170) I,X(I),Y(I),YFIT,R 110 CONTINUE C C COMPUTE RESIDUAL VECTOR NORM. C IF (M.LE.NC) GO TO 150 SIGSQ=(RNORM**2)/(M-NC) SIGFAC=sqrt(SIGSQ) write (*,220) SIGFAC write (*,200) C C COMPUTE AND PRINT COLS. OF COVARIANCE. C DO 140 J=1,NC DO 120 I=1,NC 120 COV(I)=ZERO COV(J)= ONE CALL BNDSOL (2,G,MDG,NBAND,IP,IR,COV,NC,RDUMMY) CALL BNDSOL (3,G,MDG,NBAND,IP,IR,COV,NC,RDUMMY) C C COMPUTE THE JTH COL. OF THE COVARIANCE MATRIX. DO 130 I=1,NC 130 COV(I)=COV(I)*SIGSQ 140 write (*,190) (L,J,COV(L),L=1,NC) 150 CONTINUE STOP END c ================================================================== C C DEFINE FUNCTIONS P1 AND P2 TO BE USED IN CONSTRUCTING C CUBIC SPLINES OVER UNIFORMLY SPACED BREAKPOINTS. C c ================================================================== double precision function P1(T) double precision T P1 = 0.25d0 * T**2 * T return end c ================================================================== double precision function P2(T) double precision T P2 = -(1.0d0 - T)**2 * (1.0d0 + T) * 0.75d0 + 1.0d0 return end program PROG6 c C DEMONSTRATE THE USE OF THE SUBROUTINE LDP FOR LEAST C DISTANCE PROGRAMMING BY SOLVING THE CONSTRAINED LINE DATA FITTING C PROBLEM OF CHAPTER 23. 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 ------------------------------------------------------------------ integer MX parameter(MX = 2) integer I, INDEX(3), J, L, MDE, MDGH, ME, MG, MODE, N, NP1 double precision E(4,MX), F(4), G(3,MX), G2(3,MX) double precision H(3), H2(3), ONE, RES, S(MX), SM, T(4) double precision W(4), WLDP(21), WORK(2*MX), X(MX) double precision Z(MX), ZERO, ZNORM parameter(ONE = 1.0d0, ZERO = 0.0d0) data T / 0.25d0, 0.5d0, 0.5d0, 0.8d0 / data W / 0.50d0, 0.6d0, 0.7d0, 1.2d0 / C ------------------------------------------------------------------ 110 format (' PROG6.. EXAMPLE OF CONSTRAINED CURVE FITTING'/ * 10x,'USING THE SUBROUTINE LDP.'// * 10x,'RELATED INTERMEDIATE QUANTITIES ARE PRINTED.') 120 format (/' V = ',2F10.5/(10X,2F10.5)) 130 format (/' F TILDA =',4F10.5/' S = ',2F10.5) 140 format (/' G TILDA =',2F10.5/(10X,2F10.5)) 150 format (/' H TILDA =',3F10.5) 160 format (/' Z = ',2F10.5) 170 format (/' THE COEFICIENTS OF THE FITTED LINE F(T)=X(1)*T+X(2)'/ * ' ARE X(1) = ',F10.5,' AND X(2) = ',F10.5) 180 format (/' THE CONSECUTIVE RESIDUALS ARE'/1X,4(I4,F10.5)) 190 format (/' THE RESIDUALS NORM IS ',F10.5) 200 format (/' MODE (FROM LDP) = ',I3,', ZNORM = ',F10.5) C ------------------------------------------------------------------ write (*,110) MDE=4 MDGH=3 C ME=4 MG=3 N=2 C DEFINE THE LEAST SQUARES AND CONSTRAINT MATRICES. C DO 10 I=1,ME E(I,1)=T(I) E(I,2)= ONE 10 F(I)=W(I) C G(1,1)= ONE G(1,2)= ZERO G(2,1)= ZERO G(2,2)= ONE G(3,1)=-ONE G(3,2)=-ONE C H(1)= ZERO H(2)= ZERO H(3)=-ONE C C COMPUTE THE SINGULAR VALUE DECOMPOSITION OF THE MATRIX, E. C CALL SVDRS (E, MDE, ME, N, F, 1, 1, S, WORK) C write (*,120) ((E(I,J),J=1,N),I=1,N) write (*,130) F,(S(J),J=1,N) C C GENERALLY RANK DETERMINATION AND LEVENBERG-MARQUARDT C STABILIZATION COULD BE INSERTED HERE. C C DEFINE THE CONSTRAINT MATRIX FOR THE Z COORDINATE SYSTEM. DO 30 I=1,MG DO 30 J=1,N SM= ZERO DO 20 L=1,N 20 SM=SM+G(I,L)*E(L,J) 30 G2(I,J)=SM/S(J) C DEFINE CONSTRAINT RT SIDE FOR THE Z COORDINATE SYSTEM. DO 50 I=1,MG SM= ZERO DO 40 J=1,N 40 SM=SM+G2(I,J)*F(J) 50 H2(I)=H(I)-SM C write (*,140) ((G2(I,J),J=1,N),I=1,MG) write (*,150) H2 C C SOLVE THE CONSTRAINED PROBLEM IN Z-COORDINATES. C CALL LDP (G2,MDGH,MG,N,H2,Z,ZNORM,WLDP,INDEX,MODE) C write (*,200) MODE,ZNORM write (*,160) Z C C TRANSFORM BACK FROM Z-COORDINATES TO X-COORDINATES. DO 60 J=1,N 60 Z(J)=(Z(J)+F(J))/S(J) DO 80 I=1,N SM= ZERO DO 70 J=1,N 70 SM=SM+E(I,J)*Z(J) 80 X(I)=SM RES=ZNORM**2 NP1=N+1 DO 90 I=NP1,ME 90 RES=RES+F(I)**2 RES=sqrt(RES) C COMPUTE THE RESIDUALS. DO 100 I=1,ME 100 F(I)=W(I)-X(1)*T(I)-X(2) write (*,170) (X(J),J=1,N) write (*,180) (I,F(I),I=1,ME) write (*,190) RES STOP END SUBROUTINE BVLS (A, B, BND, X, RNORM, NSETP, W, INDEX, IERR) ! Given an M by N matrix, A, and an M-vector, B, compute an ! N-vector, X, that solves the least-squares problem A *X = B ! subject to X(J) satisfying BND(1,J) <= X(J) <= BND(2,J) ! ! The values BND(1,J) = -huge(ONE) and BND(2,J) = huge(ONE) are ! suggested choices to designate that there is no constraint in that ! direction. The parameter ONE is 1.0 in the working precision. ! ! This algorithm is a generalization of NNLS, that solves ! the least-squares problem, A * X = B, subject to all X(J) >= 0. ! The subroutine NNLS appeared in 'SOLVING LEAST SQUARES PROBLEMS,' ! by Lawson and Hanson, Prentice-Hall, 1974. Work on BVLS was started ! by C. L. Lawson and R. J. Hanson at Jet Propulsion Laboratory, ! 1973 June 12. Many modifications were subsequently made. ! This Fortran 90 code was completed in April, 1995 by R. J. Hanson. ! Assumed shape arrays, automatic arrays, array operations, internal ! subroutines and F90 elemental functions are used. Comments are ! prefixed with ! More than 72 characters appear on some lines. ! The BVLS package is an additional item for the reprinting of the book ! by SIAM Publications and the distribution of the code package ! using netlib and Internet or network facilities. !INTERFACE ! SUBROUTINE BVLS (A, B, BND, X, RNORM, NSETP, W, INDEX, IERR) ! REAL(KIND(1E0)) A(:,:), B(:), BND(:,:), X(:), RNORM, W(:) ! INTEGER NSETP, INDEX(:), IERR ! END SUBROUTINE !END INTERFACE ! A(:,:) [INTENT(InOut)] ! On entry A() contains the M by N matrix, A. ! On return A() contains the product matrix, Q*A, where ! Q is an M by M orthogonal matrix generated by this ! subroutine. The dimensions are M=size(A,1) and N=size(A,2). ! ! B(:) [INTENT(InOut)] ! On entry B() contains the M-vector, B. ! On return, B() contains Q*B. The same Q multiplies A. ! ! BND(:,:) [INTENT(In)] ! BND(1,J) is the lower bound for X(J). ! BND(2,J) is the upper bound for X(J). ! Require: BND(1,J) <= BND(2,J). ! ! X(:) [INTENT(Out)] ! On entry X() need not be initialized. On return, ! X() will contain the solution N-vector. ! ! RNORM [INTENT(Out)] ! The Euclidean norm of the residual vector, b - A*X. ! ! NSETP [INTENT(Out)] ! Indicates the number of components of the solution ! vector, X(), that are not at their constraint values. ! ! W(:) [INTENT(Out)] ! An N-array. On return, W() will contain the dual solution ! vector. Using Set definitions below: ! W(J) = 0 for all j in Set P, ! W(J) <= 0 for all j in Set Z, such that X(J) is at its ! lower bound, and ! W(J) >= 0 for all j in Set Z, such that X(J) is at its ! upper bound. ! If BND(1,J) = BND(2,J), so the variable X(J) is fixed, ! then W(J) will have an arbitrary value. ! ! INDEX(:) [INTENT(Out)] ! An INTEGER working array of size N. On exit the contents ! of this array define the sets P, Z, and F as follows: ! ! INDEX(1) through INDEX(NSETP) = Set P. ! INDEX(IZ1) through INDEX(IZ2) = Set Z. ! INDEX(IZ2+1) through INDEX(N) = Set F. ! IZ1 = NSETP + 1 = NPP1 ! Any of these sets may be empty. Set F is those components ! that are constrained to a unique value by the given ! constraints. Sets P and Z are those that are allowed a non- ! zero range of values. Of these, set Z are those whose final ! value is a constraint value, while set P are those whose ! final value is not a constraint. The value of IZ2 is not returned. !... It is computable as the number of bounds constraining a component !... of X uniquely. ! ! IERR [INTENT(Out)] ! Indicates status on return. ! = 0 Solution completed. ! = 1 M <= 0 or N <= 0 ! = 2 B(:), X(:), BND(:,:), W(:), or INDEX(:) size or shape violation. ! = 3 Input bounds are inconsistent. ! = 4 Exceed maximum number of iterations. ! ! Selected internal variables: ! EPS [real(kind(one))] ! Determines the relative linear dependence of a column vector ! for a variable moved from its initial value. This is used in ! one place with the default value EPS=EPSILON(ONE). Other ! values, larger or smaller may be needed for some problems. ! Library software will likely make this an optional argument. ! ! ITMAX [integer] ! Set to 3*N. Maximum number of iterations permitted. ! This is usually larger than required. Library software will ! likely make this an optional argument. ! ITER [integer] ! Iteration counter. implicit none logical FIND, HITBND, FREE1, FREE2, FREE integer IERR, M, N integer I, IBOUND, II, IP, ITER, ITMAX, IZ, IZ1, IZ2 integer J, JJ, JZ, L, LBOUND, NPP1, NSETP integer INDEX(:) ! ! Z() [Scratch] An M-array (automatic) of working space. ! S() [Scratch] An N-array (automatic) of working space. ! real(kind(1E0)), parameter :: ZERO = 0E0, ONE=1E0, TWO = 2E0 real(kind(one)) :: A(:,:), B(:), S(size(A,2)), X(:), W(:),& Z(size(A,1)), BND(:,:) real(kind(one)) ALPHA, ASAVE, CC, EPS, RANGE, RNORM real(kind(one)) NORM, SM, SS, T, UNORM, UP, ZTEST CALL INITIALIZE ! ! The above call will set IERR. LOOPA: DO ! ! Quit on error flag, or if all coefficients are already in the ! solution, .or. if M columns of A have been triangularized. IF (IERR /= 0 .or. IZ1 > IZ2 .or. NSETP >= M) exit LOOPA ! CALL SELECT_ANOTHER_COEFF_TO!_SOLVE_FOR ! ! See if no index was found to be moved from set Z to set P. ! Then go to termination. IF ( .not. FIND ) exit LOOPA ! CALL MOVE_J_FROM_SET_Z_TO_SET_P ! CALL TEST_SET_P_AGAINST_CONSTRAINTS ! ! The above call may set IERR. ! All coefficients in set P are strictly feasible. Loop back. END DO LOOPA ! CALL TERMINATION RETURN CONTAINS ! These are internal subroutines. SUBROUTINE INITIALIZE M=size(A,1); N=size(A,2) IF (M <= 0 .or. N <= 0) then IERR = 1 RETURN END IF ! ! Check array sizes for consistency and with M and N. IF(SIZE(X) < N) THEN IERR=2 RETURN END IF IF(SIZE(B) < M) THEN IERR=2 RETURN END IF IF(SIZE(BND,1) /= 2) THEN IERR=2 RETURN END IF IF(SIZE(BND,2) < N) THEN IERR=2 RETURN END IF IF(SIZE(W) < N) THEN IERR=2 RETURN END IF IF(SIZE(INDEX) < N) THEN IERR=2 RETURN END IF IERR = 0 ! The next two parameters can be changed to be optional for library software. EPS = EPSILON(ONE) ITMAX=3*N ITER=0 ! Initialize the array index(). DO I=1,N INDEX(I)=I END DO ! IZ2=N IZ1=1 NSETP=0 NPP1=1 ! ! Begin: Loop on IZ to initialize X(). IZ=IZ1 DO IF (IZ > IZ2 ) EXIT J=INDEX(IZ) IF ( BND(1,J) <= -huge(ONE)) then IF (BND(2,J) >= huge(ONE)) then X(J) = ZERO else X(J) = min(ZERO,BND(2,J)) END IF ELSE IF ( BND(2,J) >= huge(ONE)) then X(J) = max(ZERO,BND(1,J)) else RANGE = BND(2,J) - BND(1,J) IF ( RANGE <= ZERO ) then ! ! Here X(J) is constrained to a single value. INDEX(IZ)=INDEX(IZ2) INDEX(IZ2)=J IZ=IZ-1 IZ2=IZ2-1 X(J)=BND(1,J) W(J)=ZERO ELSE IF ( RANGE > ZERO) then ! ! The following statement sets X(J) to 0 if the constraint interval ! includes 0, and otherwise sets X(J) to the endpoint of the ! constraint interval that is closest to 0. ! X(J) = max(BND(1,J), min(BND(2,J),ZERO)) else IERR = 3 RETURN END IF! ( RANGE:.) END IF ! IF ( abs(X(J)) > ZERO ) then ! ! Change B() to reflect a nonzero starting value for X(J). B(1:M)=B(1:M)-A(1:M,J)*X(J) END IF IZ=IZ+1 END DO! ( IZ <= IZ2 ) END SUBROUTINE! ( INITIALIZE ) SUBROUTINE SELECT_ANOTHER_COEFF_TO!_SOLVE_FOR ! ! 1. Search through set z for a new coefficient to solve for. ! First select a candidate that is either an unconstrained ! coefficient or else a constrained coefficient that has room ! to move in the direction consistent with the sign of its dual ! vector component. Components of the dual (negative gradient) ! vector will be computed as needed. ! 2. For each candidate start the transformation to bring this ! candidate into the triangle, and then do two tests: Test size ! of new diagonal value to avoid extreme ill-conditioning, and ! the value of this new coefficient to be sure it moved in the ! expected direction. ! 3. If some coefficient passes all these conditions, set FIND = true, ! The index of the selected coefficient is J = INDEX(IZ). ! 4. If no coefficient is selected, set FIND = false. ! FIND = .FALSE. DO IZ=IZ1,IZ2 J=INDEX(IZ) ! ! Set FREE1 = true if X(J) is not at the left end-point of its ! constraint region. ! Set FREE2 = true if X(J) is not at the right end-point of its ! constraint region. ! Set FREE = true if X(J) is not at either end-point of its ! constraint region. ! FREE1 = X(J) > BND(1,J) FREE2 = X(J) < BND(2,J) FREE = FREE1 .and. FREE2 IF ( FREE ) then CALL TEST_COEF_J_FOR_DIAG!_ELT_AND_DIRECTION_OF_CHANGE else ! Compute dual coefficient W(J). W(J)=dot_product(A(NPP1:M,J),B(NPP1:M)) ! ! Can X(J) move in the direction indicated by the sign of W(J)? ! IF ( W(J) < ZERO ) then IF ( FREE1 ) & CALL TEST_COEF_J_FOR_DIAG!_ELT_AND_DIRECTION_OF_CHANGE ELSE IF ( W(J) > ZERO ) then IF ( FREE2 ) & CALL TEST_COEF_J_FOR_DIAG!_ELT_AND_DIRECTION_OF_CHANGE END IF END IF IF ( FIND ) RETURN END DO! IZ END SUBROUTINE! ( SELECT ANOTHER COEF TO SOLVE FOR ) SUBROUTINE TEST_COEF_J_FOR_DIAG!_ELT_AND_DIRECTION_OF_CHANGE ! ! The sign of W(J) is OK for J to be moved to set P. ! Begin the transformation and check new diagonal element to avoid ! near linear dependence. ! ASAVE=A(NPP1,J) call HTC (NPP1, A(1:M,J), UP) UNORM = NRM2(A(1:NSETP,J)) IF ( abs(A(NPP1,J)) > EPS * UNORM) then ! ! Column J is sufficiently independent. Copy b into Z, update Z. Z(1:M)=B(1:M) ! Compute product of transormation and updated right-hand side. NORM=A(NPP1,J); A(NPP1,J)=UP IF(ABS(NORM) > ZERO) THEN SM=DOT_PRODUCT(A(NPP1:M,J)/NORM, Z(NPP1:M))/UP Z(NPP1:M)=Z(NPP1:M)+SM*A(NPP1:M,J) A(NPP1,J)=NORM END IF IF (abs(X(J)) > ZERO) Z(1:NPP1)=Z(1:NPP1)+A(1:NPP1,J)*X(J) ! Adjust Z() as though X(J) had been reset to zero. IF ( FREE ) then FIND = .TRUE. else ! ! Solve for ZTEST ( proposed new value for X(J) ). ! Then set FIND to indicate if ZTEST has moved away from X(J) in ! the expected direction indicated by the sign of W(J). ZTEST=Z(NPP1)/A(NPP1,J) FIND = ( W(J) < ZERO .and. ZTEST < X(J) ) .or. & ( W(J) > ZERO .and. ZTEST > X(J) ) END IF END IF ! ! If J was not accepted to be moved from set Z to set P, ! restore A(NNP1,J). Failing these tests may mean the computed ! sign of W(J) is suspect, so here we set W(J) = 0. This will ! not affect subsequent computation, but cleans up the W() array. IF ( .not. FIND ) then A(NPP1,J)=ASAVE W(J)=ZERO END IF! ( .not. FIND ) END SUBROUTINE !TEST_COEF_J_FOR_DIAG!_ELT_AND_DIRECTION_OF_CHANGE SUBROUTINE MOVE_J_FROM_SET_Z_TO_SET_P ! ! The index J=index(IZ) has been selected to be moved from ! set Z to set P. Z() contains the old B() adjusted as though X(J) = 0. ! A(*,J) contains the new Householder transformation vector. B(1:M)=Z(1:M) ! INDEX(IZ)=INDEX(IZ1) INDEX(IZ1)=J IZ1=IZ1+1 NSETP=NPP1 NPP1=NPP1+1 ! The following loop can be null or not required. NORM=A(NSETP,J); A(NSETP,J)=UP IF(ABS(NORM) > ZERO) THEN DO JZ=IZ1,IZ2 JJ=INDEX(JZ) SM=DOT_PRODUCT(A(NSETP:M,J)/NORM, A(NSETP:M,JJ))/UP A(NSETP:M,JJ)=A(NSETP:M,JJ)+SM*A(NSETP:M,J) END DO A(NSETP,J)=NORM END IF ! The following loop can be null. DO L=NPP1,M A(L,J)=ZERO END DO! L ! W(J)=ZERO ! ! Solve the triangular system. Store this solution temporarily in Z(). DO I = NSETP, 1, -1 IF (I /= NSETP) Z(1:I)=Z(1:I)-A(1:I,II)*Z(I+1) II=INDEX(I) Z(I)=Z(I)/A(I,II) END DO END SUBROUTINE! ( MOVE J FROM SET Z TO SET P ) SUBROUTINE TEST_SET_P_AGAINST_CONSTRAINTS ! LOOPB: DO ! The solution obtained by solving the current set P is in the array Z(). ! ITER=ITER+1 IF (ITER > ITMAX) then IERR = 4 exit LOOPB END IF ! CALL SEE_IF_ALL_CONSTRAINED_COEFFS!_ARE_FEASIBLE ! ! The above call sets HITBND. If HITBND = true then it also sets ! ALPHA, JJ, and IBOUND. IF ( .not. HITBND ) exit LOOPB ! ! Here ALPHA will be between 0 and 1 for interpolation ! between the old X() and the new Z(). DO IP=1,NSETP L=INDEX(IP) X(L)=X(L)+ALPHA*(Z(IP)-X(L)) END DO ! I=INDEX(JJ) ! Note: The exit test is done at the end of the loop, so the loop ! will always be executed at least once. DO ! ! Modify A(*,*), B(*) and the index arrays to move coefficient I ! from set P to set Z. ! CALL MOVE_COEF_I_FROM_SET_P_TO_SET_Z ! IF (NSETP <= 0) exit LOOPB ! ! See if the remaining coefficients in set P are feasible. They should ! be because of the way ALPHA was determined. If any are infeasible ! it is due to round-off error. Any that are infeasible or on a boundary ! will be set to the boundary value and moved from set P to set Z. ! IBOUND = 0 DO JJ=1,NSETP I=INDEX(JJ) IF ( X(I) <= BND(1,I)) then IBOUND=1 EXIT ELSE IF ( X(I) >= BND(2,I)) then IBOUND=2 EXIT END IF END DO IF (IBOUND <= 0) EXIT END DO ! ! Copy B( ) into Z( ). Then solve again and loop back. Z(1:M)=B(1:M) ! DO I = NSETP, 1, -1 IF (I /= NSETP) Z(1:I)=Z(1:I)-A(1:I,II)*Z(I+1) II=INDEX(I) Z(I)=Z(I)/A(I,II) END DO END DO LOOPB ! The following loop can be null. DO IP=1,NSETP I=INDEX(IP) X(I)=Z(IP) END DO END SUBROUTINE! ( TEST SET P AGAINST CONSTRAINTS) SUBROUTINE SEE_IF_ALL_CONSTRAINED_COEFFS!_ARE_FEASIBLE ! ! See if each coefficient in set P is strictly interior to its constraint region. ! If so, set HITBND = false. ! If not, set HITBND = true, and also set ALPHA, JJ, and IBOUND. ! Then ALPHA will satisfy 0. < ALPHA <= 1. ! ALPHA=TWO DO IP=1,NSETP L=INDEX(IP) IF (Z(IP) <= BND(1,L)) then ! Z(IP) HITS LOWER BOUND LBOUND=1 ELSE IF (Z(IP) >= BND(2,L)) then ! Z(IP) HITS UPPER BOUND LBOUND=2 else LBOUND = 0 END IF ! IF ( LBOUND /= 0 ) then T=(BND(LBOUND,L)-X(L))/(Z(IP)-X(L)) IF (ALPHA > T) then ALPHA=T JJ=IP IBOUND=LBOUND END IF! ( LBOUND ) END IF! ( ALPHA > T ) END DO HITBND = abs(ALPHA - TWO) > ZERO END SUBROUTINE!( SEE IF ALL CONSTRAINED COEFFS ARE FEASIBLE ) SUBROUTINE MOVE_COEF_I_FROM_SET_P_TO_SET_Z ! X(I)=BND(IBOUND,I) IF (abs(X(I)) > ZERO .and. JJ > 0) B(1:JJ)=B(1:JJ)-A(1:JJ,I)*X(I) ! The following loop can be null. DO J = JJ+1, NSETP II=INDEX(J) INDEX(J-1)=II call ROTG (A(J-1,II),A(J,II),CC,SS) SM=A(J-1,II) ! ! The plane rotation is applied to two rows of A and the right-hand ! side. One row is moved to the scratch array S and then the updates ! are computed. The intent is for array operations to be performed ! and minimal extra data movement. One extra rotation is applied ! to column II in this approach. S=A(J-1,1:N); A(J-1,1:N)=CC*S+SS*A(J,1:N); A(J,1:N)=CC*A(J,1:N)-SS*S A(J-1,II)=SM; A(J,II)=ZERO SM=B(J-1); B(J-1)=CC*SM+SS*B(J); B(J)=CC*B(J)-SS*SM END DO ! NPP1=NSETP NSETP=NSETP-1 IZ1=IZ1-1 INDEX(IZ1)=I END SUBROUTINE! ( MOVE COEF I FROM SET P TO SET Z ) SUBROUTINE TERMINATION ! IF (IERR <= 0) then ! ! Compute the norm of the residual vector. SM=ZERO IF (NPP1 <= M) then SM=NRM2(B(NPP1:M)) else W(1:N)=ZERO END IF RNORM=SM END IF! ( IERR...) END SUBROUTINE! ( TERMINATION ) SUBROUTINE ROTG(SA,SB,C,S) ! REAL(KIND(ONE)) SA,SB,C,S,ROE,SCALE,R ! ROE = SB IF( ABS(SA) .GT. ABS(SB) ) ROE = SA SCALE = ABS(SA) + ABS(SB) IF( SCALE .LE. ZERO ) THEN C = ONE S = ZERO RETURN END IF R = SCALE*SQRT((SA/SCALE)**2 + (SB/SCALE)**2) IF(ROE < ZERO)R=-R C = SA/R S = SB/R SA = R RETURN END SUBROUTINE !ROTG REAL(KIND(ONE)) FUNCTION NRM2 (X) ! ! NRM2 returns the Euclidean norm of a vector via the function ! name, so that ! ! NRM2 := sqrt( x'*x ) ! REAL(KIND(ONE)) ABSXI, X(:), NORM, SCALE, SSQ INTEGER N, IX N=SIZE(X) IF( N < 1)THEN NORM = ZERO ELSE IF( N == 1 )THEN NORM = ABS( X( 1 ) ) ELSE SCALE = ZERO SSQ = ONE ! DO IX = 1, N ABSXI = ABS( X( IX ) ) IF(ABSXI > ZERO )THEN IF( SCALE < ABSXI )THEN SSQ = ONE + SSQ*( SCALE/ABSXI )**2 SCALE = ABSXI ELSE SSQ = SSQ + ( ABSXI/SCALE )**2 END IF END IF END DO NORM = SCALE * SQRT( SSQ ) END IF ! NRM2 = NORM RETURN END FUNCTION SUBROUTINE HTC (P, U, UP) ! ! Construct a Householder transormation. INTEGER P REAL(KIND(ONE)) U(:) REAL UP, VNORM VNORM=NRM2(U(P:SIZE(U))) IF(U(P) > ZERO) VNORM=-VNORM UP=U(P)-VNORM U(P)=VNORM END SUBROUTINE ! HTC END SUBROUTINE ! BVLS program PROG7 ! ! DEMONSTRATE THE USE OF THE SUBROUTINE BVLS FOR LEAST ! SQUARES SOLVING WITH BOUNDS ON THE VARIABLES. ! ! The original version of this code was developed by ! Charles L. Lawson and Richard J. Hanson and published in the book ! "SOLVING LEAST SQUARES PROBLEMS." REVISED APRIL, 1995 to accompany ! reprinting of the book by SIAM. IMPLICIT NONE INTERFACE SUBROUTINE BVLS (A, B, BND, X, RNORM, NSETP, W, INDEX, IERR) REAL(KIND(1E0)) A(:,:), B(:), BND(:,:), X(:), RNORM, W(:) INTEGER NSETP, INDEX(:), IERR END SUBROUTINE END INTERFACE ! Test driver for BVLS. Bounded Variables Least Squares. ! C.L.Lawson, & R.J.Hanson, Jet Propulsion Laboratory, 1973 June 12 ! Changes made in September 1982, June 1986, October 1987. ! Conversion made to Fortran 90 by R. J. Hanson, April 1995. ! ------------------------------------------------------------------ ! Subprograms referenced: RANDOM_NUMBER, BVLS ! ------------------------------------------------------------------ integer, parameter :: MM=10, NN=10, MXCASE = 6, JSTEP=5 INTEGER I, J, ICASE, IERR, m, n, nsetp, j1, j2 integer MTAB(MXCASE), NTAB(MXCASE) real(kind(1E0)) UNBTAB(MXCASE), BNDTAB(2,NN,MXCASE) real(kind(1E0)) A(MM,NN),B(MM),X(NN),W(NN) real(kind(1E0)) A2(MM,NN),B2(MM), R(MM), D(NN),BND(2,NN) real(kind(1E0)) RNORM, RNORM2, UNBND integer INDEX(NN) data MTAB / 2, 2, 4, 5, 10, 6 / data NTAB / 2, 4, 2, 10, 5, 4 / data UNBTAB / 5 * 1.0E6, 999.0E0 / data ((BNDTAB(I,J,1),I=1,2),J=1,2)/ 1.,2., 3.,4. / data ((BNDTAB(I,J,2),I=1,2),J=1,4)/ 0,10, 0,10, 0,10, 0,10/ data ((BNDTAB(I,J,3),I=1,2),J=1,2)/ 0,100, -100,100/ data ((BNDTAB(I,J,4),I=1,2),J=1,10)/& 0,0, -.3994E0,-.3994E0, -1,1, -.3E0,-.2E0, 21,22,& -4,-3, 45,46, 100,101, 1.E6,1.E6, -1,1/ data ((BNDTAB(I,J,5),I=1,2),J=1,5)/& 0,1, -1,0, 0,1, .3E0,.4E0, .048E0,.049E0/ data ((BNDTAB(I,J,6),I=1,2),J=1,4)/& -100.,100., 999.,999., 999.,999., 999.,999. / ! ------------------------------------------------------------------ write(*,1002) DO ICASE = 1,MXCASE M = MTAB(ICASE) N = NTAB(ICASE) UNBND = UNBTAB(ICASE) DO J = 1,N BND(1,J) = BNDTAB(1,J,ICASE) BND(2,J) = BNDTAB(2,J,ICASE) END DO where(bnd(1,1:N) == UNBND) bnd(1,1:n)=-huge(1e0) where(bnd(2,1:N) == UNBND) bnd(2,1:n)= huge(1e0) write(*,'(1x/////1x,a,i3/1x)') 'Case No.',ICASE write(*,'(1x,a,i5,a,i5,a,g17.5)')& 'M =',M,', N =',N,', UNBND =',UNBND write(*,'(1X/'' Bounds ='')') DO J1 = 1, N, JSTEP J2 = MIN(J1 - 1 + JSTEP, N) write(*,*) ' ' write(*,1001) (BND(1,J),J=J1,J2) write(*,1001) (BND(2,J),J=J1,J2) END DO call random_number (b(1:m)) DO J=1,N call random_number (A(1:M,J)) END DO write(*,'(1X/'' A(,) ='')') DO J1 = 1, N, JSTEP J2 = MIN(J1 - 1 + JSTEP, N) write(*,*) ' ' DO I = 1,M write(*,1001) (A(I,J),J=J1,J2) END DO END DO ! B2(1:M)=B(1:M) A2(1:M,1:N)=A(1:M,1:N) write(*,'(1X/'' B() =''/1X)') write(*,1001) (B(I),I=1,M) ! call BVLS (A2(1:M,1:N), B2,BND,X,RNORM,NSETP,W,INDEX, IERR) ! IF(IERR > 0) THEN WRITE(*,'(1X, "ABNORMAL ERROR FLAG, IERR = ")') IERR STOP END IF write(*,'(1x/1x,a,i4)') & 'After BVLS: No. of components not at constraints =',NSETP write(*,'(1X/'' Solution vector, X() =''/1X)') write(*,1001) (X(J),J=1,N) R(1:M)=B(1:M)-matmul(A(1:M,1:N),X(1:N)) ! RNORM2=sqrt(dot_product(R(1:M),R(1:M))) ! D(1:N)=matmul(R(1:M),A(1:M,1:N)) write(*,'(1X/'' R = B - A*X Computed by the driver:''/1X)') write(*,1001) (R(I),I=1,M) write(*,'(1x,a,g17.5)') 'RNORM2 computed by the driver =',RNORM2 write(*,'(1x,a,g17.5)') 'RNORM computed by BVLS = ',RNORM ! write(*,'(1X/'' W = (A**T)*R Computed by the driver:''/1X)') write(*,1001) (D(J),J=1,N) write(*,'(1X/'' Dual vector from BVLS, W() =''/1X)') write(*,1001) (W(J),J=1,N) ! END DO ! ICASE ! 1001 FORMAT(1X,5G14.5) 1002 FORMAT(& ' TESTBVLS.. Test driver for BVLS,',& ' Bounded Variables Least Squares.'& /' If the algorithm succeeds the solution vector, X(),'/& ' and the dual vector, W(),'/& ' should be related as follows:'/& ' X(i) not at a bound => W(i) = 0'/& ' X(i) at its lower bound => W(i) .le. 0'/& ' X(i) at its upper bound => W(i) .ge. 0'/& ' except that if an upper bound and lower bound are equal then'/& ' the corresponding X(i) must take that value and W(i) may have'/& ' any value.'/1X) end program PROG7