SUBROUTINE SVA (A,MDA,M,N,MDATA,B,SING,NAMES,ISCALE,D) C C INPUT.. C C A(*,*),MDA,M,N THE ARRAY A(*,*) INITIALLY CONTAINS THE M BY C N MATRIX A OF THE LEAST SQUARES PROBLEM AX = C B. EITHER M.GE.N OR M.LT.N IS PERMITTED. THE C FIRST DIMENSIONING PARAMETER OF THE ARRAY C A(*,*) IS MDA, WHICH MUST SATISFY MDA .GE. C MAX(M,N). THE CONDITION MDA.LT.MAX(M,N) IS C CONSIDERED AN ERROR. C C MDATA THE NUMBER OF ROWS IN THE ORIGINAL LEAST C SQUARES PROBLEM, WHICH MAY HAVE HAD MORE ROWS C THAN (A B). THIS NUMBER IS USED ONLY IN C COMPUTING THE NUMBERS SG(J), J=0,...,K. C C B(*) THE ARRAY B(*) INITIALLY CONTAINS THE M-VECTOR C B OF THE LEAST SQUARES PROBLEM AX = B. C C SING(*) THE LOCATIONS, (SING(I), I=1,...,3*N), ARE C USED AS TEMPORARY WORKING SPACE. C C NAMES(*) IN LOCATION NAMES(I) THE USER MAY STORE AN C ALPHA-NUMERIC NAME AS AN IDENTIFIER FOR THE C I-TH COMPONENT OF THE SOLUTION VECTOR, (I = C 1,...,N). THESE NAMES WILL BE PRINTED WITH C THE APPROPRIATE COMPONENTS OF THE V MATRIX C AND WITH THE CANDIDATE SOLUTIONS. EACH NAME C IS PRINTED WITH AN A6 FORMAT. IF NAMES(1) HAS C THE SAME CONTENTS AS THE QUANTITY IBLANK, C DEFINED BY C C DATA IBLANK/' '/ C C THEN THE REMAINING WORDS OF THE ARRAY NAMES(*) C WILL BE IGNORED AND NO IDENTIFYING NAMES WILL C BE PRINTED. C C ISCALE,D(*) THE USER SETS ISCALE = 1,2, OR 3. C IF ISCALE = 1, THE SUBROUTINE FUNCTIONS AS C THOUGH THE SCALING MATRIX D WERE AN N BY N C IDENTITY MATRIX. IN THIS CASE THE SUBROUTINE C MAKES NO REFERENCE TO THE ARRAY D(*). C C IF ISCALE = 2, THE SUBROUTINE COMPUTES THE C LENGTH OF THE J-TH COLUMN VECTOR OF A AND SETS C D(J) EQUAL TO THE RECIPROCAL OF THIS LENGTH, C IF NONZERO, AND D(J) = 1, IF ZERO. IN THIS C CASE ALL NONZERO COLUMNS OF THE SCALED MATRIX C C = AD WILL HAVE UNIT EUCLIDEAN LENGTH. C C IF ISCALE = 3, THE SUBROUTINE WILL USE THE C USER-SUPPLIED VALUES D(J),J = 1,...,N AS THE C DIAGONAL ELEMENTS OF THE DIAGONAL SCALING C MATRIX D. IN THIS CASE THE USER MUST ASSIGN C VALUES TO D(J), J = 1,...,N, AND THE C SUBROUTINE WILL NOT MODIFY THESE VALUES. FOR C EXAMPLE, THE USER MIGHT SET D(J) EQUAL TO THE C A PRIORI STANDARD DEVIATION OF THE J-TH C COMPONENT OF THE SOLUTION VECTOR IF SUCH C INFORMATION IS KNOWN. AS A FURTHER EXAMPLE, C THE USER MIGHT SET D(J) EQUAL TO THE C RECIPROCAL OF THE MEAN OF THE M COMPONENTS OF C THE J-TH COLUMN VECTOR S UNCERTAINTY. THIS C TYPE OF SCALING HAS THE EFFECT OF MAKING THE C UNCERTAINTY OF EVERY COLUMN THE SAME MAGNITUDE C IN THE MATRIX C. C C OUTPUT.. C C A(*,*) ON OUTPUT, THE J-TH CANDIDATE SOLUTION WILL BE C STORED IN THE FIRST N LOCATIONS OF THE J-TH C COLUMN OF A(*,*) FOR J = 1,...,MIN(M,N). C C B(*) ON OUTPUT THE M-VECTOR G IS STORED IN THE C ARRAY B(*). C C SING(*) ON OUTPUT THE SINGULAR VALUES OF THE SCALED C MATRIX C ARE STORED IN SING(I),I = 1,..., C MIN(M,N). C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER(I-N) C DIMENSION A(MDA,N), B(M), SING(01),D(N) C SING(3*N) CHARACTER*6 NAMES(N) LOGICAL TEST C DZERO=0.D0 ONE=1. ZERO=0. IF (M.LE.0 .OR. N.LE.0) RETURN IF(MAX(M,N).GT.MDA) RETURN NP1=N+1 WRITE (6,270) WRITE (6,260) M,N,MDATA GO TO (60,10,10), ISCALE C C APPLY COLUMN SCALING TO A C 10 DO 50 J=1,N A1=D(J) GO TO (20,20,40), ISCALE 20 SB=DZERO 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 40 DO 50 I=1,M 50 A(I,J)=A(I,J)*A1 WRITE (6,280) ISCALE,(D(J),J=1,N) GO TO 70 60 CONTINUE WRITE (6,290) 70 CONTINUE C C OBTAIN SING. VALUE DECOMP. OF SCALED MATRIX. C C********************************************************** CALL SVDRS (A,MDA,M,N,B,1,1,SING) C********************************************************** C C PRINT THE V MATRIX. C CALL MFEOUT (A,MDA,N,N,NAMES,1) IF (ISCALE.EQ.1) GO TO 90 C C REPLACE V BY D*V IN THE ARRAY A() C DO 80 I=1,N DO 80 J=1,N 80 A(I,J)=D(I)*A(I,J) C C G NOW IN B ARRAY. V NOW IN A ARRAY. C C C OBTAIN SUMMARY OUTPUT. C 90 CONTINUE WRITE (6,220) C C COMPUTE CUMULATIVE SUMS OF SQUARES OF COMPONENTS OF C G AND STORE THEM IN SING(I), I=MINMN+1,...,2*MINMN+1 C SB=DZERO MINMN=MIN(M,N) MINMN1=MINMN+1 IF (M.EQ.MINMN) GO TO 110 DO 100 I=MINMN1,M 100 SB=SB+B(I)**2 110 SING(2*MINMN+1)=SB DO 120 JJ=1,MINMN J=MINMN+1-JJ SB=SB+B(J)**2 JS=MINMN+J 120 SING(JS)=SB A3=SING(MINMN+1) A4=SQRT(A3/DBLE(MAX(1,MDATA))) WRITE (6,230) A3,A4 C NSOL=0 C C C DO 160 K=1,MINMN IF (SING(K).EQ.ZERO) GO TO 130 NSOL=K PI=B(K)/SING(K) A1=ONE/SING(K) A2=B(K)**2 A3=SING(MINMN1+K) A4=SQRT(A3/DBLE(MAX(1,MDATA-K))) TEST=SING(K).GE.100..OR.SING(K).LT..001 IF (TEST) WRITE (6,240) K,SING(K),PI,A1,B(K),A2,A3,A4 IF (.NOT.TEST) WRITE(6,250)K,SING(K),PI,A1,B(K),A2,A3,A4 GO TO 140 130 WRITE (6,240) K,SING(K) PI=ZERO 140 DO 150 I=1,N A(I,K)=A(I,K)*PI 150 IF (K.GT.1) A(I,K)=A(I,K)+A(I,K-1) 160 CONTINUE C C COMPUTE AND PRINT VALUES OF YNORM AND RNORM. C WRITE (6,300) J=0 YSQ=ZERO GO TO 180 170 J=J+1 YSQ=YSQ+(B(J)/SING(J))**2 180 YNORM=SQRT(YSQ) JS=MINMN1+J RNORM=SQRT(SING(JS)) YL=-1000. IF (YNORM.GT.0.) YL=LOG10(YNORM) RL=-1000. IF (RNORM.GT.0.) RL=LOG10(RNORM) WRITE (6,310) J,YNORM,RNORM,YL,RL IF (J.LT.NSOL) GO TO 170 C C COMPUTE VALUES OF XNORM AND RNORM FOR A SEQUENCE OF VALUES OF C THE LEVENBERG-MARQUARDT PARAMETER. C IF (SING(1).EQ.ZERO) GO TO 210 EL=LOG10(SING(1))+ONE EL2=LOG10(SING(NSOL))-ONE DEL=(EL2-EL)/20. TEN=10. ALN10=LOG(TEN) WRITE (6,320) DO 200 IE=1,21 C COMPUTE ALAMB=10.**EL ALAMB=EXP(ALN10*EL) YS=0. JS=MINMN1+NSOL RS=SING(JS) 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=-1000. IF (RNORM.GT.ZERO) RL=LOG10(RNORM) YL=-1000. IF (YNORM.GT.ZERO) YL=LOG10(YNORM) WRITE (6,330) ALAMB,YNORM,RNORM,EL,YL,RL EL=EL+DEL 200 CONTINUE C C PRINT CANDIDATE SOLUTIONS. C 210 IF (NSOL.GE.1) CALL MFEOUT (A,MDA,N,NSOL,NAMES,2) RETURN 220 FORMAT (42H0 INDEX SING. VALUE P COEF ,48HRECIP. S. 1V. G COEF G**2 ,39H C.S.S. 2 N.S.R.C.S.S.) 230 FORMAT (1H ,5X,1H0,88X,1PE15.4,1PE17.4) 240 FORMAT (1H ,I6,E12.4,1P(E15.4,4X,E15.4,4X,E15.4,4X,E15.4,4X,E15.4, 12X,E15.4)) 250 FORMAT (1H ,I6,F12.4,1P(E15.4,4X,E15.4,4X,E15.4,4X,E15.4,4X,E15.4, 12X,E15.4)) 260 FORMAT (5H0M = ,I6,8H, N = ,I4,12H, MDATA = ,I8) 270 FORMAT (45H0SINGULAR VALUE ANALYSIS OF THE LEAST SQUARES,42H PROBL 1EM, A*X=B, SCALED AS (A*D)*Y=B .) 280 FORMAT (19H0SCALING OPTION NO.,I2,18H. D IS A DIAGONAL,46H MATRIX 1 WITH THE FOLLOWING DIAGONAL ELEMENTS../(5X,10E12.4)) 290 FORMAT (50H0SCALING OPTION NO. 1. D IS THE IDENTITY MATRIX./1X) 300 FORMAT (6H0INDEX,13X,28H YNORM RNORM,14X,28H LOG1 10(YNORM) LOG10(RNORM)/1X) 310 FORMAT (1H ,I4,14X,2E14.5,14X,2F14.5) 320 FORMAT (54H0NORMS OF SOLUTION AND RESIDUAL VECTORS FOR A RANGE OF, 144H VALUES OF THE LEVENBERG-MARQUARDT PARAMETER,9H, LAMBDA./1H0,4X 2,42H LAMBDA YNORM RNORM,42H LOG10(LAMBDA) 3LOG10(YNORM) LOG10(RNORM)) 330 FORMAT (5X,3E14.5,3F14.5) END SUBROUTINE SVDRS (A,MDA,MM,NN,B,MDB,NB,S) C C INPUT.. C C A(*,*),MDA,M,N THE ARRAY A(*,*) IS DOUBLY SUBSCRIPTED WITH C FIRST DIMENSIONING PARAMETER EQUAL TO MDA. C THE ARRAY A(*,*) INITIALLY CONTAINS THE M BY C N MATRIX A. EITHER M.GE.N OR M.LT.N IS C PERMITTED. THE FIRST DIMENSIONING PARAMETER C MDA MUST SATISFY MDA.GE.MAX(M,N). THE C CONDITION MDA.LT.MAX(M,N) IS CONSIDERED AN C ERROR. C C B(*),MDB,NB NB DENOTES THE NUMBER OF COLUMN VECTORS IN THE C MATRIX B. IF NB = 0 THE ARRAY B(*) WILL NOT C BE REFERENCED BY THIS SUBROUTINE NOR WILL ITS C CONTENTS BE MODIFIED. IF NB.GE.2, THE ARRAY B C SHOULD BE DOUBLE SUBSCRIPTED WITH FIRST C DIMENSIONING PARAMETER EQUAL TO MDB.GE.M. IF C NB = 1, THEN B WILL BE USED AS A SINGLY C SUBSCRIPTED ARRAY OF LENGTH M. IN THIS LATTER C CASE THE VALUE OF MDB IS ARBITRARY BUT IT C SHOULD BE SET TO A DEFINED VALUE SUCH AS C MDB = M. THE CONTENTS OF THE ARRAY B(*) MUST C CONTAIN THE MATRIX B. THE CONDITION NB.GT.1 C AND. MDB.LT.M IS CONSIDERED AN ERROR. C C S(*) THE ARRAY S(*) IS USED AS 3*N LOCATIONS OF C TEMPORARY WORKING SPACE BY THE SUBROUTINE. C C OUTPUT.. C C A(*,*) ON OUTPUT THE ARRAY A(*,*) CONTAINS THE N BY N C ORTHOGONAL MATRIX V. C C B(*) ON OUTPUT THE ARRAY B(*),(OR B(*,*)), CONTAINS C THE PRODUCT OF THE TRANSPOSE OF THE M BY M C MATRIX U AND THE RIGHTHAND SIDE B. C C S(*) ON OUTPUT THE FIRST N LOCATIONS OF S(*) C CONTAIN THE DIAGONAL TERMS OF THE M BY N C DIAGONAL MATRIX S. THESE NONNEGATIVE TERMS C ARE ORDERED FROM LARGEST TO SMALLEST. IF C M.LT.N, ENTRIES M+1,...,N OF S(*) ARE ZERO. C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER(I-N) C DIMENSION A(MDA,NN),B(MDB,NB),S(NN,3) C ZERO=0. ONE=1. C C BEGIN.. SPECIAL FOR ZERO ROWS AND COLS. C C PACK THE NONZERO COLS TO THE LEFT C N=NN IF (N.LE.0.OR.MM.LE.0) RETURN IF(MAX(MM,NN).GT.MDA) RETURN C IF(NB.GT.1.AND.MM.GT.MDB) RETURN C J=N 10 CONTINUE DO 20 I=1,MM IF (A(I,J)) 50,20,50 20 CONTINUE C C COL J IS ZERO. EXCHANGE IT WITH COL N. C IF (J.EQ.N) GO TO 40 DO 30 I=1,MM 30 A(I,J)=A(I,N) 40 CONTINUE 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=MM 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,S(J,3),A(J+1,1),MDA,1,M-J) 170 CONTINUE C C COPY THE BIDIAGONAL MATRIX INTO THE ARRAY S() FOR QRBD. C IF (N.EQ.1) GO TO 190 DO 180 J=2,N S(J,1)=A(J,J) 180 S(J,2)=A(J-1,J) 190 S(1,1)=A(1,1) C NS=N IF (M.GE.N) GO TO 200 NS=M+1 S(NS,1)=ZERO S(NS,2)=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,S(I,3),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,1),S(1,2),NS,A,MDA,N,B,MDB,NB) C GO TO (240,310), IPASS 240 CONTINUE IF (NS.GE.N) GO TO 260 NSP1=NS+1 DO 250 J=NSP1,N 250 S(J,1)=ZERO 260 CONTINUE IF (N.EQ.NN) RETURN NP1=N+1 C MOVE RECORD OF PERMUTATIONS C AND STORE ZEROS DO 280 J=NP1,NN S(J,1)=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,NN I=S(K,1) S(K,1)=ZERO DO 290 J=1,NN 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 310 CONTINUE C RETURN END C SUBROUTINE QRBD (IPASS,Q,E,NN,V,MDV,NRV,C,MDC,NCC) C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 C TO APPEAR IN ^SOLVING LEAST SQUARES PROBLEMS^, PRENTICE-HALL, 1974 C C MODIFIED AT SANDIA LABS, MAY 1977, TO INCREASE EFFICIENCY C THROUGH THE USE OF THE BASIC LINEAR ALGEBRA MODULE SROT , C WHICH PERFORMS THE GIVENS ROTATION APPLICATION STEP. C C QR ALGORITHM FOR SINGULAR VALUES OF A BIDIAGONAL MATRIX. 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 SUBROUTINE QRBD (IPASS,Q,E,NN,V,MDV,NRV,C,MDC,NCC) C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER(I-N) C LOGICAL WNTV ,HAVERS,FAIL DIMENSION Q(NN),E(NN),V(MDV,NN),C(MDC,NCC) ZERO=0. ONE=1. TWO=2. 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)) 50,25,50 C C SINCE Q(K) IS SMALL WE WILL MAKE A SPECIAL PASS TO C TRANSFORM E(K) TO ZERO. C 25 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 C THE FOLLOWING TWO STATEMENTS HAVE BEEN REPLACED BY A CALL TO SROT. C DO 30 J=1,NRV C 30 CALL G2 (CS,SN,V(J,I),V(J,K)) CALL SROT (NRV,V(1,I),1,V(1,K),1,CS,SN) 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)) 55,100,55 55 IF(DIFF(DNORM+Q(L-1),DNORM)) 60,70,60 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)) 75,100,75 75 CALL G1 (Q(I),F,CS,SN,Q(I)) IF (.NOT.HAVERS) GO TO 90 C THE FOLLOWING TWO STATEMENTS HAVE BEEN REPLACED BY A CALL TO SROT. C DO 80 J=1,NCC C 80 CALL G2 (CS,SN,C(I,J),C(L-1,J)) CALL SROT (NCC,C(I,1),MDC,C(L-1,1),MDC,CS,SN) 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.LT.ZERO) GO TO 110 T=F+G GO TO 120 110 T=F-G 120 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 (.NOT.WNTV) GO TO 140 C C ACCUMULATE ROTATIONS (FROM THE RIGHT) IN ^V^ C THE FOLLOWING TWO STATEMENTS HAVE BEEN REPLACED BY A CALL TO SROT. C DO 130 J=1,NRV C 130 CALL G2 (CS,SN,V(J,I-1),V(J,I)) CALL SROT (NRV,V(1,I-1),1,V(1,I),1,CS,SN) 140 CALL G1 (F,H,CS,SN,Q(I-1)) F=CS*G+SN*Y X=-SN*G+CS*Y IF (.NOT.HAVERS) GO TO 160 C THE FOLLOWING TWO STATEMENTS HAVE BEEN REPLACED BY A CALL TO SROT. C DO 150 J=1,NCC C 150 CALL G2 (CS,SN,C(I-1,J),C(I,J)) CALL SROT (NCC,C(I-1,1),MDC,C(I,1),MDC,CS,SN) 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 FAIL=.TRUE. C .. C CUTOFF FOR CONVERGENCE FAILURE. ^NQRS^ WILL BE 2*N USUALLY. 170 IF (Z.GE.ZERO) GO TO 190 Q(K)=-Z IF (.NOT.WNTV) GO TO 190 DO 180 J=1,NRV 180 V(J,K)=-V(J,K) 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 (.NOT.HAVERS) GO TO 250 DO 240 J=1,NCC T=C(I-1,J) C(I-1,J)=C(K,J) 240 C(K,J)=T 250 IF (.NOT.WNTV) GO TO 270 DO 260 J=1,NRV T=V(J,I-1) V(J,I-1)=V(J,K) 260 V(J,K)=T 270 CONTINUE C END OF ORDERING ALGORITHM. C IF (FAIL) IPASS=2 RETURN END SUBROUTINE G1(A,B,COS,SIN,SIG) C C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 C TO APPEAR IN ^SOLVING LEAST SQUARES PROBLEMS^, PRENTICE HALL, 1974 C C COMPUTE ORTHOGONAL ROTAION MATRIX C COMPUTE MATRIX (C,S) SO THAT (C,S)(A) = (SQRT(A**2+B**2) C (-S,C) SO THAT (-S,C)(B) = (0) C COMPUTE SIG = SQRT(A**2+B**2) C SIG IS COMPUTED TO ALLOW FOR THE POSSIBILITY THAT C SIG MAY BE IN THE SAME LOCATION AS A OR B C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER(I-N) C ZERO = 0.0 ONE = 1.0 IF(ABS(A).LE.ABS(B)) GO TO 10 XR = B/A YR = SQRT(ONE+XR**2) COS = SIGN(ONE/YR,A) SIN = COS*XR SIG = ABS(A)*YR RETURN 10 IF(B) 20,30,20 20 XR = A/B YR = SQRT(ONE+XR**2) SIN = SIGN(ONE/YR,B) COS = SIN*XR SIG = ABS(B)*YR RETURN 30 SIG = ZERO COS = ZERO SIN = ONE RETURN END SUBROUTINE G2(COS,SIN,X,Y) C C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1972 DEC 15 C TO APPEAR IN ^SOLVING LEAST SQUARES PROBLEMS^, PRENTICE-HALL, 1974 C APPLY THE ROTATION COMPUTED BY G1 TO (X,Y) C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER(I-N) C XR = COS*X + SIN*Y Y = -SIN*X + COS*Y X = XR RETURN END FUNCTION DIFF(X,Y) C C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUNE 7 C TO APPEAR IN ^SOLVING LEAST SQUARES PROBLEMS^, PRENTICE-HALL, 1974 C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER(I-N) C DIFF = X - Y RETURN END SUBROUTINE LDP(G,MDG,M,N,H,X,XNORM,W,INDEX,MODE) C C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1974 MAR 1 C TO APPEAR IN ^SOLVING LEAST SQUARES PROBLEMS^, PRENTICE HALL, 1974 C C LEAST DISTANCE PROGRAMMING C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER(I-N) C INTEGER INDEX(M) DIMENSION G(MDG,N),H(M),X(N),W(1) ZERO = 0.0 ONE = 1.0 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 A 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 C COPY G(T) INTO FORST 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) IF = IW + 1 C C STORE N ZEROS FOLLOWED BY A ONE INTO F. C 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 IDUAL = IY + M C CALL NNLS(W,NP1,NP1,M,W(IF),W(IY),RNORM,W(IDUAL),W(IZ),INDEX, 1 MODE) C C USE THE FOLLOWING RETURN IF UNSUCCESFUL IN NNLS C 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 C HERE WE ARE USING THE SOLUTION VECTOR Y C 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 C HERE WE ARE USING THE SOLUTION VECTOR Y C 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 C SUCCESSFUL RETURN C 110 MODE = 1 RETURN C C ERROR RETURN C 120 MODE = 2 RETURN 130 MODE = 4 RETURN END SUBROUTINE NNLS(A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE) C C C.L.LAWSON AND R.J.HANSON,JET PROPULSION LABORATORY, 1973 JUNE 15 C TO APPEAR IN ^SOLVING LEAST SQUARES PROBLEMS^, PRENTICE HALL,1974 C C *****NON-NEGATIVE LEAST SQUARES***** C C GIVEN AN M BY N MATRIX, A, AND AN M VECTOR, B, COMPUTE AN N C VECTOR, X, WHICH SOLVES THE LEAST SQUARES PROBLEM C C A*X = B SUBJECT TO X >= 0, AND 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 THE PRODUCT C MATRIX, Q*A, WHERE Q IS A M BY M ORTHOGONAL C MATRIX GENERATED IMPLICITLY BY THIS SUBROUTINE. C B() ON ENTRY B() CONTAINS THE M-VECTOR, B. ON EXIT B() C CONTAINS 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 EUCLIDIAN NORN OF THE RESIDUAL C VECTOR C W() AN N-ARRAY OF WORKING SPACE. ON EXIT W() WILL CONTAIN THE C DUAL SOLUTION VECTOR. W WILL SATISFY W(I) = 0. FOR ALL C I IN THE SET P AND W(I) .LE. 0 FOR ALL I IN THE 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) = SETP 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 IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER(I-N) C DIMENSION A(MDA,N),B(M),X(N),W(N),ZZ(M),INDEX(N) ZERO = 0.0 ONE = 1.0 TWO = 2.0 FACTOR = 0.01 C MODE = 1 IF(M.GT.0.AND.N.GT.0) GO TO 10 MODE = 2 RETURN 10 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 C *****MAIN LOOP BEGINS HERE***** C 30 CONTINUE C C QUIT IF ALL COEFFICIENTS ARE ALREADY IN THE SOLUTION OF IF M C COLUMNS 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) 50 W(J) = SM C C FIND LARGEST POSITIVE W(J) C 60 WMAX = ZERO DO 70 IZ=IZ1,IZ2 J = INDEX(IZ) IF(W(J).LE.WMAX) GO TO 70 WMAX = W(J) IZMAX = IZ 70 CONTINUE C C IF WMAX .LE. 0 GO TO TERMINATION C THIS INDICATES SATISFACTION OF THE KUHN-TUCKER CONDITIONS C IF(WMAX) 350,350,80 80 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.EQ.0) GO TO 100 DO 90 L=1,NSETP 90 UNORM = UNORM + A(L,J)**2 100 UNORM = SQRT(UNORM) IF(DIFF(UNORM+ABS(A(NPP1,J))*FACTOR,UNORM)) 130,130,110 C C SOL J IS SUFFICIENTLY INDEPENDENT. COPY B INTO ZZ, UPDATE ZZ AND C SOLVE FOR ZTEST ( = PROPOSED NEW VALUE FOR X(J)). C 110 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) 130,130,140 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 130 A(NPP1,J) = ASAVE W(J) = ZERO GO TO 60 C C THE INDEX J = INDEX(IZ) HAS BEEN SELECTED TO BE MOVED FROM THE 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 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.GT.IZ2) GO TO 170 DO 160 JZ=IZ1,IZ2 JJ = INDEX(JZ) 160 CALL H12(2,NSETP,NPP1,M,A(1,J),1,UP,A(1,JJ),1,MDA,1) 170 CONTINUE C IF(NSETP.EQ.M) GO TO 190 DO 180 L=NPP1,M 180 A(L,J) = ZERO 190 CONTINUE C W(J) = ZERO C C SOLVE THE TRIANGULAR SYSTEM C STORE THE SOLUTION TEMPORARILY IN ZZ() C ASSIGN 200 TO NEXT GO TO 400 200 CONTINUE C C *****SECONDARY LOOP BEGINS HERE***** C C ITERATION COUNTER C 210 ITER = ITER + 1 IF(ITER.LE.ITMAX) GO TO 220 MODE = 3 WRITE(6,500) GO TO 350 220 CONTINUE 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)) 230,230,240 C 230 T = -X(L)/(ZZ(IP)-X(L)) IF(ALPHA.LE.T) GO TO 240 ALPHA = T JJ = IP 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) 250 X(L) = X(L) + ALPHA*(ZZ(IP)-X(L)) 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 X(I) = ZERO C IF(JJ.EQ.NSETP) GO TO 290 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) CALL G2(CC,SS,A(J-1,L),A(J,L)) 270 CONTINUE 280 CALL G2(CC,SS,B(J-1),B(J)) 290 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 NON-POSITIVE 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)) 260,260,300 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) ASSIGN 320 TO NEXT GO TO 400 320 CONTINUE GO TO 210 C C *****END OF SECONDARY LOOP***** C 330 DO 340 IP=1,NSETP I = INDEX(IP) 340 X(I) = ZZ(IP) C C ALL NEW COEFFICIENTS ARE POSITIVE. LOOP BACK TO BEGINNING. C 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 SM = ZERO IF(NPP1.GT.M) GO TO 370 DO 360 I=NPP1,M 360 SM = SM + B(I)**2 GO TO 390 370 DO 380 J=1,N 380 W(J) = ZERO 390 RNORM = SQRT(SM) RETURN C 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 DO 430 L=1,NSETP IP = NSETP + 1 - L IF(L.EQ.1) GO TO 420 DO 410 II=1,IP 410 ZZ(II) = ZZ(II) - A(II,JJ)*ZZ(IP+1) 420 JJ = INDEX(IP) 430 ZZ(IP) = ZZ(IP)/A(IP,JJ) GO TO NEXT,(200,320) 500 FORMAT(35H0 NNLS QUITTING ON ITERATION COUNT ) END SUBROUTINE SEQHT(A,N,N1,W,B,T,L) C C ALGORITHM SEQHT (SEQUENTIAL HOUSEHOLDER TRIANGULARIZATION) C FROM LAWSON,C.L. AND R.J.HANSON (1973) C SOLVING LEAST SQUARES PROBLEMS C PRENTICE HALL,INC., NEW JERSEY C C W(N) VECTOR OF COEFFICIENTS FOR PRESENT CASE C B KNOWN CONSTANT FOR THE PRESENT CASE C T ULP (Q=I+U*U(T)/(S*ULP)) C C ON THE FINAL CALL Q[A B] = A = R D C 0 E C 0 0 C (DIMENSIONED N1 BY N 1 1) (R IS UPPER TRIANGULAR) C C LS R X = D N1 >= N+2 C 0 E C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER(I-N) C DIMENSION A(N1,1),W(N) INTEGER P C C INITIALLY (FIRST CALL) L IS SET TO ZERO C STEP 1 L IS INITIALLY SET TO 0 C STEP 2 SUBROUTINE IS CALLED NCASES TIMES C STEP 3 M(I) = 1 ALWAYS C P = L + 1 C C STEP 4 C DO 10 I=1,N 10 A(P,I) = W(I) A(P,N+1) = B C C STEP 5 C NP1 = MIN(N+1,P-1) IF(NP1.EQ.0) GO TO 30 DO 20 I=1,NP1 20 CALL H12(1,I,MAX(I+1,L+1),P,A(1,I),1,T,A(1,I+1),1,N1,N+1-I) 30 CONTINUE C C STEP 6 C L = MIN(N+1,P) C C STEP 7 UPON THE FINAL CALL HOUSEHOLDER TRIANGULARIZATION IS C COMPLETE C RETURN END SUBROUTINE SROT (N,SX,INCX,SY,INCY,C,S) C C APPLIES THE GIVENS ROTATION C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER(I-N) C DIMENSION SX(1),SY(1) C IF (N.LE.0) RETURN C IF (INCX.EQ.1.AND.INCY.EQ.1) GO TO 20 C IX=1 IY=1 IF (INCX.LT.0) IX= (-N+1)*INCX+1 IF (INCY.LT.0) IY= (-N+1)*INCY+1 C DO 10 I=1,N STEMP = C*SX(IX) + S*SY(IY) SY(IY)= C*SY(IY) - S*SX(IX) SX(IX)= STEMP IX= IX + INCX IY= IY + INCY 10 CONTINUE RETURN C 20 DO 30 I=1,N STEMP = C*SX(I) + S*SY(I) SY(I) = C*SY(I) - S*SX(I) SX(I) = STEMP 30 CONTINUE RETURN END