sqp.for

非线性回归问题SQP解法

*********************************************************************
* THIS PROGRAM IS FOR SOLVING THE NONLINEAR PROGRAM PROBLEM         *
*	MINMIZE		F(X)											  *
*	SUBJECT TO	CI(X).EQ.0,I=1,2...,ME							  *
*				CI(X).GE.0,I=ME+1,...,M							  *
* where X is a N-vector,F(X) and Ci(X) are twicecontinupusly		  *
* differentiable funtions											  *
*********************************************************************
	PARAMETER(NMAX=5,MEMAX=5,MIMAX=10)
	PARAMETER(MMAX=MEMAX+MIMAX)
	PARAMETER(NAMAX=NMAX+2*MEMAX+MIMAX,MAMAX=3*MEMAX+2*MIMAX)
	DOUBLE PRECISION X(NMAX),U(MAMAX),OFUN,OGRA(NAMAX),CFUN(MMAX),
     *		CGRA(NAMAX,MAMAX),B(NAMAX,NAMAX),D1(NAMAX),D2(NMAX),
     *		PK(NMAX),AK(NAMAX,MMAX),PFUN,PPRA,QPOBJ,RHS(MAMAX),
     *		CGRA1(NAMAX,MMAX),UMAX,NEWX(NMAX),OLDPF,U2(MMAX),S(NMAX),
     *		Y(NMAX),W0(NAMAX),W1(NAMAX),W2(NAMAX),W3(NAMAX),
     *		W4(NAMAX),W5(NAMAX),WL(NAMAX,NAMAX),WQ(NAMAX,NAMAX),
     *		WR(NAMAX,NAMAX),NM(NAMAX,NAMAX),SCL(MAMAX),
     *		INITPP,EPS1,EPS2,EPS3,BOUND,RHO
	INTEGER QPINDX,QPIT,JSET(MAMAX)
	LOGICAL VALID,GRAD,FEASBL,CORR,BASIC,MJ(MAMAX),SCALE
	DATA EPS1,EPS2,EPS3,BOUND/1.0D-3,1.0D-3,1.0D-6,1.0D+20/
	DATA INITPP,RHO/1.0D0,1.0D+1/
	DATA GRAD,SCALE/.TRUE.,.TRUE./
*----------Read in problem and initial data from file.
	CALL INPUT(N,ME,MI,M,NMAX,MEMAX,MIMAX,X,VALID)
	NA=N+2*ME+MI
	MA=3*ME+MI
	MEA=ME
	IF(.NOT.VALID)THEN
		WRITE(6,100)
100		FORMAT('*INPUT FILE INVALID')
		STOP
	ENDIF
*-----------Initialize.
	PPRA=INITPP
	ITMAX=MIN0(200,MAX0(100,20*(N+M)))
	IT=0
	NQP=0
	ITQP=0
	NF=0
	NG=0
	DO 1 I=1,N
		DO 1 J=1,N
			IF(I.EQ.J)THEN
				B(I,J)=1.0D0
			ELSE
				B(I,J)=0.0
			ENDIF
1	CONTINUE
	CALL FVAL(N,ME,MI,M,X,OFUN,CFUN)
	NF=NF+1
	IF (GRAD)THEN
		CALL GVAL(N,ME,MI,M,NAMAX,X,OGRA,CGRA)
		NG=NG+1
	ELSE
		CALL DIF(N,ME,MI,M,NAMAX,X,OGRA,CGRA,W2,W3,W4,NF)
	ENDIF
	CALL PVAL(ME,MI,M,OFUN,CFUN,PPRA,PFUN)
*-----------Iteration.
*-----------Print out intermediate results.
2	CALL PRNT(N,ME,MI,M,X,U,OFUN,CFUN,PPRA,PFUN,IT)
*-----------Check termination criteria.
	CALL CONTES(N,ME,M,NAMAX,X,U,PFUN,OLDPF,CFUN,
     *			OGRA,CGRA,W0,BOUND,EPS1,EPS2,INDEX)
	IF (INDEX.LE.1.OR.INDEX.EQ.4)GOTO 8
	IT=IT+1
*-----------Solve basic QP subproblem.
	DO 3 I=1,M
		RHS(I)=-CFUN(I)
