sub qpsolv.for

非线性回归问题SQP解法

      SUBROUTINE QPSOLV(N,ME,M,NAMAX,A,B,C,G,X,V,L,Q,W,AS,DV,ZV,RV,
     *		R,NM,JSET,MJ,QPINDX,IT,SCALE,MAMAX,SCL)
	DOUBLE PRECISION A(NAMAX,M),B(M),C(N),G(NAMAX,N),X(N),V(M),
	*        R(NAMAX,N),NM(NAMAX,N),L(NAMAX,N),W(N),Q(NAMAX,N),AS(N),
     *		DV(N),ZV(N),RV(M),US,ZERO,SCL(MAMAX),OSCL,ZMAX,T1,T2,T		
      INTEGER JSET(M),QPINDX
	LOGICAL MJ(M),SCALE,TRANS
	DATA ZERO/1.0D-8/
	IF (SCALE)THEN
		CALL SCAL(N,M,NAMAX,MAMAX,G,C,A,B,OSCL,SCL,ZERO)
	ENDIF
*----------Find the initial solution.
	CALL INITIA(N,ME,M,NAMAX,A,B,C,G,X,V,MJ,JSET,NJ,R,NM,L,W,Q,
     *			QPINDX,ZERO)
	IF (QPINDX.NE.0) RETURN
	IT=0
	ITMAX=10*(M+N)
*----------Check optimality.
1	CALL CHKOPT(N,ME,M,NAMAX,A,B,X,MJ,JS,US,OPTM,ZERO)
	IF (JS.EQ.0) THEN
		QPINDX=0
		IF (SCALE)THEN
			CALL QPOUT(N,M,NAMAX,MAMAX,G,C,A,B,V,OSCL,SCL)
		ENDIF
		RETURN
	ENDIF
	DO 2 J=1,N
		AS(J)=A(J,JS)
2	CONTINUE
*----------Update the current solution.
3	CONTINUE
	IT=IT+1
	IF (IT.GT.ITMAX) THEN
		QPINDX=4
		WRITE(6,100)ITMAX
100		FORMAT(//'NUMBER OF ITERATION EXCEEDS ITERATION LIMIT',I5)
		RETURN
	ENDIF
*----------Compute the vector DV
	DO 4 I=1,N
		DV(I)=0.0
		DO 4 J=1,N
			DV(I)=DV(I)+NM(J,I)*AS(J)
4	CONTINUE
*----------Compute the vectors ZV and RV.
	DO 5 J=1,N
		ZV(J)=0.0
		DO 5 I=NJ+1,N
			ZV(J)=ZV(J)+NM(J,I)*DV(I)
5	CONTINUE
	DO 6 I=1,NJ
		RV(I)=DV(I)
6	CONTINUE
	TRANS=.FALSE.
	CALL BWDSUB(NJ,NMAX,R,TRANS,RV)
*----------Check if (a)ZV<>0 or (b)ZV=0.
	ZMAX=0.0
	DO 7 J=1,N
		ZMAX=DMAX1(DABS(ZV(J)),ZMAX)
7	CONTINUE
	IF (ZMAX.LT.ZERO)GOTO 15
*----------Case(a):ZV<>0.
*----------Compute T1.
	CALL RTEST(M,ME,NJ,V,RV,JSET,T1,JK,JP,ZERO)
*----------Compute T2
	TS=0.0
	DO 8 J=1,N
		TS=TS+AS(J)*ZV(J)
8	CONTINUE
	T2=-US/TS
	IF (T1.LT.T2)GOTO 11
*----------Case(a-1):T1>=T2.
*----------Constrain JS is added to JSET(full step)
	DO 9 J=1,N
		X(J)=X(J)+t2*ZV(J)
9	CONTINUE
	DO 10 I=1,NJ
		V(JSET(I))=V(JSET(I))-T2*RV(I)
10	CONTINUE
	V(JS)=V(JS)+T2
	WRITE(6,110)JS
110	FORMAT('(ADD CONSTRAINT ',I3,':FULL STEP)')
	CALL ADD(N,NMAX,M,NJ,MJ,JSET,R,NM,DV,JS,ZERO)
	GOTO 1
*----------Case(a-2):T1<t2.
*----------Constraint JK is dropped from JSET (partial step).
11	DO 12 J=1,N
		X(J)=X(J)+T1*ZV(J)
12	CONTINUE
	DO 13 I=1,NJ
		V(JSET(I))=V(JSET(I))-T1*RV(I)
13	CONTINUE
	V(JS)=V(JS)+T1
	US=-B(JS)
	DO 14 J=1,N
		US=US+AS(J)*X(J)
14	CONTINUE
	WRITE(6,120)JK
120	FORMAT('(DROP CONSTRAIN',I3,':PARTIAL STEP)')
	CALL DROP(N,NMAX,M,NJ,MJ,JSET,R,NM,JK,JP,ZERO)
	GOTO 3
*----------Case(b):ZV=0.
*----------Check if RV(I)<=0 for active inequality contraints.
15	CALL RTEST(M,ME,NJ,V,RV,JSET,T,JK,JP,ZERO)
	IF (JK.NE.0)GOTO 16
*----------Case(b-1).
*----------Problem has turned out to be infeasible.
	QPINDX=3
	WRITE(6,130)
130	FORMAT(/'*PROBLEM IS INFEASTIBLE*')
	RETURN
*----------Case(b-2).
*----------Contraint JK is dropped from JSET (partial step).
16	DO 17 I=1,NJ
		V(JSET(I))=V(JSET(I))-T*RV(I)
17	CONTINUE
	V(JS)=T
	WRITE(6,120)JK
	CALL DROP(N,NMAX,M,NJ,MJ,JSET,R,NM,JK,JP,ZERO)
	GOTO 3
	END