sub bupd.for

非线性回归问题SQP解法

	SUBROUTINE BUPD(N,ME,MI,M,NAMAX,X,NEWX,OGRA,CGRA,
     *			U,B,S,Y,BS,LG,LGNEW,W3,W4,W5,GRAD,IT,NG)
      DOUBLE PRECISION X(N),NEWX(N),OGRA(N),CGRA(NAMAX,M),U(M),
	*        B(NAMAX,M),LG(N),LGNEW(N),S(N),Y(N),BS(N),W3(N),W4(N),
     *        W5(N),WMAX,W,SY,SBS,THETA,THETA0,EPS
	LOGICAL GRAD
	DATA EPS/1.0D-5/
	DO 1 J=1,N
		LG(J)=OGRA(J)
		DO 1 I=1,M
			LG(J)=LG(J)-U(I)*CGRA(J,I)
1	CONTINUE
	IF (GRAD)THEN
		CALL GVAL(N,ME,MI,M,NAMAX,NEWX,OGRA,CGRA)
		NG=NG+1
	ELSE
		CALL DIF(N,ME,MI,M,NAMAX,NEWX,OGRA,CGRA,W3,W4,W5,NF)
	ENDIF
	NN=5*N
	IF (MOD(IT,NN).EQ.0)THEN
		DO 2 I=1,N
			DO 2 J=1,N
				IF (I.EQ.J)THEN
					B(I,J)=1.0D0
				ELSE
					B(I,J)=0
				ENDIF
2		CONTINUE
		GOTO 10
	ENDIF
	DO 3 J=1,N
		LGNEW(J)=OGRA(J)
		DO 3 I=1,M
			LGNEW(J)=LGNEW(J)-U(I)*CGRA(J,I)
3	CONTINUE
	DO 4 J=1,N
		S(J)=NEWX(J)-X(J)
		Y(J)=LGNEW(J)-LG(J)
4	CONTINUE
	DO 5 I=1,N
		BS(I)=0.0
		DO 5 J=1,N
			BS(I)=BS(I)+B(I,J)*S(J)
5	CONTINUE
	WMAX=0.0
	DO 6 J=1,N
		W=Y(J)-BS(J)
		IF (DABS(W).GT.WMAX) WMAX=DABS(W)
6	CONTINUE
	IF (WMAX.LT.EPS)GOTO 10
	SY=0.0
	SBS=0.0
	DO 7 J=1,N
		SY=SY+S(J)*Y(J)
		SBS=SBS+S(J)*BS(J)
7	CONTINUE
	IF (DSQRT(SBS).LT.EPS)GOTO 10
	IF (SY.LT.0.2D0*SBS)THEN
		THETA=0.8D0*SBS/(SBS-SY)
		THETA0=1.0D0-THETA
		SY=0.0
		DO 8 J=1,N
			Y(J)=THETA*Y(J)+THETA0*BS(J)
			SY=SY+S(J)*Y(J)
8		CONTINUE
	ENDIF
	C1=1.0D0/SY
	C2=1.0D0/SBS
	DO 9 I=1,N
		DO 9 J=1,I
			B(I,J)=B(I,J)+C1*Y(I)*Y(J)-C2*BS(I)*BS(J)
			B(J,I)=B(I,J)
9	CONTINUE
10	DO 11 J=1,N
		X(J)=NEWX(J)		
11	CONTINUE
	END