3	CONTINUE
	CALL QPSOLV(N,ME,M,NAMAX,CGRA,RHS,OGRA,B,D1,U,WL,WQ,W0,W1,W2,W3,
     *		W4,WR,NM,JSET,MJ,QPINDX,QPIT,SCALE,MAMAX,SCL)
	NQP=NQP+1
	ITQP=ITQP+QPIT
	IF (QPINDX.EQ.1.OR.QPINDX.EQ.2.OR.QPINDX.EQ.4)THEN
		INDEX=6
		GOTO 8
	ENDIF
	IF (QPINDX.EQ.3)THEN
*-----------Basic QP subproblem is infeasible.
		FEASBL=.FALSE.
		CALL AUXQP(N,ME,MI,M,NA,MEA,MA,NAMAX,CGRA,RHS,OGRA,B,PPRA)
		CALL QPSOLV(NA,MEA,MA,NAMAX,CGRA,RHS,OGRA,B,D1,U,WL,WQ,W0,W1,
     *			W2,W3,W4,WR,NM,JSET,MJ,QPINDX,QPIT,SCALE,MAMAX,SCL)
		NQP=NQP+1
		ITQP=ITQP+QPIT
		IF (QPINDX.GT.0)THEN
			INDEX=6
			GOTO 8
		ENDIF
		BASIC=.TRUE.
		GOTO 5
	ELSE
	FEASBL=.TRUE.
	ENDIF
*-----------Check if second-order correction is necessary.
	CALL CHECK(N,ME,MI,M,NAMAX,X,CFUN,CGRA,CGRA1,
     *			D1,AK,W0,W1,W2,W3,W4,W5,GRAD,CORR,NF,NG)
	IF (.NOT.CORR)THEN
		BASIC=.TRUE.
		GOTO 5
	ENDIF
*-----------Solve modified QP subproblem.
	DO 4 J=1,N
		PK(J)=OGRA(J)
		DO 4 I=1,M
			PK(J)=PK(J)+5.0D-1*U(I)*(CGRA1(J,I)-CGRA(J,I))
4	CONTINUE
	CALL QPSOLV(N,ME,M,NAMAX,AK,RHS,PK,B,D2,U2,WL,WQ,W0,W1,W2,W3,
     *			W4,WR,NM,JSET,MJ,QPINDX,QPIT,SCALE,MAMAX,SCL)
	NQP=NQP+1
	ITQP=ITQP+QPIT
	IF (QPINDX.EQ.0)THEN
		BASIC=.FALSE.
	ELSE
		BASIC=.TRUE.
	ENDIF
*-----------Update penalty parameter.
5	IF (.NOT.FEASBL)GOTO 7
	UMAX=0.0
	DO 6 I=1,M
		UMAX=DMAX1(DABS(U(I)),UMAX)
6	CONTINUE
	IF (UMAX.GT.PPRA)THEN
		PPRA=UMAX+RHO
		CALL PVAL(ME,MI,M,OFUN,CFUN,PPRA,PFUN)
	ENDIF
*-----------Perform one-dimensional search.
7	CALL QPVAL(N,NA,NAMAX,OGRA,B,D1,PPRA,QPOBJ,FEASBL)
	CALL LSEARC(N,ME,MI,M,X,NEWX,D1,D2,BASIC,OFUN,CFUN,
     *			QPOBJ,PPRA,PFUN,OLDPF,EPS3,INDEX,NF)
*-----------Update current solution and approximate Hessian matrix.
	CALL BUPD(N,ME,MI,M,NAMAX,X,NEWX,OGRA,CGRA,U,B,S,Y,
     *			W0,W1,W2,W3,W4,W5,GRAD,IT,NG)
      IF (INDEX.EQ.5.OR.IT.GT.ITMAX)GOTO 8
	GOTO 2
*-----------Print out computatonal results.
8	CALL PVAL(ME,MI,M,OFUN,CFUN,PPRA,PFUN)
	CALL OUTPUT(N,ME,MI,M,NAMAX,X,U,OFUN,CFUN,OGRA,CGRA,W0,
     *			PPRA,PFUN,INDEX,QPINDX,IT,NQP,NF,NG,ITQP)
	STOP
	END