SUBROUTINE LDP (G,MDG,M,N,H,X,XNORM,W,INDEX,MODE) LDP00100 C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1974 MAR 1 LDP00200 C TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974LDP00300 C LDP00400 C ********** LEAST DISTANCE PROGRAMMING ********** LDP00500 C LDP00600 INTEGER INDEX(M) LDP00700 DIMENSION G(MDG,N), H(M), X(N), W(1) LDP00800 ZERO=0. LDP00900 ONE=1. LDP01000 IF (N.LE.0) GO TO 120 LDP01100 DO 10 J=1,N LDP01200 10 X(J)=ZERO LDP01300 XNORM=ZERO LDP01400 IF (M.LE.0) GO TO 110 LDP01500 C LDP01600 C THE DECLARED DIMENSION OF W() MUST BE AT LEAST (N+1)*(M+2)+2*M. LDP01700 C LDP01800 C FIRST (N+1)*M LOCS OF W() = MATRIX E FOR PROBLEM NNLS. LDP01900 C NEXT N+1 LOCS OF W() = VECTOR F FOR PROBLEM NNLS. LDP02000 C NEXT N+1 LOCS OF W() = VECTOR Z FOR PROBLEM NNLS. LDP02100 C NEXT M LOCS OF W() = VECTOR Y FOR PROBLEM NNLS. LDP02200 C NEXT M LOCS OF W() = VECTOR WDUAL FOR PROBLEM NNLS. LDP02300 C COPY G**T INTO FIRST N ROWS AND M COLUMNS OF E. LDP02400 C COPY H**T INTO ROW N+1 OF E. LDP02500 C LDP02600 IW=0 LDP02700 DO 30 J=1,M LDP02800 DO 20 I=1,N LDP02900 IW=IW+1 LDP03000 20 W(IW)=G(J,I) LDP03100 IW=IW+1 LDP03200 30 W(IW)=H(J) LDP03300 IF=IW+1 LDP03400 C STORE N ZEROS FOLLOWED BY A ONE INTO F.LDP03500 DO 40 I=1,N LDP03600 IW=IW+1 LDP03700 40 W(IW)=ZERO LDP03800 W(IW+1)=ONE LDP03900 C LDP04000 NP1=N+1 LDP04100 IZ=IW+2 LDP04200 IY=IZ+NP1 LDP04300 IWDUAL=IY+M LDP04400 C LDP04500 CALL NNLS (W,NP1,NP1,M,W(IF),W(IY),RNORM,W(IWDUAL),W(IZ),INDEX, LDP04600 * MODE) LDP04700 C USE THE FOLLOWING RETURN IF UNSUCCESSFUL IN NNLS.LDP04800 IF (MODE.NE.1) RETURN LDP04900 IF (RNORM) 130,130,50 LDP05000 50 FAC=ONE LDP05100 IW=IY-1 LDP05200 DO 60 I=1,M LDP05300 IW=IW+1 LDP05400 C HERE WE ARE USING THE SOLUTION VECTOR Y.LDP05500 60 FAC=FAC-H(I)*W(IW) LDP05600 C LDP05700 IF (DIFF(ONE+FAC,ONE)) 130,130,70 LDP05800 70 FAC=ONE/FAC LDP05900 DO 90 J=1,N LDP06000 IW=IY-1 LDP06100 DO 80 I=1,M LDP06200 IW=IW+1 LDP06300 C HERE WE ARE USING THE SOLUTION VECTOR Y.LDP06400 80 X(J)=X(J)+G(I,J)*W(IW) LDP06500 90 X(J)=X(J)*FAC LDP06600 DO 100 J=1,N LDP06700 100 XNORM=XNORM+X(J)**2 LDP06800 XNORM=SQRT(XNORM) LDP06900 C SUCCESSFUL RETURN. LDP07000 110 MODE=1 LDP07100 RETURN LDP07200 C ERROR RETURN. N .LE. 0. LDP07300 120 MODE=2 LDP07400 RETURN LDP07500 C RETURNING WITH CONSTRAINTS NOT COMPATIBLE.LDP07600 130 MODE=4 LDP07700 RETURN LDP07800 END LDP07900