SUBROUTINE LSEQIEQ(E,N,NN,W,X,G,M,MM,H,C,M1,D,F,B,Y1,Y,GWORK, $WKAREA,WKLDP,INDEX,MODE,TAU,HH,IP,GG,KRANK,NCASES,XY,YY,INITL,LAST $,XSQR,YSQR,PRINT,DEBUG,IUNIT) C C LINEAR LEAST SQUARES WITH: C C LINEAR EQUALITY AND C LINEAR INEQUALITY CONSTRAINTS C C SOLVE \\EX-F\\ SUBJECT TO : GX >= H C CX = D C WITH SEQUENTIAL ROW BY ROW ACCUMULATION OF E, AND PSEUDORANK C DETERMINATION USING THE INCLUSION CRITERIA VARIABLE : TAU. C C VARIABLE(PROBLEM DIMENSION) [VARIABLE(ACTUAL DIMENSION)] C C E(N+2,N+1) [E(NN,NN)] CREATED BY NCASES CALLS TO LSEQIEQ. IT EQUAL C THE QR DECOMPOSITION OF [E:F] I.E. R D C 0 E C 0 0 C (NOTATION OF LAWSON AND HANSON,1973) C R IS N BY N (OR REDUCED BY THE NUMBER OF EQUALITY C CONSTRAINTS. C :WHEN LAST = .TRUE. THE COVARIANCE MATRIX C IS OUTPUT IF THE PROBLEM IS FULL RANK. C W(N) [W(NN)] :SET OF COEFFICIENTS FOR THE ITH CASE. C :WHEN LAST = .FALSE., THE ITH CASE COEFFICIENTS ARE C RETURNED, MODIFIED FOR THE EQUALITY CONSTRAINTS. C :WHEN LAST = .TRUE., DIAGONAL ELEMENTS OF THE INVERSE C CORRELATION MATRIX,IF THE PROBLEM IS FULL RANK. C B :CONSTANT FOR THE ITH CASE. I.E. W*X = B (ITH CASE). C :WHEN LAST = .FALSE., RETURNED, MODIFIED FOR THE C EQUALITY CONSTRAINTS. C :WHEN LAST = .TRUE., RESIDUAL SUM OF SQUARES. C X(N) [X(NN)] :FINAL SOLUTION VECTOR. C G(M,N) [G(MM,NN)] :MATRIX OF INEQUALITY CONSTRAINT C COEFFICIENTS. C :THE QR DECOMPOSITION OF THE CORRESPONDING C MATRIX REDUCED BY THE NUMBER OF EQUALITY CONSTRAIN C H(M) [H(MM)] :INEQUALITY CONSTRAINT CONSTANTS. C :REDUCED BY THE QR DECOMPOSITION OF C. C C(M1,N) [C(NN,MM)] :EQUALITY CONSTRAINT MATRIX. C :CORRESPONDING QR DECOMPOSED MATRIX. C D(M1) [D(NN)] :EQUALITY CONSTRAINT CONSTANTS. C :RETURNED UNMODIFIED. C F(MM) :FINAL CONSTANT VECTOR FOR THE QR DECOMPOSED PROBLEM C Y1(NN) INTERMEDIATE SOLUTION VECTOR. C Y(NN) INTERMEDIATE SOLUTION VECTOR. C GWORK([MM,NN]) WORKSPACE NEEDED TO COMPUTE THE MODIFIED INEQUALITY C CONSTRAINT MATRIX TO SEND TO LDP . C WKAREA() [WKAREA(2*NN)] WORKSPACE. C WKLDP[((NN+1)*(MM+2)+2*MM)] WORKSPACE FOR LDP. C INDEX[(MM)] WORKSPACE PASSED TO NNLS. C MODE :INTEGER PARAMETER INDICATING TERMINATION CONDITION. C MODE = 0 NORMAL EXECUTION. C 1 EQUALITY CONSTRAINT MATRIX IS SINGULAR. C 2 MATRIX R OF THE QR DECOMPOSITION OF THE DATA MATRIX IS C SINGULAR (PROBLEM IS NOT FULL RANK). C 3 BAD DIMENSION CONDITION IN LDP. C 4 MAXIMUM NUMBER OF ITERATIONS (3*M) EXCEEDED IN NNLS. C 5 INEQUALITY CONSTRAINTS ARE INCOMPATIBLE. C C TAU :RELATIVE TOLERANCE (DATA UNCERTAINTY). C :ABSOLUTE TOLERANCE PARAMETER FOR DETERMINATION C OF PSEUDORANK. SET EQUAL TO * C \\E\\. C HH(NM1),GG(NM1) [HH(NN),GG(NN)] WORKSPACE FOR HOLDING THE C HOUSEHOLDER PIVOT ELEMENTS DURING PSEUDORANK C DETERMINATION. C IP(NM1) [IP(NN)] PERMUTATON MATRIX FOR THE COLUMNS OF E. C COLUMN INTERCHANGE ARISES DURING THE DETERMINATION C OF PSEUDORANK. C KRANK :PSEUDORANK OF THE QR DECOMPOSED PROBLEM. C NCASES :NUMBER OF OBSERVATIONS. C INITL :SET = .TRUE. ON FIRST CALL. C MUST BE SET TO .FALSE. FOR ALL SUBSEQUENT CALLS. C LAST :SET = .TRUE. ON LAST CALL. C I.E.THIS INITIALIZES THE LEAST SQUARES PROBLEM. C XSQR [NN] :INITIALIZE TO ZERO. C :SUM OF THE SQUARES OF OBSERVATIONS OF THE C INDEPENDENT VARIABLES MODIFIED FOR THE C EQUALITY CONSTRAINTS. C YSQR :INITIALIZE TO ZERO. C :SUM OF THE SQUARES OF OBSERVATIONS OF THE C DEPENDENT VARIABLE MODIFIED FOR THE EQUALITY C CONSTRAINTS. C PRINT :LOGICAL VARIABLE TO ALLOW OR SUPPRESS PRINTED C OUTPUT FROM THE PACKAGE. C DEBUG :LOGICAL VARIABLE TO ALLOW OR SUPPRESS EXTENSIVE C INTERMEDIATE RESULT OUTPUT. C IUNIT :NUMBER OF THE LOGICAL UNIT TO WHICH OUTPUT WILL C BE DIRECTED. C C ****************************************************************** C C THE ROUTINE LOGIC IS DERIVED FROM: C C THE OUTLINE CHAPT 23 OF: C LAWSON,C.L. AND R.J.HANSON (1973) C SOLVING LEAST SQUARES PROBLEMS C PRENTICE HALL, INC. NEW JERSEY C C MARK S. GHIORSO VERSION 1.0 JANUARY 1981 C UNIVERSITY OF WASHINGTON C REVISION 2.0 JULY 1981 C LAWRENCE BERKELEY LAB C UNIVERSITY OF CALIFORNIA C REVISION 3.0 OCTOBER 1981 C UNIVERSITY OF WASHINGTON C REVISION 3.1 OCTOBER 1982 C UNIVERSITY OF WASHINGTON C (CONVERSION TO ANSI STANDARD FORTRAN IV) C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER(I-N) C DIMENSION E(NN,NN), W(NN), X(NN), G(MM,NN), H(MM), C(NN,NN), D(NN) $, F(MM), Y1(NN), Y(NN), XY(MM), GWORK(MM,NN), XSQR(NN), HH(NN), GG $(NN), IP(NN) DIMENSION WKAREA(1), WKLDP(1) LOGICAL INITL,LAST,PRINT,DEBUG CHARACTER*6 NAMES(1) C IF (LAST) GO TO 120 IF (.NOT.INITL) GO TO 70 NCASES = 0 MODE = 0 IF (DEBUG) WRITE (IUNIT,10) 10 FORMAT (////,1X,37H*****LSEQIEQ IS RUNNING IN DEBUG MODE,//) L = 0 IF (M1.EQ.0) GO TO 70 C C PERFORM THE QR DECOMPOSITION OF C C DO 20 I=1,M1 20 CALL H12 (1,I,I+1,N,C(I,1),NN,Y(I),C(I+1,1),NN,1,M1-I) C C C IS NOW M1 BY M1, LOWER TRIANGULAR C IF (DEBUG) CALL MATOUT (2,C,M1,NN,N,'C ','QR DECOMPOSED EQUAL $ITY CONSTRAINT MATRIX.',IUNIT) IF (DEBUG) CALL MATOUT (1,Y,1,1,M1,'Y ','VECTOR U OF THE ABOV $E QR DECOMPOSITION.',IUNIT) C C SOLVE THE SYSTEM C*Y1 = D C IF (DIFF(C(1,1),0.0).EQ.0.0) GO TO 1000 Y1(1) = D(1)/C(1,1) IF (M1.EQ.1) GO TO 50 DO 40 I=2,M1 JJ = I-1 Y1(I) = D(I) DO 30 J=1,JJ 30 Y1(I) = Y1(I)-C(I,J)*Y1(J) IF (DIFF(C(I,I),0.0).EQ.0.0) GO TO 1000 40 Y1(I) = Y1(I)/C(I,I) C C Y1(M1) NOW CONTAINS THE LINEAR CONSTRAINT SOLUTION VECTOR C IF (DEBUG) CALL MATOUT (1,Y1,1,1,M1,'Y1 ','LINEARLY CONSTRAINE $D SOLUTION VECTOR.',IUNIT) C C DECOMPOSE THE INEQUALITY CONSTRAINT MATRIX IN PLACE C 50 DO 60 I=1,M1 60 CALL H12 (2,I,I+1,N,C(I,1),NN,Y(I),G(1,1),MM,1,M) C C G NOW CONTAINS [G1(M,M1)G2(M,N-M1)] C IF (DEBUG) CALL MATOUT (2,G,M,MM,N,'G ','G1(M,M1)G2(M,N-M1)-M $ODIFIED INEQUAL CONSTR MATRIX.',IUNIT) C C ****************************************************************** C C ENTRY POINT FOR NCASES .NE.1 C 70 NCASES = NCASES+1 IF (M1.EQ.0) GO TO 100 C C MODIFY THE OBSERVATION VECTOR C DO 80 I=1,M1 80 CALL H12 (2,I,I+1,N,C(I,1),NN,Y(I),W(1),1,1,1) C IF (DEBUG) CALL MATOUT (1,W,1,1,N,'W ','MODIFIED INPUT OBSERV $ATION VECTOR (LC).',IUNIT) C C MODIFY THE CONSTANT VECTOR C DO 90 I=1,M1 90 B = B-W(I)*Y1(I) C C ACCUMULATE THE SUM OF SQUARES VECTOR C 100 DO 110 I=1,N IF (I.LE.M1) GO TO 110 XSQR(I) = XSQR(I)+W(I)**2 110 CONTINUE YSQR = YSQR+B**2 C C UPDATE THE QR DECOMPOSITION OF E:F C CALL SEQHT (E,N-M1,NN,W(M1+1),B,T,L) C C RETURN POINT TO RETRIEVE ANOTHER CASE C RETURN C C AT THIS POINT Q[E:F] = R(N,N) F1(N) N C 0 E 1 C 0 0 1 C C N 1 C C R IS UPPER TRIANGULAR R X = D C 0 E C C NOTE:THE TOLERANCE OF EACH VARIABLE IN THE EQUATION IS GIVEN BY C E(I,I), I=1,N-M1. C 120 NM1 = N-M1 C IF (DEBUG) CALL MATOUT (2,E,NM1+2,NN,NM1+1,'E ','[R/0:F1/E] M $ATRIX OF MOD OBS VECTORS (LC).',IUNIT) C C ZERO THE LOWER SUBDIAGONAL ELEMENTS OF R C DO 130 I=2,NM1 IJ = I-1 DO 130 J=1,IJ 130 E(I,J) = 0.0 IF (PRINT) WRITE (IUNIT,140) NM1 140 FORMAT (1H1,19HSUBROUTINE LSEQIEQ:,//,20X,41HCOMPUTED UNNORMALIZED $ TOLERANCES FOR THE ,I5,9H REDUCED ,10HVARIABLES:,//) IF (PRINT) WRITE (IUNIT,150) (E(I,I),I=1,NM1) 150 FORMAT (20X,5G15.6) TEMP1 = ABS(E(NM1+1,NM1+1)) IF (PRINT) WRITE (IUNIT,160) TEMP1 160 FORMAT (//,20X,40HNORM OF THE RESIDUALS FOR THE UNBOUNDED , $8HPROBLEM:,G15.6) IF (PRINT) WRITE (IUNIT,170) 170 FORMAT (//,20X,27HANOVA (UNBOUNDED VARIABLES),//,35X,7HSUM OF , $7HSQUARES,1X,14HDEG OF FREEDOM,5X,7HF (SIG),/) XSQR(N+1) = YSQR IDF = NCASES-NM1-1 IF (IDF.LE.0) GO TO 200 ENM1SQ = E(NM1+1,NM1+1)**2/YSQR B = E(NM1+1,NM1+1)**2/FLOAT(IDF) FTEST = FLOAT(IDF)*(1.0/ENM1SQ-1.0)/FLOAT(NM1) IF (FTEST.GT.0.0) CALL MDFD (FTEST,NM1,IDF,SIG,IER) TEMP1 = YSQR*(1.0-ENM1SQ) TEMP2 = YSQR*ENM1SQ IF (PRINT) WRITE (IUNIT,180) TEMP1,NM1,FTEST,SIG,TEMP2,IDF 180 FORMAT (20X,11HREGRESSION:,3X,G13.6,2X,I5,10X,G10.3,2H (,F6.3,1H), $/,20X,10HRESIDUALS:,4X,G13.6,2X,I5) RR = SQRT(1.0-ENM1SQ) TEMP1 = RR**2-FLOAT(NM1)*(1.0-RR**2)/FLOAT(IDF) IF (PRINT) WRITE (IUNIT,190) RR,TEMP1 190 FORMAT (/,20X,37HMULTIPLE CORRELATION COEFFICIENT (R):,F10.5,/,20X $,19HADJUSTED R SQUARED:,F10.5) C C PERFORM A SINGLE VALUE DECOMPOSITION OF THE UNBOUNDED PROBLEM C 200 DO 210 I=1,NM1 DO 210 J=1,NM1 210 GWORK(I,J) = E(I,J) DO 220 I=1,NM1 220 WKLDP(I) = E(I,NM1+1) WKLDP(NM1+1) = E(NM1+1,NM1+1) ISCALE = 1 C C A(,),MDA,M,N = [GWORK(NM1,NM1),MM,NM1,NM1] C MDATA = [NCASES] C B() = [WKLDP(1---NM1+1)] C SING() = [WKLDP(NM1+2---4*NM1+1)] C NAMES() = [WKLDP(4*NM1+2---5*NM1+1)] C ISCALE C D() = [WKLDP(5*NM1+2---6*NM1+1)] C C LAWSON AND HANSON SUBROUTINE *SVA* IS USED FOR THE ANALYSIS C C CALL SVA(A,MDA,M,N,MDATA,B,SING,NAMES,ISCALE,D) C NAMES(1) = ' ' CALL SVA (GWORK,MM,NM1,NM1,NCASES,WKLDP(1),WKLDP(NM1+2),NAMES, $ISCALE,WKLDP(5*NM1+2)) C C SET THE ABSOLUTE TOLERANCE PARAMETER C = RELATIVE TOLERANCE (DATA UNCERTAINTY)*\\E\\ C = TAU* TAU = TAU*ABS(WKLDP(NM1+2)) C C WE WILL NOW PERFORM A PSEUDORANK DETERMINATION USING THE METHOD C OUTLINED IN ALGORITHM HFTI FROM CHAPT 14 (PROBLEM LS WITH C POSSIBLY DEFICIENT PSEUDORANK) IN C C.L. LAWSON AND R.J. HANSON (1974) C SOLVING LEAST SQUARES PROBLEMS C PRENTICE HALL, INC., ENGLEWOOD CLIFFS, N.J. C C K = 0 FACTOR = .001 DO 300 J=1,NM1 IF (J.EQ.1) GO TO 240 C C UPDATE SQUARED COLUMN LENGTHS AND FIND LMAX C LMAX = J DO 230 L=J,NM1 HH(L) = HH(L)-E(J-1,L)**2 IF (HH(L).GT.HH(LMAX)) LMAX = L 230 CONTINUE IF (DIFF(HMAX+FACTOR*HH(LMAX),HMAX)) 240,240,270 C C COMPUTE SQUARED COLUMN LENGTHS AND FIND LMAX C 240 LMAX = J DO 260 L=J,NM1 HH(L) = 0.0 DO 250 I=J,NM1 250 HH(L) = HH(L)+E(I,L)**2 IF (HH(L).GT.HH(LMAX)) LMAX = L 260 CONTINUE HMAX = HH(LMAX) C C LMAX HAS BEEN DETERMINED C DO COLUMN INTERCAHNGES IF NEEDED C 270 CONTINUE IP(J) = LMAX IF (IP(J).EQ.J) GO TO 290 DO 280 I=1,NM1 TMP = E(I,J) E(I,J) = E(I,LMAX) 280 E(I,LMAX) = TMP HH(LMAX) = HH(J) C C COMPUTE THE JTH TRANSFORMATION C 290 CALL H12 (1,J,J+1,NM1,E(1,J),1,HH(J),E(1,J+1),1,NN,NM1-J) 300 CALL H12 (2,J,J+1,NM1,E(1,J),1,HH(J),E(1,NM1+1),1,1,1) C IF (.NOT.DEBUG) GO TO 320 NM3 = 1 310 NM4 = MIN0(NM3+8,NM1) WRITE (IUNIT,330) (J,IP(J),J=NM3,NM4) NM3 = NM3+9 IF (NM3.GT.NM1) GO TO 320 GO TO 310 C 320 CONTINUE 330 FORMAT (/,1X,15(1H*),/,1X,33H*DEBUG OUTPUT*PERMUTATION MATRIX, $29H FOR COLUMN INTERCHANGES OF E,/,1X,15(1H*),//,14X,9(1X,3HIP(,I3 $,3H)= ,I3)) IF (DEBUG) CALL MATOUT (1,HH,1,1,NM1,'HH ','HOUSEHOLDER PIVOT $ELEMENTS FROM COLUMN INTER (LC).',IUNIT) IF (DEBUG) CALL MATOUT (2,E,NM1+2,NN,NM1+1,'E ','[R/0:F1/E] M $ATRIX AFTER COLUMN INTER (PR).',IUNIT) C C DETERMINE THE PSEUDORANK, KRANK, USING THE TOLERANCE, TAU C DO 340 J=1,NM1 IF (ABS(E(J,J)).LE.TAU) GO TO 350 340 CONTINUE KRANK = NM1 GO TO 360 C 350 KRANK = J-1 360 KP1 = KRANK+1 C C COMPUTE THE NORMS OF THE RESIDUAL VECTORS C TMP = 0.0 IF (KP1.GT.NM1) GO TO 380 DO 370 I=KP1,NM1 370 TMP = TMP+E(I,NM1+1)**2 380 RNORM = SQRT(TMP) C C SPECIAL FOR PSEUDORANK EQUALS ZERO C IF (KRANK.GT.0) GO TO 400 DO 390 I=1,NM1 390 X(I) = 0.0 GO TO 900 C C IF THE PSEUDORANK IS LESS THAN NM1, COMPUTE HOUSEHOLDER C DECOMPOSITION OF THE FIRST KRANK ROWS C 400 IF (KRANK.EQ.NM1) GO TO 420 DO 410 II=1,KRANK I = KP1-II 410 CALL H12 (1,I,KP1,NM1,E(I,1),NN,GG(I),E,NN,1,I-1) C IF (DEBUG) CALL MATOUT (1,GG,1,1,KRANK,'GG ','HOUSEHOLDER PIVO $T ELEMENTS (PR PROBLEM).',IUNIT) IF (DEBUG) CALL MATOUT (2,E,KRANK,NN,KRANK,'E ','[R] MATRIX ( $RANK DEFICIENT PROBLEM).',IUNIT) C 420 CONTINUE C C REORDER THE LINEARLY MODIFIED INEQUALITY CONSTRAINT MATRIX, G, C ACCORDING TO THE COLUMN INTERCHANGES USED IN CONSTRUCTING RK(T) C DO 440 JJ=1,NM1 J = JJ IF (IP(J).EQ.J) GO TO 440 L = IP(J) LM1 = L+M1 JM1 = J+M1 DO 430 I=1,M TMP = G(I,LM1) G(I,LM1) = G(I,JM1) 430 G(I,JM1) = TMP 440 CONTINUE C IF (DEBUG) CALL MATOUT (2,G,M,MM,N,'G ','REORD (CI/PR) MODIFI $ED (LC) INEQ CON MATRIX.',IUNIT) C C IF THE PSEUDORANK IS LESS THAN NM1, APPLY HOUSEHOLDER DECOMPOSITIO C TO THE FIRST KRANK ROWS OF THE INEQUALITY CONTRAINT MATRIX G C C IF (KRANK.EQ.NM1) GO TO 460 DO 450 II=1,KRANK I = KP1-II 450 CALL H12 (2,I,KP1,NM1,E(I,1),NN,GG(I),G,MM,1,M) C IF (DEBUG) CALL MATOUT (2,G,M,MM,N,'G ','MOD G (AS ABOVE) DEC $OMP FOR PSEUDORANK.',IUNIT) C C OUTPUT INFO ON THE PSEUDORANK SECTION C 460 IF (PRINT) WRITE (IUNIT,470) KRANK,RNORM,TAU 470 FORMAT (//,20X,11HPSEUDORANK:,I5,5X,14HRESIDUAL NORM:,G12.5,5X, $29HABSOLUTE TOLERANCE PARAMETER:,G10.3,//,20X,36HREORGANIZED UNNOR $MALIZED TOLERANCES:,//) IF (PRINT) WRITE (IUNIT,150) (E(I,I),I=1,KRANK) C COMPUTE THE QUANTITY KR(-1) AND STORE IT IN GWORK C C COPY E(KRANK,KRANK) INTO GWORK C DO 480 I=1,KRANK DO 480 J=1,KRANK 480 GWORK(I,J) = E(I,J) C C ZERO THE LOWER DIAGONAL OF GWORK C DO 490 I=1,KRANK IF (KRANK.EQ.I) GO TO 500 II = 1+I DO 490 J=II,KRANK 490 GWORK(J,I) = 0.0 500 CONTINUE C C INVERT GWORK IN PLACE = KR(-1) C D1 = -1. D2 = 0.0 CALL LINV3F (GWORK,DUMMY,1,KRANK,MM,D1,D2,WKAREA,IER) IF (IER.EQ.130) GO TO 1010 C IF (DEBUG) CALL MATOUT (2,GWORK,KRANK,MM,KRANK,'GWORK ','THE MATR $IX K*R(-1).',IUNIT) C C FORM G2*K*R(-1) IN GWORK IN PLACE C DO 530 I=1,M DO 510 J=1,KRANK JM1 = J+M1 510 XY(J) = G(I,JM1) DO 520 J=1,KRANK JM1 = J+M1 G(I,JM1) = 0.0 DO 520 K=1,KRANK 520 G(I,JM1) = G(I,JM1)+XY(K)*GWORK(K,J) 530 CONTINUE C C COPY G2*K*R(-1) BACK INTO GWORK C DO 540 I=1,M DO 540 J=1,KRANK JM1 = J+M1 540 GWORK(I,J) = G(I,JM1) C IF (DEBUG) CALL MATOUT (2,GWORK,M,MM,KRANK,'GWORK ','THE MATRIX G $2*K*R(-1) (THE LDP MATRIX).',IUNIT) C C MODIFY THE INEQUALITY CONSTRAINT VECTOR TO ACCOUNT FOR THE C EQUALITY CONSTRAINTS : H-G1*Y1 - G2*K*R(-1)*F1 C DO 570 I=1,M F(I) = H(I) IF (M1.EQ.0) GO TO 560 DO 550 J=1,M1 550 F(I) = F(I)-G(I,J)*Y1(J) 560 DO 570 J=1,KRANK 570 F(I) = F(I)-GWORK(I,J)*E(J,NM1+1) C IF (DEBUG) CALL MATOUT (1,F,1,1,M,'F ','H-G1*Y1-G2*K*R(-1)*F1 $ OF THE LDP PROBLEM.',IUNIT) C C NOTE:F1 WAS QR DECOMPOSED AS THE E(I,NM1+1) ELEMENT C C WE NOW HAVE THE LDP PROBLEM MINIMIZE \\Z\\ SUBJECT TO GWORK*Z>=F C CALL LDP (GWORK,MM,M,KRANK,F,X,XNORM,WKLDP,INDEX,MODE) C C UPON RETURN X CONTAINS THE LDP SOLUTION C C CHECK ERROR CONDITIONS C IF (PRINT) WRITE (IUNIT,580) XNORM 580 FORMAT (///,20X,33HVALUE OF XNORM RETURNED FROM LDP:,G15.6) IF (PRINT) WRITE (IUNIT,590) 590 FORMAT (/,20X,42HCOMPUTED DUAL VECTOR FOR THE NNLS PROBLEM:,//) J = (KRANK+1)*(M+2)+M+1 L = J+M-1 IF (PRINT) WRITE (IUNIT,150) (WKLDP(I),I=J,L) GO TO (600,1020,1030,1040), MODE 600 MODE = 0 C IF (DEBUG) CALL MATOUT (1,X,1,1,KRANK,'X ','SOLUTION VECTOR F $ROM THE LDP PROBLEM.',IUNIT) C C NOW WE MUST BACK TRANSFORM THE VARIABLES C C FORM X = Z + F1 [= ITH ELEMENT E(I,NM1+1)] C 610 DO 620 I=1,KRANK 620 X(I) = X(I)+E(I,NM1+1) C C SOLVE FOR Y3 : RK(T)*Y3 = X RK(T) IS UPPER TRIANGULAR C DO 650 II=1,KRANK SM = 0.D0 I = KRANK+1-II IF (I.EQ.KRANK) GO TO 640 JJ = I+1 DO 630 J=JJ,KRANK 630 SM = SM+E(I,J)*X(J) 640 SM1 = SM 650 X(I) = (X(I)-SM1)/E(I,I) C IF (DEBUG) CALL MATOUT (1,X,1,1,KRANK,'X ','VECTOR Y3:R*K(T)* $Y3 = SOLUTION FROM LDP.',IUNIT) C C COMPLETE COMPUTATION OF THE SOLUTION VECTOR C IF (KRANK.EQ.NM1) GO TO 680 DO 660 J=KP1,NM1 660 X(J) = 0.0 DO 670 I=1,KRANK 670 CALL H12 (2,I,KP1,NM1,E(I,1),NN,GG(I),X,1,1,1) C IF (DEBUG) CALL MATOUT (1,X,1,1,NM1,'X ','VECTOR Y3 EXPANDED $TO FULL RANK.',IUNIT) C C REORDER THE SOLUTION VECTOR TO COMPENSATE FOR THE COLUMN C INTERCHANGES C 680 DO 690 JJ=1,NM1 J = NM1+1-JJ IF (IP(J).EQ.J) GO TO 690 L = IP(J) TMP = X(L) X(L) = X(J) X(J) = TMP 690 CONTINUE C IF (DEBUG) CALL MATOUT (1,X,1,1,NM1,'X ','VECTOR Y3 REORD FOR $ COLUMN INTER.',IUNIT) C C X IS THE SOLUTION IF THERE ARE NO EQUALITY CONSTRAINTS C C REPLACE THE MATRIX R OF THE QR DECOMPOSITION OF E (STORED AS C E(NM1,NM1) WITH THE INVERSE CORRELATION MATRIX. EXECUTE C ALGORITHM COV (12.12) OF: C LAWSON,C.L. AND R.J. HANSON C "SOLVING LEAST SQUARES PROBLEMS" C PRENTICE HALL, 1974 P.281 C IF (KRANK.LT.NM1.OR.NCASES.LE.(NM1+1)) GO TO 900 C C PERMUTE THE SUM OF SQUARES VECTOR INTO XY C DO 700 I=1,NM1 IM1 = I+M1 700 XY(I) = SQRT(XSQR(IM1)/FLOAT(NCASES-1)) DO 710 I=1,NM1 IF (IP(I).EQ.I) GO TO 710 L = IP(I) TMP = XY(L) XY(L) = XY(I) XY(I) = TMP 710 CONTINUE DO 720 I=1,NM1 DO 720 J=1,I 720 E(J,I) = E(J,I)/XY(I) NM2 = NM1+1 XSQRT = SQRT(YSQR/FLOAT(NCASES-1)) DO 730 I=1,NM2 730 E(I,NM2) = E(I,NM2)/XSQRT DO 740 J=1,NM2 740 E(J,J) = 1.0/E(J,J) IF (NM2.EQ.1) GO TO 770 DO 760 I=1,NM1 IP1 = I+1 DO 760 J=IP1,NM2 JM1 = J-1 SM = 0.D0 DO 750 L=I,JM1 750 SM = SM+E(I,L)*E(L,J) 760 E(I,J) = -SM*E(J,J) 770 DO 790 I=1,NM2 DO 790 J=I,NM2 SM = 0.D0 DO 780 L=J,NM2 780 SM = SM+E(I,L)*E(J,L) 790 E(I,J) = SM IF (NM2.LT.2) GO TO 860 DO 850 II=2,NM2 I = NM2+1-II IF (IP(I).EQ.I) GO TO 850 K = IP(I) TMP = E(I,I) E(I,I) = E(K,K) E(K,K) = TMP IF (I.EQ.1) GO TO 810 DO 800 L=2,I TMP = E(L-1,I) E(L-1,I) = E(L-1,K) 800 E(L-1,K) = TMP 810 IP1 = I+1 KM1 = K-1 IF (IP1.GT.KM1) GO TO 830 DO 820 L=IP1,KM1 TMP = E(I,L) E(I,L) = E(L,K) 820 E(L,K) = TMP 830 IF (K.EQ.NM2) GO TO 850 KP1 = K+1 DO 840 L=KP1,NM2 TMP = E(I,L) E(I,L) = E(K,L) 840 E(K,L) = TMP 850 CONTINUE 860 CONTINUE C IF (DEBUG) CALL MATOUT (2,E,NM1,NN,NM1,'E ','INV CORREL MATRI $X [A(T)*A](-1).',IUNIT) IF (M1.EQ.0) GO TO 880 DO 870 I=1,M1 870 W(I) = 0.0 880 DO 890 I=1,NM1 IM1 = I+M1 890 W(IM1) = E(I,I) C 900 IF (M1.EQ.0) GO TO 940 C C PUT X INTO THE LOWER PART OF THE Y1 VECTOR TO FORM Y1:X C DO 910 I=1,NM1 IM1 = I+M1 910 Y1(IM1) = X(I) DO 920 I=1,N 920 X(I) = Y1(I) C IF (DEBUG) CALL MATOUT (1,X,1,1,N,'X ','SOLUTION VECTOR PRIOR $ TO LC REORGANIZATION.',IUNIT) C C NOW BACK TRANSFORM Y1 INTO X USING THE QR DECOMPOSITION OF C C DO 930 II=1,M1 I = M1+1-II 930 CALL H12 (2,I,I+1,N,C(I,1),NN,Y(I),X(1),1,1,1) C C X NOW CONTAINS THE SOLUTION VECTOR TO THE PROBLEM C 940 CONTINUE IF (.NOT.PRINT) RETURN WRITE (IUNIT,950) 950 FORMAT (//,5X,29HLSEQIEQ - X = SOLUTION VECTOR,/,15X,34HSE = STAND $ARD ERROR F = PARTIAL F,/) DO 980 I=1,N IF (NM1.GT.KRANK.OR.M1.GT.0.OR.(NM1+1).GT.NCASES) GO TO 970 TEMP1 = 0.0 TEMP2 = 0.0 IF (XSQR(I).NE.0.0) TEMP1 = SQRT(B*W(I)/XSQR(I)) IF (W(I).NE.0.0) TEMP2 = XSQR(I)*X(I)**2/(B*W(I)) WRITE (IUNIT,960) I,X(I),TEMP1,TEMP2 960 FORMAT (5X,2HX(,I3,4H) = ,G25.16,5X,5HSE = ,G25.16,5X,4HF = ,G25. $16) GO TO 980 C 970 WRITE (IUNIT,960) I,X(I) 980 CONTINUE C C OUTPUT THE UNBOUNDED INVERSE CORRELATION MATRIX DIAGONAL ELEMENTS C IF ((NM1+1).GT.NCASES) RETURN IF ((KRANK.EQ.NM1).AND.PRINT) WRITE (IUNIT,990) 990 FORMAT (//,20X,42HDIAGONAL ELEMENTS OF THE UNBOUNDED INVERSE, $20H CORRELATION MATRIX:,//) IF ((KRANK.EQ.NM1).AND.PRINT) WRITE (IUNIT,150) (E(I,I),I=1,NM1) RETURN C C PROCESS ERROR REGIONS C C ****************************************************************** C C EQUALITY CONSTRAINT MATRIX IS SINGULAR C 1000 MODE = 1 RETURN C C MATRIX R OF THE QR DECOMPOSITION OF THE DATA MATRIX IS SINGULAR C 1010 MODE = 2 RETURN C C BAD DIMENSION CONDITION IN LDP C 1020 MODE = 3 RETURN C C MAXIMUM NUMBER (3*M) OF ITERATIONS EXCEEDED IN NNLS C 1030 MODE = 4 GO TO 610 C C INEQUALITY CONSTRAINTS ARE INCOMPATIBLE C 1040 MODE = 5 RETURN END