sub initia.for

非线性回归问题SQP解法

	SUBROUTINE INITIA(N,ME,M,NAMAX,A,B,C,G,X,V,MJ,JSET,NJ,R,NM,
     *			L,W,Q,QPINDX,ZERO)
	DOUBLE PRECISION A(NAMAX,M),B(M),C(N),G(NAMAX,N),L(NAMAX,N),W(N),
     *		Q(NAMAX,N),R(NAMAX,N),NM(NAMAX,N),X(N),V(M),ZERO
	INTEGER JSET(M),QPINDX
	LOGICAL MJ(M),TRANS
	QPINDX=0
*----------Compute the Cholesky decomposition of G.
	CALL CHODEC(N,NAMAX,G,L,QPINDX,ZERO)
	IF (QPINDX.EQ.1)THEN
		WRITE(6,100)
100		FORMAT(//'*MATRIX G IS NOT POSITIVE DEFINE*')
		RETURN
	ENDIF
*----------Initialize the index set J.
	NJ=ME
	DO 1 I=1,M
		MJ(I)=.FALSE.
1	CONTINUE
	DO 2 I=1,ME
		JSET(I)=I
		MJ(I)=.TRUE.
2	CONTINUE
*----------Compute BsubJ=(inverse of L)*AsubJ.
*----------(matrix BsubJ is stored in array R.)
	DO 5 I=1,NJ
		DO 3 J=1,N
			W(J)=A(J,I)
3		CONTINUE
		TRANS=.FALSE.
		CALL FWDSUB(N,NAMAX,L,TRANS,W)
		DO 4 J=1,N
			R(J,I)=W(J)
4		CONTINUE
5	CONTINUE
*----------Compute the QR decomposition of matrix BsubJ.
	CALL QRDEC(N,NJ,NAMAX,R,Q,JRANK,ZERO)
	IF (JRANK.EQ.0)GOTO 6
	WRITE(6,110)JRANK
110	FORMAT(//'*THE EQUALITY CONSTRAINS ARE NOT INDEPENDENT*'
     *		/'*RANK DEFICIENCY DETECTED AT COLUM ',I5,'OF A')
	QPINDX=2
	RETURN
*----------Compute matrix NM.
6	DO 9 I=1,NJ
		DO 7 J=1,N
			W(J)=Q(J,I)
7		CONTINUE
		TRANS=.TRUE.
		CALL BWDSUB(N,NAMAX,L,TRANS,W)
		DO 8 J=1,N
			NM(J,I)=W(J)
8		CONTINUE
9	CONTINUE
	DO 12 I=1,N-NJ
		DO 10 J=1,N
			W(J)=Q(J,NJ+I)
10		CONTINUE
		TRANS=.TRUE.
		CALL BWDSUB(N,NAMAX,L,TRANS,W)
		DO 11 J=1,N
			NM(J,I+NJ)=W(J)
11		CONTINUE
12	CONTINUE
*----------Compute the initial solution.
	IF (NJ.EQ.0)THEN
		DO 13 J=1,N
			X(J)=0.0
13		CONTINUE
		GOTO 16
	ENDIF
	DO 14 I=1,NJ
		W(I)=B(I)
14	CONTINUE
	TRANS=.TRUE.
	CALL FWDSUB(NJ,NAMAX,R,TRANS,W)
	DO 15 J=1,N
		X(J)=0.0
		DO 15 I=1,NJ
			X(J)=X(J)+NM(J,I)*W(I)
15	CONTINUE
16	DO 17 I=NJ+1,N
		W(I)=0.0
		DO 17 J=1,N
			W(I)=W(I)+NM(J,I)*C(J)
17	CONTINUE
	DO 18 J=1,N
		DO 18 I=NJ+1,N
			X(J)=X(J)-NM(J,I)*W(I)
18	CONTINUE
*----------Compute the initial mnltipliers V.
	DO 19 I=1,M
		V(I)=0.0
19	CONTINUE
	IF (NJ.EQ.0)RETURN
	DO 20 I=1,NJ
		V(I)=B(I)
20	CONTINUE
	TRANS=.TRUE.
	CALL FWDSUB(NJ,NAMAX,R,TRANS,V)
	TRANS=.FALSE.
	CALL BWDSUB(NJ,NAMAX,R,TRANS,V)
	DO 21 I=1,NJ
		W(I)=0.0
		DO 21 J=1,N
		W(I)=W(I)+NM(J,I)*C(J)
21	CONTINUE
	TRANS=.FALSE.
	CALL BWDSUB(NJ,NAMAX,R,TRANS,W)
	DO 22 I=1,NJ
		V(I)=V(I)+W(I)
22	CONTINUE
	RETURN
	END