ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Solving a Nonlinear problem by A Simple Decomposition c based $SQP$ Algorithm c c c Mehdi Lachiheb, Faculte des Sciences de Gabes, c Cite Erriadh 6072, Gabes, Tunisie. c lachihebm@yahoo.ca c c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c min f(x) c gi(x)<= 0 i=1,mc c xi>= bi i=1,nx c c nx: number of variables c mc: number of constraints c m1: number of inquality constraints c m2: number of equality constraints c it: number of iterations c XL1:the lower bounds of the variables; c XK1: The starting point c dr: Accurency of the solution c bakcom.dat: file of results c badcom.dat: file of data, the starting point c stokcom.fd: file contains: integer mc,m1,m2,nx1,nx c parameter (nx= ,nx1=nx+1,m1= ,m2= ,mc= ) c Sd: the solution of quadratic problem c XK2:the solution of Nonlinear problem c itmin: the lower number where we use the procedure of approximation c extreme point set, itmin>1 c nbi: total number of extreme points c npe: total number of pivots c nitmax: maximal number of iterations program SEQUENTIELLE USE MSIMSL USe portlib $debug include"stokcom.fd" integer it,code,nu1,nbi,iposv(mc),npe,MU1(nx1,nx+1),itpq, + MNN,IOUT,IFAIL,IPRINT,LWAR,LIWAR,IWAR(nx),NACT,IACT(mc+nx), + itmin, nitmax real*8 GRADF(nx),GRADC(mc,nx),GRADC1(mc,nx),HF1(nx,nx),FO,FC(mc), + XK1(nx),XK2(nx),SD(nx),LK3(mc+nx),GRADL1(nx),temp, + GRADL2(nx),DL(nx),DX(nx),Xscale(nx),C(nx),B(mc), + epsfcn,CP1(mc+nx),CP2(mc+nx),alp,test,ST(3),dr,XXL(nx), + XXU(nx),BXL(nx+m1),BXU(nx+m1),BXL1(nx+m1),testg,XU(nx), + XL(nx),LK1(mc+2*nx),eps,rrri,rrrs,WAR(9000000),AQ(mc+nx,nx), + BQ(mc+nx),DIAG REAL(4) Te,VVA(2) external OBJE,CONTR COMMON /CMACHE/EPS eps=1.d-4 nitmax=60 itmin=400 Lwar=900000 Liwar=nx temp=timef() it=1 itpq=1 nbi=0 npe=0 dr=1.D-4 c starting point and hessian matrix call INIT(XK1,HF1,LK3,XXL,XXU,Xscale,epsfcn) c objectif function call OBJE(nx,XK1,FO) c contrainte function call CONTR(mc,nx,XK1,FC) c Gradient of objectif function call DCDGRD(OBJE,nx,XK1,xscale,epsfcn,GRADF) open(14,file='bakcom.dat',status='old') rewind(14) c Jacobian matrix of constaint functions c call DFDJAC(CONTR,mc,nx,XK1,xscale,FC,epsfcn,GRADC,mc) 99 continue c Formulation of quadratic problem call FSPQ(HF1,GRADF,GRADC,XK1,FC,XXL,XXU,BXL,BXU,C,B) c résolution du problème quadratique call RSPQ(HF1,C,GRADC,B,BXL,BXL1,BXU,MU1,SD,LK3,nu1,code,it,nbi, + npe,itmin) do 201 i=1,nx+m1 BXL1(i)=BXL(i) 201 continue call DCP(CP1,CP2,LK3,it) call DALPHA(ST,alp) call TRAX(XK1,XK2,SD,alp) if (TEST(SD,nx).LE.dr) goto 98 call GRADLAG(GRADF,GRADC,LK3,GRADL1) call OBJE(nx,XK2,FO) call CONTR(mc,nx,XK2,FC) call DCDGRD(OBJE,nx,XK2,xscale,epsfcn,GRADF) call DFDJAC(CONTR,mc,nx,XK2,xscale,FC,epsfcn,GRADC,mc) call GRADLAG(GRADF,GRADC,LK3,GRADL2) call CDVV(XK1,XK2,DX) call CDVV(GRADL1,GRADL2,DL) call STOCV(XK1,XK2,CP1,CP2) call MAJH(HF1,DX,DL) it=it+1 c maximal number of iteration if(IT.LE.nitmax) goto 99 98 continue call OBJE(nx,XK2,FO) call CONTR(mc,nx,XK2,FC) write(14,*) 'number of iteration ',it write(14,99997) XK2 write(14,*)'objectif function', FO write(14,99998) FC write(14,*)'number of extreme points',nbi write(14,*)'number of pivots',npe temp=timef() write(14,*)'temp',temp Te=DTIME(VVA) CALL UMACH (2, NOUT) WRITE (14,*) 'CPU time (seconds) = ',Te WRITE (14,*) ' time of user (seconds) = ',VVA(1) WRITE (14,*) 'system time (seconds) = ',VVA(2) close(14) 99997 FORMAT (' La solution est', /, + 3(5X,5F10.6,/),/) 99998 FORMAT (' constraint function', /, + 3(5X,5F10.6,/),/) end c ****************************************************************** c c initialisation of (LK1),XK1,HF1,Xscale,epsfcn c c ****************************************************************** subroutine INIT(XK1,HF1,LK3,XXL,XXU,Xscale,epsfcn) include"stokcom.fd" real*8 XK1(nx),HF1(nx,nx),LK3(mc+nx),Xscale(nx),epsfcn, + XXL(nx),XXU(nx) epsfcn=0.D0 open(12,file='badcom.dat',status='old') rewind(12) read(12,*)(XK1(i),i=1,nx) close(12) do 20 i=1,nx do 10 j=1,nx HF1(i,j)=0.0D0 10 continue HF1(i,i)=1.0D0 Xscale(i)=1.D0 XXL(i)=0.1d0 XXU(i)=1.D30 20 continue do 15 i=1,mc+nx LK3(i)=0.D0 15 continue return end c ********************************************************************** c c ****************determination of lagrangien gradient*************** c c ********************************************************************** subroutine GRADLAG(GRADF,GRADC,LK3,GRADL1) include"stokcom.fd" real*8 GRADF(nx),GRADC(mc,nx),LK3(mc+nx),GRADL1(nx),rs do 15 i=1,nx rs=0.D0 do 14 j=1,mc rs=rs+(LK3(j)*GRADC(j,i)) 14 continue GRADL1(i)=rs+GRADF(i)+LK3(mc+i) 15 continue return end c ************************************************************** c formulation of quadratic problem c ************************************************************** subroutine FSPQ(HF1,GRADF,GRADC,XK1,FC,XXL,XXU,BXL,BXU,C,B) include"stokcom.fd" real*8 HF1(nx,nx),GRADF(nx),GRADC(mc,nx),FC(mc),C(nx), + B(mc),zs,XK1(nx),XXL(nx),XXU(nx),BXL(nx+m1),BXU(nx+m1) do 41 i=1,nx C(i)=GRADF(i) if (XXL(i).eq.-1.D30) then BXL(i)=-1.D2 else BXL(i)=XXL(i)-XK1(i) endif if (XXU(i).eq.1.D30) then BXU(i)=1.D2 else BXU(i)=XXU(i)-XK1(i) endif 41 continue do 46 i=1,mc B(i)=-FC(i) 46 continue return end c **************************************************************** c translation of solution c **************************************************************** subroutine trans(SD,XK1,nx) integer nx real*8 SD(nx),XK1(nx) do 45 i=1,nx SD(i)=SD(i)+(0.1D0-XK1(i)) 45 continue return end c **************************************************************** c Computation of Lagrange Multipliers by simplexe method c **************************************************************** subroutine CMDL(HF1,GRADF,GRADC,XK1,XXL,XXU,FC,SD,LK3) include"stokcom.fd" integer Irtype(nx),k real*8 HF1(nx,nx),GRADF(nx),GRADC(mc,nx),FC(mc),SD(nx), + dw,A(nx,mc+2*nx),zd,BU1(nx),BU2(nx),XLB(mc+2*nx),XUB(mc+2*nx), + C(mc+2*nx),obj,DSOl(nx),fw,fk,LK3(mc+nx),LK4(mc+2*nx),XK1(nx), + XXL(nx),XXU(nx) dw=1.D-5 do 74 i=1,m1 XLB(i)=0.D00 XUB(i)=-1.D30 C(i)=0.D0 fw=0.D0 do 13 j=1,nx fw=fw+(GRADC(i,j)*SD(j)) 13 continue fk=fw+fc(i) if(DABS(FC(i)+fw).GE.dw) then XUB(i)=0.D0 endif 74 continue do 72 i=m1+1,mc XUB(i)=-1.D30 XLB(i)=1.D30 C(i)=0.D0 72 continue do 76 i=1,mc do 71 j=1,nx A(j,i)=GRADC(i,j) 71 continue 76 continue do 21 i=1,nx zd=0.D0 do 22 j=1,nx zd=zd+(HF1(i,j)*SD(j)) 22 continue BU1(i)=-zd-GRADF(i) BU2(i)=BU1(i) Irtype(i)=0 21 continue do 11 i=1,nx K=mc+i C(k)=0.D0 A(i,k)=1.D0 if(DABS(XXL(i)-SD(i)).GE.dw) then XLB(k)=0.D0 else XLB(k)=1.D30 endif XUB(k)=0.D0 11 continue do 12 i=1,nx K=mc+nx+i C(k)=1.D0 A(i,k)=1.D0 XUB(k)=-1.D30 XLB(k)=0.D0 12 continue call DDLPRS(nx,mc+2*nx,A,nx,BU1,BU2,C,Irtype,XLB,XUB,obj,LK4,DSOL) if (IERCD().EQ.3) goto 90 do 14 i=1,nx+mc LK3(i)=LK4(i) 14 continue 90 write(6,*)'non' return end c **************************************************************** c determination of penalite(CP2) coefficient c **************************************************************** subroutine DCP(CP1,CP2,LK3,it) include"stokcom.fd" integer it real*8 CP1(mc+nx),CP2(mc+nx),LK3(mc+nx),zs if(it.EQ.1) then do 12 i=1,mc CP2(i)=DABS(LK3(i)) 12 continue else do 15 i=1,mc zs=5.D-1*(DABS(LK3(i))+CP1(i)) if(zs.GE.DABS(LK3(i))) then CP2(i)=zs else CP2(i)=DABS(LK3(i)) endif 15 continue endif return end subroutine IQFP(XK1,CP2,SD,ST) include"stokcom.fd" real*8 XK1(nx),SD(nx),CP2(mc+nx),VC1(mc),ca1,ca(3),X1(nx),T(3,3), + ST(3) write(6,*)'XK1',xk1 do 21 i=1,nx X1(i)=XK1(i)+(1.1D0*SD(i)) 21 continue call OBJE(nx,X1,ca1) call CONTR(mc,nx,X1,VC1) ca(1)=ca1 do 15 i=1,m1 if(VC1(i).GE.0) then ca(1)=ca(1)+(CP2(i)*VC1(i)) endif 15 continue do 97 i=m1+1,mc ca(1)=ca(1)+(CP2(i)*DABS(VC1(i))) 97 continue do 16 i=1,nx X1(i)=XK1(i)+(0.9D0*SD(i)) 16 continue call OBJE(nx,X1,ca1) call CONTR(mc,nx,X1,VC1) ca(2)=ca1 do 17 i=1,m1 if(VC1(i).GE.0.D0) then ca(2)=ca(2)+(CP2(i)*VC1(i)) endif 17 continue do 94 i=m1+1,mc ca(2)=ca(2)+(CP2(i)*DABS(VC1(i))) 94 continue do 19 i=1,nx X1(i)=XK1(i)+SD(i) 19 continue call OBJE(nx,X1,ca1) call CONTR(mc,nx,X1,VC1) ca(3)=ca1 do 18 i=1,m1 if(VC1(i).GE.0.D0) then ca(3)=ca(3)+(CP2(i)*VC1(i)) endif 18 continue do 92 i=m1+1,mc ca(3)=ca(3)+(CP2(i)*DABS(VC1(i))) 92 continue T(1,1)=1.D0 T(2,1)=1.D0 T(3,1)=1.D0 T(1,2)=1.1D0 T(2,2)=0.9D0 T(3,2)=1.D0 T(1,3)=1.1D0*1.1D0 T(2,3)=0.9D0*0.9D0 T(3,3)=1.D0 call DLSLRG(3,T,3,ca,1,ST) return end c *************************************************************** c determination of alpha (Xk+1 =Xk + a.d) c *************************************************************** subroutine DALPHA(ST,alp) real*8 ST(3),alp c if (DABS(ST(3)).LE.0.d-8) then alp=1.D0 c else c alp=-ST(2)/(2.D0*ST(3)) c endif return end c ************************************************************ c compuation Xk+1 c ************************************************************ subroutine TRAX(XK1,XK2,SD,alp) include"stokcom.fd" real*8 XK1(nx),XK2(nx),SD(nx),alp do 12 i=1,nx XK2(i)=XK1(i)+(alp*SD(i)) 12 continue return end c ************************************************************** c stok vectors Xk+1,CPk+1,GRADL c ************************************************************** subroutine STOCV(XK1,XK2,CP1,CP2) include"stokcom.fd" real*8 XK1(nx),XK2(nx),CP1(mc+nx),CP2(mc+nx) do 45 i=1,nx XK1(i)=XK2(i) 45 continue do 46 i=1,mc CP1(i)=CP2(i) 46 continue return end c ******************************************************************* c difference of tow vectors c ******************************************************************* subroutine CDVV(V1,V2,V3) include"stokcom.fd" real*8 V1(nx),V2(nx),V3(nx) do 75 i=1,nx V3(i)=V2(i)-V1(i) 75 continue return end c ****************************************************************** c Adapted the Hessian matrix c ****************************************************************** subroutine MAJH(HF1,DX,DL) include"stokcom.fd" real*8 HF1(nx,nx),DX(nx),DL(nx),V1(nx),V2(nx),a1,a2,a3,a4,O,DCPVV external DCPVV write(6,*)'MISE A JOUR' a1=DCPVV(DX,DL,nx,nx) call DCPMV(HF1,DX,V1,nx,nx,nx) a2=DCPVV(DX,V1,nx,nx) c ******* calcule de O ************* if(a1.GE.(2.D-1*a2)) then O=1.D0 else O=(8.D-1*a2)/(a2-a1) endif c ******* calcule de V2=Odx+(1-O)Bdl ******* do 47 i=1,nx V2(i)=(O*DL(i))+((1.D0-O)*V1(i)) 47 continue a1=DCPVV(DX,V2,nx,nx) do 18 i=1,nx do 17 j=1,nx a3=V2(i)*V2(j) a4=V1(i)*V1(j) HF1(i,j)=HF1(i,j)-(a4/a2)+(a3/a1) 17 continue 18 continue write(6,*)'FIN MISE A JOUR' c CALL UMACH (2, NOUT) return end c ****************************************************************** c convergence test c ****************************************************************** real*8 function test(SD,nx) integer nx real*8 SD(nx) test=DABS(SD(1)) do 41 i=2,nx if(DABS(SD(i)).GE.test) then test=DABS(SD(i)) endif 41 continue return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc C Solving the quadratic problem ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine RSPQ(B,S,A,B1,XLB,BXL1,XUB,MU1,XS2,MG,nu1,code,it,nbi, + npe,itmin) $debug USE PORTLIB include"stokcom.fd" real*8 M(nx1,nx1),B(nx,nx),C(nx1),S(nx),E1(nx+m1),U(nx1), + B1(mc),XLB(nx+m1),XUB(nx+m1),XS1(nx),XS2(nx),DS(mc),B2(mc+1), + XS3(nx),MU(nx1,nx),CSM(nx+mc),obj,rs,DCPVV,X(nx1),W(nx), + A1(mc,nx),dmx,DS1(mc+1),BL(mc),BU(mc),XSO(nx+m1), + A(mc,nx),E(nx),CSM1(nx),BXL1(nx+m1),xxr,SMM(nx+1), + MG(mc+nx),temp,MB(mc+1,mc),rrri,Cc(nx+mc),SM(mc+1),dvaf1, + dvaf2 integer nm,nu,nv,Irtype(mc),Irtype1(mc+1),bor,ok,code,nq,nim + ,nu1,it,nbi,Izrov(nx+mc),IPOSV(mc),ICASE,itn,ns1,ns2, + nl1,nll1,L1(mc),LL1(mc),iposeq(m1),npe,iin,ind, + iposve(m1),MU1(nx1,nx+1),Vice(m2),cpi,Lp1(mc),Lc(mc), + Lp2(mc),nlc,MBs(mc+1,mc),Nvar(nx),As(mc+1,nx),isdp,itmin external DCPVV c nu1: final number of extreme points c c nu : number of extreme points c if isdp=0 use the primal simplex if isdp=1 use the dual simplex isdp=0 MP=mc+2 NP=nx+1 ind=1 cpi=0 do 91 i=1,m2 BU(i)=B1(i) BL(i)=B1(i) rrri=0.d0 do j=1,nx rrri=rrri-(A(i,j)*XLB(j)) enddo SM(m1+i)=B1(i)+rrri irtype(i)=0 91 continue ns1=0 ns2=0 do 92 i=m2+1,mc BU(i)=B1(i) BL(i)=Bu(i) rrri=B1(i) do j=1,nx rrri=rrri-(A(i,j)*XLB(j)) enddo if( rrri.GE.0.d0) then ns1=ns1+1 SM(ns1)=rrri iposeq(ns1)=i iposve(i-m2)=ns1 do j=1,nx A1(ns1,j)=A(i,j) enddo else SM(m1-ns2)=-rrri iposeq(m1-ns2)=i iposve(i-m2)=m1-ns2 do j=1,nx A1(m1-ns2,j)=-A(i,j) enddo ns2=ns2+1 endif irtype(i)=0 92 continue Irtype1(m1+m2+1)=3 code=0 dmx=1.D-5 do i=1,m2 do j=1,nx A1(i+m1,j)=A(i,j) enddo enddo do 103 j=1,m1 XLB(nx+j)=0.d0 XUB(nx+j)=-1.D30 103 continue do i=m1+1,mc VICE(i-m1)=0 if (SM(i).LE.0.d0) then VIce(i-m1)=1 do j=1,nx+1 A1(i,j)=-A1(i,j) enddo SM(i)=-SM(i) endif enddo if (it.LE.itmin) then call DCPMV(B,XS2,E,nx,nx,nx) do 61 i=1,nx E(i)=E(i)+S(i) Cc(i)=E(i) 61 continue do i=nx+1,nx+mc Cc(i)=0.d0 enddo itn=0 call simdual(A1,As,MB,MBs,Cc,SM,SMM,Iposv,L1,lp1,Izrov,LL1,nll1, + nl1,nl3,LC,Lp2,Nvar,nlc,Icase,itn,ns1,ns2,m2,mc,nx,npe,isdp) c temp=timef() c write(*,*)'temp',temp c call DDLPRS(mc,nx+m1,A1,mc,BL,BU,E1,Irtype,XLB,XUB,obj,XSO,DS) c WRITE(*,*) 'SIMPLEX ENDED ' c write(*,*)'ICASE',icase c if (IERCD().EQ.3) goto 110 nbi=nbi+1 do i=1,nx+m1 XSO(i)=XLB(i) enddo if(isdp.eq.0) then do i=1,nl1 if (iposv(i).LE.nx) XSO(iposv(i))=SM(i)+XLB(iposv(i)) enddo else do i=1,nl1 ii=L1(i) if (iposv(ii).LE.nx) XSO(iposv(ii))=SMM(i)+XLB(iposv(ii)) enddo endif iin=0 do 115 i=1,nl3 if (izrov(i).le.(nx+m1)) then iin=iin+1 if(izrov(i).LE.nx) then MU1(1,iin+1)=Izrov(i) else MU1(1,iin+1)=(iposeq(izrov(i)-nx)+nx)-m2 endif endif 115 continue MU1(1,1)=iin do 104 i=1,nx XS1(i)=XSO(i) 104 continue c ****** initialisation de matrice M ************ nm=1 M(1,1)=0.D0 do 22 i=1,nx rs=0.D0 do 23 j=1,nx rs=rs+B(j,i)*XS1(j) 23 continue M(1,1)=M(1,1)+rs*XS1(i) 22 continue c *********** initialisation de matrice MU ******** nu=1 nu1=1 nv=0 do 25 i=1,nx MU(1,i)=XS1(i) 25 continue c ********** initialisation de vecteur C ********* C(1)=0.D0 do 26 i=1,nx C(1)=C(1)+XS1(i)*S(i) 26 continue else nm=0 do 120 k=1,nu1 do i=1,nx+mc CSM(i)=0.d0 enddo iin=MU1(k,1) do i=1,iin if(MU1(k,i+1).LE.nx) then CSM(MU1(k,i+1))=1.d0 else j=iposve((MU1(k,i+1)-nx)) CSM(j+nx)=1.d0 endif enddo if (k.eq.1) then itn=0 else itn=1 endif call simdual(A1,As,MB,MBs,CSM,SM,SMM,Iposv,L1,lp1,Izrov,LL1,nll1, + nl1,nl3,LC,Lp2,Nvar,nlc,Icase,itn,ns1,ns2,m2,mc,nx,npe,isdp) do i=1,nx+m1 XSO(i)=XLB(i) enddo if(isdp.eq.0) then do i=1,nl1 if (iposv(i).LE.nx) XSO(iposv(i))=SM(i)+XLB(iposv(i)) enddo else do i=1,nl1 ii=L1(i) if (iposv(ii).LE.nx) XSO(iposv(ii))=SMM(i)+XLB(iposv(ii)) enddo endif nbi=nbi+1 if (IERCD().EQ.3) goto 110 iin=0 do i=1,nl3 if (izrov(i).le.(nx+m1)) then iin=iin+1 if(izrov(i).LE.nx) then MU1(k,iin+1)=Izrov(i) else MU1(k,iin+1)=(iposeq(izrov(i)-nx)+nx)-m2 endif endif enddo MU1(k,1)=iin do 123 i=1,nx XS1(i)=XSO(i) MU(k,i)=XS1(i) 123 continue call CMAAVE(M,B,MU,XS1,X,W,k,nm,nx) call CCAAVE(C,S,XS1,k,nx,nm) 120 continue endif nim=2 100 continue call initialise (U,nu,nx,nx1) dvaf2=dvaf1 call resoudre(M,C,U,dvaf1,nu,nm,bor) if (bor.EQ.0) then code=1 goto 110 else call CCSM(B,MU,U,S,XS2,CSM,nm) do i=1,mc CSM(nx+i)=0.d0 enddo itn=0 c ****************************************************** dmx=1.D-5 itn=1 do i=1,nx CSM(i)=CSM(i) enddo call simdual(A1,As,MB,MBs,CSM,SM,SMM,Iposv,L1,lp1,Izrov,LL1,nll1, + nl1,nl3,LC,Lp2,Nvar,nlc,Icase,itn,ns1,ns2,m2,mc,nx,npe,isdp) nbi=nbi+1 do i=1,m1+nx XSO(i)=XLB(i) enddo if(isdp.eq.0) then do i=1,nl1 if (iposv(i).LE.nx) XSO(iposv(i))=SM(i)+XLB(iposv(i)) enddo else do i=1,nl1 ii=L1(i) if (iposv(ii).LE.nx) XSO(iposv(ii))=SMM(i)+XLB(iposv(ii)) enddo endif obj=0.d0 do i=1,nx obj=obj+(XSO(i)*CSM(i)) enddo do 108 i=1,nx XS1(i)=XSO(i) CSM1(i)=CSM(i) 108 continue if (IERCD().EQ.3) goto 110 c WRITE(*,*) '****SIMPLEX ENDED' c if (IERCD().EQ.1) then c call ajouter(A,A1,CSM,B1,B2,nx,mc) c nq=mc+1 c call DDLPRS(nq,nx,A1,nq,B2,B2,CSM,Irtype1,XLB,XUB,obj,XS3,DS1) c B2(mc+1)=-110000.D0 c call DDLPRS(nq,nx,A1,nq,B2,B2,CSM,Irtype1,XLB,XUB,obj,XS1,DS1) c k=nm+1 c ok=0 c call normee(XS1,XS3) c WRITE (24,99999) (XS1(I),I=1,nx) c call EVNU(M,MU,U,C,nu,nv,nx,nm) c call CMAAVE(M,B,MU,XS1,X,W,k,nm,nx) c call CCAAVE(C,S,XS1,k,nx,nm) c call AVEMU(MU,XS1,nm,nx,nu,nv,k,ok) if((DABS(obj-DCPVV(CSM1,XS2,nx,nx)).LE.dmx).or.(ind.eq.0)) then c if((DABS(obj-DCPVV(CSM1,XS2,nx,nx)).LE.dmx).or.(ind.eq.0).or. c + (DABS(dvaf2-dvaf1).le.1.d-10)) then write(14,*)'obj',obj,DCPVV(CSM1,XS2,nx,nx) call EVNU(M,MU1,MU,U,C,nu,nv,nx,m1,nm) nu1=nu goto 110 else call EVNU(M,MU1,MU,U,C,nu,nv,nx,m1,nm) k=nu+1 ok=1 call CMAAVE(M,B,MU,XS1,X,W,k,nm,nx) call CCAAVE(C,S,XS1,k,nx,nm) call AVEMU(MU,XS1,nm,nx,nu,nv,k,ok) iin=0 do i=1,nx if (izrov(i).le.(nx+m1)) then iin=iin+1 if(izrov(i).LE.nx) then MU1(nu,iin+1)=Izrov(i) else MU1(nu,iin+1)=(iposeq(izrov(i)-nx)+nx)-m2 endif endif enddo MU1(nu,1)=iin nu1=nu endif endif write(14,*)'obj',obj,DCPVV(CSM1,XS2,nx,nx) nim=nim+1 goto 100 c110 call CMLMDIV(A,CSM,Iposeq,Iposv,MG,nx,mc, c + m1,m2) c 110 Call CMLMDT(A1,Minv,Vnp,CSM,VICE,Iposeq,Iposv,MG,nx,mc,m1,m2, c + ns1,ns2,cpi) c110 Call CMLMD(A1,IPOSV,IZROV,VICE,iposeq,MG,CSM,nx,mc,m1,m2,ns1,ns2) 110 continue if(mc.LE.nx) then do i=1,nl1 LC(i)=L1(i) enddo nlc=nl1 do i=1,nl1 k=L1(i) if(iposv(k).GT.nx) then kk=iposv(k)-nx do j=1,nlc if(kk.eq.LC(j)) goto 63 enddo 63 continue nlc=nlc-1 do jj=j,nlc Lc(jj)=Lc(jj+1) enddo endif enddo write(6,*)'nlc',nlc endif call CMLMD(MB,MBs,IPOSV,IZROV,VICE,iposeq,MG,CSM,nl1, + L1,Lc,nlc,nx,mc,m1,m2,ns1,ns2) do i=1,nx if ( (dABS (XLB(i)- XS2(i)).LE.dmx).or.(dABS (XUB(i)- XS2(i)) + .LE.dmx)) then rrri=0.d0 do j=1,mc rrri= rrri+ (A(j,i)*MG(j)) enddo MG(mc+i)=-CSM(i)-rrri else rrri=0.d0 do j=1,mc rrri= rrri+ (A(j,i)*MG(j)) enddo endif enddo c Call CMLMD(A1,IPOSV,IZROV,VICE,iposeq,MG,CSM,nx,mc,m1,m2,ns1,ns2) c do i=m2+1,mc c rrri=0.d0 c do j=1,nx c rrri=rrri+(A(i,j)*XS2(j)) c enddo c write(24,*)'GRFC(',i,')=',rrri-b1(i),MG(i) c enddo c do i=1,nx c rrri=0.d0 c do j=1,mc c rrri=rrri+(MG(j)*A(j,i)) c enddo c write(24,*)'GRADL',rrri+CSM(i) c enddo write(14,*)'nbi=',nbi,'npe',npe end c c c c ******************************************************** c produt of tow vectors c ******************************************************** real*8 function DCPVV(XS1,XS2,nx,nx1) integer nx,nx1 real*8 XS1(nx1),XS2(nx1) DCPVV=0.D0 do 10 i=1,nx DCPVV=DCPVV+XS1(i)*XS2(i) 10 continue return end c c c c ********************************************************** c computation of H ( min Hx Ax.GE.B1) c ********************************************************** c c subroutine CCSM(B,MU,U,S,XS2,CSM,nm) include"stokcom.fd" integer nm real*8 B(nx,nx),MU(nx1,nx),XS2(nx),S(nx),CSM(nx+mc),U(nx1), + CSM1(nx) call DCPMV(MU,U,XS2,nm,nx,nx1) call DCPMV(B,XS2,CSM1,nx,nx,nx) do 40 i=1,nx CSM(i)=CSM1(i)+S(i) 40 continue return end c ******************************************************** c product vector and matrix c ******************************************************** c subroutine DCPMV(M,XS1,XS2,nm,nx,n) integer nm,nx,n real*8 M(n,nx),XS1(n),XS2(nx),ds do 100 i=1,nx ds=0.D0 do 10 j=1,nm ds=ds+M(j,i)*XS1(j) 10 continue XS2(i)=ds 100 continue return end c ***************************************************************** c c change the vecor C with ellimination an extreme point c c c ***************************************************************** c subroutine CCAEVE(C,i,nm,nx) integer nm,nx,i real*8 C(nx+1) do 14 j=i,nm-1 C(j)=C(j+1) 14 continue return end c **************************************************** c change the matrix M with eliminatin an extreme point c **************************************************** c subroutine CMAEVE(M,k,nm,nx) integer k,nm,nx,l real*8 M(nx+1,nx+1) do 45 i=k,nm-1 l=i+1 do 44 j=1,nm M(j,i)=M(j,l) 44 continue 45 continue do 48 i=k,nm-1 l=i+1 do 47 j=1,nm M(i,j)=M(l,j) 47 continue 48 continue nm=nm-1 return end c ********************************************************* c ellimination some vectors c change the matix (trMU.B.MU) c and the vector C c ********************************************************* subroutine EVNU(M,MU1,MU,U,C,nu,nv,nx,m1,nm) integer nu,nm,nx,nv,n,kl,np,nq,m1,MU1(nx+1,nx+1) real*8 M(nx+1,nx+1),MU(nx+1,nx),U(nx+1),C(nx+1),ds ds=1.D-12 kl=nu n=nm np=0 do 14 i=1,n if (DABS(U(i)).LE.ds) then nq=i-np call CCAEVE(C,nq,nm,nx) call CMAEVE(M,nq,nm,nx) if(i.LE.kl) then nu=nu-1 else nv=nv-1 endif do 12 j=nq,nm do 11 k=1,nx MU(j,k)=MU(j+1,k) 11 continue do 112 k=1,nx+1 MU1(j,k)=MU1(j+1,k) 112 continue 12 continue np=np+1 endif 14 continue end c c c ************************************************* c change MU with adding an extreme point c ************************************************* c subroutine AVEMU(MU,XS1,nm,nx,nu,nv,k,ok) integer k,nm,nx,nu,nv,ok real*8 MU(nx+1,nx),XS1(nx) if (ok.eq.1) then nu=nu+1 else nv=nv+1 endif do 25 i=nm,k+1,-1 do 24 j=1,nx MU(i,j)=MU(i-1,j) 24 continue 25 continue do 26 i=1,nx MU(nu,i)=XS1(i) 26 continue return end c c c c c ************************************************* c change M with adding an extreme point c ************************************************* c subroutine CMAAVE(M,B,MU,XS1,X,W,k,nm,nx) integer k,nm,nx real*8 M(nx+1,nx+1),B(nx,nx),MU(nx+1,nx),XS1(nx),X(nx+1),W(nx), + sr,var,DCPVV external DCPVV call DCPMV(B,XS1,W,nx,nx,nx) do 51 i=1,nm sr=0.D0 do 50 j=1,nx sr=sr+MU(i,j)*W(j) 50 continue X(i)=sr 51 continue do 61 i=nm,k,-1 do 60 j=1,nm var= M(i,j) M(i+1,j)=var 60 continue 61 continue do 71 i=nm,k,-1 do 70 j=1,nm+1 var= M(j,i) M(j,i+1)=var 70 continue 71 continue nm=nm+1 do 80 i=1,k-1 M(k,i)=X(i) M(i,k)=X(i) 80 continue M(k,k)=DCPVV(W,XS1,nx,nx) do 81 i=k+1,nm M(k,i)=X(i-1) M(i,k)=X(i-1) 81 continue return end c c c *********************************************************** c change C with adding an extreme point c *********************************************************** subroutine CCAAVE(C,S,XS1,k,nx,nm) c implicit real*8(D) integer k,nm,nx,l real*8 C(nx+1),XS1(nx),S(nx),DCPVV c external DCPVV do 14 i=nm,k+1,-1 l=i-1 C(i)=C(l) 14 continue C(K)=DCPVV(S,XS1,nx,nx) return end c ************************************ subroutine ajouter(A,A1,CSM,B1,B2,nx,mc) integer mc,nx real*8 A(mc,nx),A1(mc,nx),CSM(nx),B1(mc),B2(mc+1) do 10 i=1,nx do 11 j=1,mc A1(j,i)=A(j,i) 11 continue A1(mc+1,i)=CSM(i)*1 10 continue do 15 i=1,mc B2(i)=B1(i) 15 continue B2(mc+1)=-100000.D0 return end c **************************************** subroutine normee(XS1,XS3) include"stokcom.fd" real*8 XS1(nx),XS3(nx),z z=0.D0 do 16 i=1,nx XS1(i)=XS1(i)-XS3(i) if (XS1(i).LE.1D-7) then XS1(i)=0.D0 endif z=z+(XS1(i)*XS1(i)) 16 continue write(*,*)'z',z do 17 i=1,nx XS1(i)=XS1(i)/DSQRT(z) 17 continue return end c ******************************************** subroutine initialise (U,nu,nm,nx1) integer nm,nx1,nu real*8 U(nx1),rs rs=1.d0/nu do 14 i=1,nu U(i)=rs 14 continue do 15 i=nu+1,nm U(i)=1.D0 15 continue return end c ************************************* c computation the lagrange multipliers c *************************************** subroutine CMLMD(MB,MBs,IPOSV,IZROV,VICE,iposeq,MG,CSM,nl1, + L1,Lc,nlc,nx,mc,m1,m2,nn1,nn2) Real*8 MB(mc+1,mc), MG(mc+nx),CSM(nx+mc),rrin Integer Iposv(mc),Izrov(nx+mc),iposeq(m1),nx,mc,nne,nn1,nn2,m2, + VICE(m2),MBs(mc+1,mc),L1(mc),Lc(mc),nlc,m1 do i=1,nx+mc MG(i)=0.d0 enddo do i=1,nlc ii=Lc(i) rrin=0.d0 c write(6,*)'ii',ii do j=1,MBs(1,ii) k=MBs(j+1,ii) c write(6,*)'kkkk',k kk=L1(k) kkk=iposv(k) if(kkk.LE.nx) then rrin=rrin+(MB(k,ii)*CSM(kkk)) endif enddo if (ii.LE.m1) then if(ii.LE.nn1) then MG(iposeq(ii))=-rrin else MG(iposeq(ii))=rrin endif else if(VICE(ii-m1).eq.1) then MG(ii)=-rrin else MG(ii)=rrin endif endif enddo return end c ********************************************************** c computation the lagrange multipliers by c transformations c ********************************************************** subroutine CMLMDT(A1,Minv,Vnp,CSM,VICE,Iposeq,Iposv,MG,nx,mc, + m1,m2,ns1,ns2,cpi) integer nx,mc,m1,m2,ns1,ns2,cpi,Iposeq(m1),iposv(mc), + Vnp(10000),Vice(m2) Real*8 MG(mc+nx),Minv(mc,10000),CSM(nx+mc),MG1(mc),rrri, + A1(mc,nx) c write(14,*)'iposv',iposv,'*****cpi',cpi do i=1,mc if(IPOSV(i).LE.nx) then MG1(i)=CSM(iposv(i)) else MG1(i)=0.d0 endif enddo do i=1,mc+nx MG(i)=0.d0 enddo rrri=0.d0 do j=1,mc if(IPOSV(j).LE.nx) then rrri=rrri+MG1(j)*Minv(j,cpi) endif enddo write(*,*)'cpi',cpi write(*,*)'vnp',vnp(cpi) MG1(Vnp(cpi))=rrri do i=cpi-1,1,-1 rrri=0.d0 do j=1,mc rrri=rrri+(MG1(j)*Minv(j,i)) enddo MG1(Vnp(i))=rrri enddo do i=1,mc if(IPOSV(i).GT.nx) then MG1(iposv(i)-nx)=0.d0 endif enddo do i=1,m2 if(VIce(i).eq.1) then MG(i)=MG1(i+m1) else MG(i)=-MG1(i+m1) endif enddo do i=1,ns1 MG(iposeq(i))=-MG1(i) enddo do i=ns1+1,m1 MG(iposeq(i))=MG1(i) enddo return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccccccccccc computaton thelagange mutipliers by the invese of basis matrix ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine CMLMDIV(A,CSM,Iposeq,Iposv,MG,nx,mc, + m1,m2) integer nx,mc,m1,m2,Iposeq(m1),iposv(mc),nbca,Cac(mc),k, + xac(nx),l Real*8 MG(mc+nx),CSM(nx+mc),MG1(nx),rrri, + A(mc,nx),MI(nx,nx),BI(nx),rri c write(14,*)'iposv',iposv,'*****cpi',cpi do i=1,nx xac(i)=1 enddo do i=1,mc if(IPOSV(i).le.nx) then xac(iposv(i))=0 endif enddo nbca=mc do i=1,mc Cac(i)=1 enddo c write(14,*)'k',k do i=1,mc if(IPOSV(i).GT.nx) then k=iposv(i)-nx if(k.LE.m1) then Cac(Iposeq(k))=0 else Cac(k-m1)=0 endif nbca=nbca-1 endif enddo write(*,*)'nbca',nbca k=0 do i=1,nx if(xac(i).EQ.0) then k=k+1 BI(k)=CSM(i) endif enddo k=0 do i=1,mc if(Cac(i).EQ.1) then k=k+1 l=0 do j=1,nx if(xac(j).EQ.0) then l=l+1 MI(l,k)=A(i,j) endif enddo endif enddo call RSEL(MI,BI,MG1,nx,nbca) do i=1,mc+nx MG(i)=0.d0 enddo k=0 do i=1,mc if(Cac(i).EQ.1) then k=k+1 MG(i)=-MG1(k) c write(14,*)'mul',MG(i) endif enddo return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccc solution of linear system cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc Subroutine RSEL(M,V,X,nx1,n) Real*8 M(nx1,nx1),V(nx1),A(n,n),b(n),X(nx1),dmax Integer i,j,k,ip,n,nx1 do i=1,n do j=1,n A(i,j)=M(i,j) enddo b(i)=V(i) enddo do k=1,n-1 dmax=dabs(A(k,k)) ip=k do i=k+1,n if (dabs(A(i,k)).GT.dmax) then dmax=dabs(A(i,k)) ip=i endif enddo if (dmax.LE.1.d-15) then write(14,*)'Erreur' return endif do i=k,n dmax=A(ip,i) A(ip,i)=A(k,i) A(k,i)=dmax enddo dmax=b(ip) b(ip)=b(k) b(k)=dmax dmax=1.d0/A(k,k) do i=k,n A(k,i)=A(k,i)*dmax enddo b(k)=b(k)*dmax do i=k+1,n dmax=A(i,k) do j=k,n A(i,j)=A(i,j)-(dmax*A(k,j)) enddo b(i)=b(i)-(dmax*b(k)) enddo enddo if (dabs(A(n,n)).LE.1.d-15) then write(14,*)'erreur',A(n,n) return endif X(n)=b(n)/A(n,n) do i=n-1,1,-1 dmax=0.d0 do j=i+1,n dmax=dmax+(A(i,j)*X(j)) enddo X(i)=b(i)-dmax enddo return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC subroutine resoudre(M,C,U,dvaf,nu,nm,bor) $debug include"stokcom.fd" integer nu,nm,bor,nbit,arr,NACT,IACT(nu),ires real*8 m(nx1,nx1),c(nx1),u(nx1),U1(nx1),r,re,h(nx1,nx1), + d1(nx1),d2(nx1),d(nx1),df1,df2,w,dmx,dam,dz,uye, + dma,ds,dhx,DCPVV,don,DERIVF,e(nx1),SABD,SD,dnx, + BQ(nu,nu),AQ(nu+1,nu),CQ(nu),BsQ(nu+1),DIAG,SOL(nu), + ALamda(nu),f(nx1),dvaf external don,SABD,SD external DLSLSF,derivf,dcpvv ires=1 if (ires.eq.1) then do i=1,nu do j=1,nu BQ(i,j)=M(i,j) AQ(i+1,j)=0.d0 enddo CQ(i)=C(i) BSQ(i+1)=0 AQ(1,i)=1 AQ(1+i,i)=1 enddo BSQ(1)=1 c solving the quadratic problem by IMSL CALL DQPROG (nu,nu+1,1,AQ,nu+1, BSQ,CQ,BQ,nu, DIAG, + SOL, NACT, IACT, ALAMDA) do i=1,nu U(i)=sol(i) c write( 14,*)'sol(',i,')=',sol(i) enddo else c Solving the quadratic problem by Barrier method do 20 i=1,nu e(i)=1.D0 20 continue do 21 i=nu+1,nm e(i)=0.D0 21 continue ds=1.D-10 dnx=1.D-10 dz=1D-10 re=1.0D-1 r=1.D-1 bor=1 nbit=0 arr=0 df1=don(m,u,c,r,nm) 10 continue do 3 i=1,nm U1(i)=U(i) 3 continue nbit=nbit+1 df2=df1 r=r*re call hessien(m,h,u,r,nm) call gradient(m,u,d,c,r,nm) c write(*,*)'debut' call DLSLSF(nm,h,nx1,d,d1) call DLSLSF(nm,h,nx1,e,d2) c write(*,*)'fin' call wcalc(d1,d2,nu,w) call dcalc(d,d1,d2,w,nm) if ((SABD(d,nm).LE.dnx).or.(DABS(SD(d,nu)).GE.dz)) goto 7 call cdmxx(u,d,nm,dmx) call MDFC(m,u,d,dmx,c,r,dma,nm) call DDMA(m,c,d,u,dmx,dam,r,nm) c write (24,*) 'Dmx',dmx c write (24,*) 'Dma',dma c write (24,*) 'Dam',dam call change(u,d,dmx,dma,dam,nm) df1= don(m,u,c,r,nm) call bornee(u,df1,nu,nm,bor) if(bor.EQ.0) goto 1 7 continue if (dabs(df2).LT.1) then dhx=1 else dhx=dabs(df2) endif uye=1.D0 do 78 i=1,nm if (abs(U(i)).LE.uye) then uye= abs(U(i)) endif 78 continue c do 31 i=1,nu c write(24,*)'U(',i,')',U(i) c 31 continue if (((dabs(df2-df1)/dhx).LE.ds).and.(nu.GE.2)) then call securute(M,C,U,nu,arr) endif endif call DCPMV(m,u,f,nm,nx1,nx1) dvaf=(DCPVV(f,u,nm,nx1)/2.D0)+DCPVV(c,u,nm,nx1) 1 return end c ********************************** subroutine wcalc(d1,d2,nu,w) include"stokcom.fd" integer nu real*8 d1(nx1),d2(nx1),w,dz,dj dz=0.D0 dj=0.D0 do 10 i=1,nu dz=dz+d1(i) dj=dj+d2(i) 10 continue w=-dz/dj return end c *********************************************** subroutine hessien(m,h,u,r,nm) include"stokcom.fd" integer nm real*8 m(nx1,nx1),u(nx1),h(nx1,nx1),r do 14 i=1,nm do 13 j=1,nm h(j,i)=m(j,i) 13 continue 14 continue do 12 k=1,nm h(k,k)=h(k,k)+(r/(u(k)*u(k))) 12 continue return end c *********************************** subroutine gradient(m,u,g,c,r,nm) include"stokcom.fd" integer nm real*8 m(nx1,nx1),u(nx1),g(nx1),c(nx1),r c call DCPMV(m,u,g,nm,nx1,nx1) do 22 i=1,nm g(i)=0.d0 do 23 j=1,nm g(i)=g(i)+(m(i,j)*u(j)) 23 continue 22 continue do 12 i=1,nm g(i)=-(g(i)+c(i)-(r/u(i))) 12 continue return end c ****************************************** subroutine dcalc(d,d1,d2,w,nm) include"stokcom.fd" integer nm real*8 d(nx1),d1(nx1),d2(nx1),w do 10 i=1,nm d(i)=d1(i)+(w*d2(i)) 10 continue return end c *********************************************** real*8 function don(m,u,c,r,nm) include"stokcom.fd" integer nm real*8 m(nx1,nx1),u(nx1),c(nx1),f(nx1),r,dhe,DCPVV external DCPVV dhe=0.D0 call DCPMV(m,u,f,nm,nx1,nx1) do 14 i=1,nm dhe=dhe+DLOG(u(i)) 14 continue don=(DCPVV(f,u,nm,nx1)/2.D0)+DCPVV(c,u,nm,nx1)-(r*dhe) return end c *************************************** c c subroutine cdmxx(u,d,nm,dmx) include"stokcom.fd" integer nm real*8 u(nx1),d(nx1),dmx,dcs,dcf dcf=1.D30 do 12 i=1,nm if (d(i).LT.0.D0) then dcs=-u(i)/d(i) if (dcs.LT.dcf) then dcf=dcs endif endif 12 continue dmx=dcf return end c c c ************************************************* c subroutine DDMA(m,c,d,u,dmx,dam,r,nm) include"stokcom.fd" real*8 m(nx1,nx1),c(nx1),d(nx1),u(nx1),dmx,dam,r,er,ez,ee,ef,es, + ew,dea,deb integer nm er=0.D0 ez=0.D0 ee=1.D0 ef=1.D0 do 17 i=1,nm es=0.D0 ew=0.D0 do 18 j=1,nm es=es+m(j,i)*u(j) ew=ew+m(j,i)*(u(j)+dmx*d(j)) 18 continue c ee=ee+dlog(u(i)) ee=ee*u(i) ef=ef*(u(i)+dmx*d(i)) es=es/2.D0 ew=ew/2.D0 er=er+(es+c(i))*u(i) ez=ez+(ew+c(i))*(u(i)+dmx*d(i)) 17 continue c deb=(ez-er+(r*(ee-DLOG(dmx))))/dmx c deb=(ez-er)/dmx ef=(ee-ef)/dmx dea=er deb=(ez-dea)/dmx if(deb.EQ.0.D0) then deb=1.D-10 endif dam=(ee/ef)+(r/deb) c dam=dmx-(r/deb) return end c **************************************************************** real*8 function DERIVF(m,c,d,u,r,v,nm) include"stokcom.fd" integer nm real*8 m(nx1,nx1),c(nx1),u(nx1),d(nx1),f(nx1),v,r,dz,dx,dy,dww, + DCPVV external DCPVV dz=0.D0 do 14 i=1,nm dz=dz+(d(i)/(u(i)+(v*d(i)))) 14 continue call DCPMV(m,d,f,nm,nx1,nx1) dx=DCPVV(f,d,nm,nx1) dy=DCPVV(f,u,nm,nx1) dww=DCPVV(c,d,nm,nx1) DERIVF=dy+dww+v*dx-r*dz return end c *********************************************************** subroutine MDFC(m,u,d,dmx,c,r,dma,nm) include"stokcom.fd" integer nm real*8 m(nx1,nx1),u(nx1),d(nx1),c(nx1),b,r,a,z,ds,wq,dma,dmx, + DERIVF external DERIVF ds=1.D-8 a=0.D0 b=9.9D-1*dmx 11 continue z=(b+a)/2.D0 wq=DERIVF(m,c,d,u,r,z,nm) if (DERIVF(m,c,d,u,r,a,nm)*DERIVF(m,c,d,u,r,b,nm).GE.0.d0) + then if(DERIVF(m,c,d,u,r,a,nm).GE.0.d0) then dma=a else dma=b endif else if (DERIVF(m,c,d,u,r,a,nm)*wq.GE.0.d0) then a=z dma=a if (dabs(b-a).GT.ds) goto 11 else b=z dma=b if (dabs(b-a).GT.ds) goto 11 endif return end c ****************************************************** c c c subroutine change(u,d,dmx,dma,dam,nm) include"stokcom.fd" integer nm real*8 u(nx1),d(nx1),dmx,dma,dam,dz,d0 dz=dmx*9.99D-1 if ((dma.LT.dam).and.(dam.LT.dmx).and.(dam.LT.1.D0)) then d0=dam elseif (dz.LT.dma) then d0=dz else d0=dma endif do 14 i=1,nm u(i)=u(i)+d0*d(i) 14 continue return end c ************************************************************* real*8 function SABD(d,nm) include"stokcom.fd" integer nm real*8 d(nx1) SABD=0.D0 do 12 i=1,nm SABD=SABD+DABS(d(i)) 12 continue return end c ************************************************************* real*8 function SD(d,nu) include"stokcom.fd" integer nu real*8 d(nx1) SD=0.D0 do 12 i=1,nu SD=SD+d(i) 12 continue return end c *************************************************************** subroutine bornee(u,df1,nu,nm,bor) include"stokcom.fd" integer bor,nu,nm real*8 u(nx1),df1 bor=1 do 14 i=nu+1,nm if(u(i).GE.1.D15) then if(DABS(df1).GE.1.D15) then bor=0 else u(i)=1.D0 endif endif 14 continue return end c c c c c subroutine securute(M,C,U,nu,arr) include"stokcom.fd" integer nu,arr,k,l,q real*8 m(nx1,nx1),u(nx1),RM(nx1,nx1),C(nx1),RC(nx1),ren,w, + d1(nx1),d2(nx1),e(nx1),rs,rz,RU(nx1),vin ren=1.D-5 k=0 c U(nu)=1 do 33 i=1,nu if (U(i).GE.ren) then k=k+1 RC(k)=C(i) do 34 j=1,nu RM(k,j)=M(i,j) 34 continue endif 33 continue do 35 i=1,k e(i)=1.d0 35 continue q=0 do 41 i=1,nu if (U(i).GE.ren) then q=q+1 do 42 j=1,k RM(j,q)=RM(j,i) 42 continue endif 41 continue c write(*,*)'debut1' call DLSLSF(k,RM,nx1,RC,d1) c write(*,*)'fin1' c write(*,*)'debut2' call DLSLSF(k,RM,nx1,e,d2) c write(*,*)'fin2' rz=0.d0 rs=0.d0 do 36 i=1,k rz=rz+d1(i) rs=rs+d2(i) 36 continue w=(1+rz)/rs do 37 i=1,k RU(i)=w*d2(i)-d1(i) 37 continue k=0 do 38 i=1,nu if (U(i).GE.ren) then k=k+1 U(i)=RU(k) else U(i)=0.d0 endif 38 continue arr=50 return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Solving quadratic by goldfarb idnani dual method c =================================== c dualcom.inc is afile contain c implicit none c include 'dualpar.inc' c logical proto, lz, spd, lag c integer case, output , a , q c double precision u c common/drivpa/proto, lz, output, case, spd, lag, c x u(mx),a(mx) , q c character*15 outdat c common/case1/outdat c double precision sum , tau c common/case2/sum(4), tau c double precision vek1p, vek2p, vek1m, vek2m c common/case3/ vek1p(nv), vek2p(nv), vek1m(mx), vek2m(mx) c integer iq c double precision epsmac, rnorm c c common/parm/ epsmac,rnorm,iq c double precision xj,d,r,np c common /up/ xj(nv,nv),d(nv),r(mx,mx),np(nv) c save c ==================================================== c dualpar.inc is a file contain c maximum dimension primal = 300, dual = 500 c integer nv, mx c c parameter(nv = 300, mx = 500) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine driver(g, ci, ce, g0, ci0, ce0, x,n,m,p) include 'dualcom.inc' integer i,n,m,p external dual Real*8 g(n,n), ci(n,m), ce(n,p), g0(n) Real*8 ci0(m), ce0(p), x(n) character*15 RES1 c c----------------- c parameters for basic test c----------------- case = 1 proto = .false. spd = .true. lag = .false. c outdat= RES1 output=6 open(10,file='TEST.RES',status='old') rewind(10) c write(*,*)'la solution est' c open(10,file=outdat, status='unknown') c c-------------------- c start der iteration c-------------------- c c call dual(g, ci, ce, g0, ci0, ce0, x,n,m,p) close(10) end c/////////////////////////////////////////////////////////////// c///////////////////////////////////////////////////////////////////// subroutine dual(g, ci, ce, g0, ci0, ce0, x,n,m,p) include 'dualcom.inc' c parameters stored in common-areas c input c integer output: channel number for messages c logical proto: if true, any useful intermediate information is stored c in the file 10 (resname) c logical spd : if false, g isn't positive definite c logical lag : if true, ce , the matrix of gradients of equlity constraints c doesn't have full rank c c double precision g,ci,ce,g0,ci0,ce0 : matrices and vectors which describe c the problem parameters c a (of length q) gives the indices in the working set ( negative for equality constraints) c and u is the vector of multpliers integer i,j,l,k,ip,me,mi,n,m,p Real*8 g(n,n), ci(n,m), ce(n,p), g0(n) Real*8 ci0(m), ce0(p), x(n) double precision diagg(nv),bigbd,t1,t2,psi,c1,c2,yl,zz,su,ss double precision f,t,s(mx),z(nv),vr(mx) double precision xold(nv),uold(mx) integer iai(mx),iaexcl(mx),aold(mx),ifail double precision skal1 external skal1 cc bigbd=1.d20 me=p mi=m ifail=0 q=0 c*** this stands for "infinity" c1=0.d0 do i=1,n diagg(i)=g(i,i) c1=c1+diagg(i) enddo c2=0.d0 call chol(g,n,g,ifail) c write(*,*)'test2',ifail if(ifail .ne. 0) then spd = .false. return endif cc initialize matrices "j" (that is xj) and r do i=1,n if(i .le.n) d(i)=0.d0 do j=1,n r(i,j)=0.d0 enddo enddo rnorm=1.d0 do i=1,n if(i .gt. 1) d(i-1)=0.d0 d(i)=1.d0 call leftel(g,d,z,yl,n) do j=1,n xj(i,j)=z(j) enddo c2=c2+z(i) enddo cc c1*c2 gives a rough estimate for cond(g) cc compute unconstrained minimizer x of f and f_min call cholso(g,n,g0,x) do i=1,n x(i)=-x(i) enddo f=0.5*skal1(1,n,g0,x) if (proto) then write(*,fmt='('' unconstrained solution : '')') write(10,fmt='(5(d14.7,1x))') (x(i),i=1,n) endif cc add equality constraints to the working set iq=0 do i=1,me if (proto) then write(10,fmt='(i4,''-th equality constraint: add '')') f i endif do j=1,n np(j)=ce(j,i) enddo call zup(z,n) if ( iq .ne. 0 ) call rup(vr) t2 = 0.d0 zz = skal1(1,n,z,z) if ( zz .ne. 0.d0 ) 1 t2=(-skal1(1,n,np,x)-ce0(i))/skal1(1,n,z,np) do k =1,n x(k)=x(k)+t2*z(k) enddo u(iq+1)=t2 do k=1,iq u(k)=u(k)-t2*vr(k) enddo f=f+t2**2*.5d0*skal1(1,n,z,np) if (proto) then write(10,*) ' vector z' write(10,fmt='(5(d14.7,1x))') (z(k),k=1,n) write(10,*) ' vector u' write(10,fmt='(5(d14.7,1x))') (u(k),k=1,iq+1) write(10,*) ' vector x' write(10,fmt='(5(d14.7,1x))') (x(k),k=1,n) endif a(i)=-i call addc(ifail,n) if ( ifail .ne. 0 ) then if (proto) then write(*,*) ' equality constraints linearly dependent ' endif lag = .true. return endif enddo cc set iai=k\a do i=1,mi iai(i)=i enddo cc step 1 50 continue do i=me+1,iq ip=a(i) iai(ip)=0 enddo cc s(x)=ci(trans)*x+ci0 ss=0.d0 psi=0.d0 ip=0 do i=1,mi iaexcl(i)=1 su=0.d0 do j=1,n su=su+ci(j,i)*x(j) enddo su=su+ci0(i) s(i)=su psi=psi+min(0.d0,su) enddo if (proto) then write(10,*) f '***********************************************************' write(10,*) ' vector g(x)' write(10,fmt='(3(i4,1x,d14.7,1x))') (i,s(i),i=1,mi) write(10,*) ' sum of infeasibilities =',psi endif if ( dabs(psi) .le.dble(mi)*epsmac*c1*c2*100.d0 ) then q=iq return endif do i=1,iq uold(i)=u(i) aold(i)=a(i) enddo do i=1,n xold(i)=x(i) enddo 60 continue do i=1,mi if(s(i) .lt. ss .and. iai(i) .ne. 0 .and. iaexcl(i) .ne.0)then ss=s(i) ip=i endif enddo if(ss .ge. 0.d0) then q=iq return endif do i=1,n np(i)=ci(i,ip) enddo u(iq+1)=0.d0 a(iq+1)=ip if (proto) then write(10,*) ' introduce restriction p=',ip write(10,*) ' working set ' write(10,fmt='(10(i5,1x))') (a(i),i=1,iq+1) write(10,*) ' gradient of restriction p' write(10,fmt='(5(d14.7,1x))') (np(i),i=1,n) endif 100 continue cc step 2a if (proto) then write(10,*) f '------------------------------------------------------------' endif call zup(z,n) if(iq .ne. 0) call rup(vr) if (proto) then write(10,*) ' vector d ' write(10,fmt='(5(d14.7,1x))') (d(i),i=1,n) write(10,*) ' primal correction z' write(10,fmt='(5(d14.7,1x))') (z(i),i=1,n) endif if ( iq .gt. 0)then if (proto) then write(10,*) ' dual correction v ' write(10,fmt='(5(d14.7,1x))') (vr(i),i=1,iq) endif endif l=0 t1=bigbd do k=me+1,iq if(vr(k) .gt. 0.d0) then if(u(k)/vr(k) .lt. t1) then t1=u(k)/vr(k) l=a(k) endif endif enddo cc |z|=0? zz=skal1(1,n,z,z) if(zz .ne. 0.d0 ) then t2=-s(ip)/skal1(1,n,z,np) else t2=bigbd endif t=dmin1(t1,t2) if (proto) then write(10,fmt='('' stepsizes: t1='',d14.7, x '' t2='',d14.7)') t1,t2 endif if(t .ge. bigbd ) then write(*,*) ' dual unbounded ' q=iq return endif if(t2 .ge. bigbd ) then do k=1,iq u(k)=u(k)+t*(-vr(k)) enddo u(iq+1)=u(iq+1)+t iai(l)=l call delc(l,n,m,p) if (proto) then write(10,*) ' partial step in the duals only ' write(10,*) ' working set ' write(10,fmt='(10(i4,1x))') (a(i),i=1,iq) write(10,*) ' vector u' write(10,fmt='(5(d14.7,1x))')( u(i),i=1,iq) endif goto 100 endif do k=1,n x(k)=x(k)+t*z(k) enddo f=f+t*skal1(1,n,z,np)*(0.5*t+u(iq+1)) if (proto) then write(10,*) ' current solution x ' write(10,fmt='(5(d14.7,1x))') (x(i),i=1,n) write(10,*) ' primal function value ',f endif do k=1,iq u(k)=u(k)+t*(-vr(k)) enddo u(iq+1)=u(iq+1)+t if(t .eq.t2) then call addc(ifail,n) if ( ifail .ne. 0 ) then iaexcl(ip)=0 call delc(ip,n,m,p) do i=1,m iai(i)=i enddo do i=1,iq a(i)=aold(i) iai(a(i))=0 u(i)=uold(i) enddo do i=1,n x(i)=xold(i) enddo goto 60 else iai(ip)=0 endif if (proto) then write(10,*) ' full step ' write(10,*) ' working set ' write(10,fmt='(10(i4,1x))') (a(i),i=1,iq) write(10,*) ' vector u' write(10,fmt='(5(d14.7,1x))') (u(i),i=1,iq) endif goto 50 endif iai(l)=l call delc(l,n,m,p) if (proto) then write(10,*) ' partial step ' write(10,*) ' working set ' write(10,fmt='(10(i4,1x))') (a(i),i=1,iq) write(10,*) ' vector u ' write(10,fmt='(5(d14.7,1x))') (u(i),i=1,iq) endif su=0.d0 do k=1,n su=su+ci(k,ip)*x(k) enddo s(ip)=su+ci0(ip) goto 100 end c\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ subroutine zup(z,n) c update z include 'dualcom.inc' integer n double precision z(nv) c*** local integer i,j double precision su cc d = j(trans) *np do i=1,n su=0.d0 do j=1,n su=su+xj(j,i)*np(j) enddo d(i)=su enddo cc bestimmung von z do i=1,n z(i)=0.d0 do j=iq+1,n z(i)=z(i)+xj(i,j)*d(j) enddo enddo return end c\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ subroutine rup(rv) c update r include 'dualcom.inc' c*** local integer i,j double precision s double precision rv(mx) do i=iq,1,-1 s=0.d0 do j=i+1,iq s=s+r(i,j)*rv(j) enddo rv(i)=(d(i)-s)/r(i,i) enddo return end c\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ subroutine delc(l,n,m,p) include 'dualcom.inc' integer l,n,m,p c*** local integer i,j,k,qq,me double precision cc,ss,h,c1,s1,xny,t1,t2,dsq1 character*40 line1,line2 external dsq1 data line1/' matrix r '/ data line2/' matrix j '/ if (proto) then write(10,fmt='('' delete restriction '',i3)') l endif me=p do i=me+1,iq if(a(i) .eq. l) then qq=i goto 10 endif enddo 10 continue do i=qq,iq-1 a(i)=a(i+1) u(i)=u(i+1) do j=1,n r(j,i)=r(j,i+1) enddo enddo 20 continue a(iq)=a(iq+1) u(iq)=u(iq+1) a(iq+1)=0 u(iq+1)=0.d0 do j=1,iq r(j,iq)=0.d0 enddo iq=iq-1 if(iq .eq. 0) goto 100 do j=qq,iq cc=r(j,j) ss=r(j+1,j) h=dsq1(cc,ss) if(h .eq. 0.d0) goto 90 c1=cc/h s1=ss/h r(j+1,j)=0.d0 if ( c1 .lt. 0.d0 ) then r(j,j)=-h c1=-c1 s1=-s1 else r(j,j)=h endif xny=s1/(1.d0+c1) do k=j+1,iq t1=r(j,k) t2=r(j+1,k) r(j,k)=t1*c1+t2*s1 r(j+1,k)=xny*(t1+r(j,k))-t2 enddo do k=1,n t1=xj(k,j) t2=xj(k,j+1) xj(k,j)=t1*c1+t2*s1 xj(k,j+1)=xny*(xj(k,j)+t1)-t2 enddo 90 continue enddo 100 continue if (proto) then if ( iq .ne. 0 ) f call matdru(r,mx,mx,iq,iq,line1,10,.false.,.true.) call matdru(xj,nv,nv,n,iq,line2,10,.true.,.false.) endif return end c\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ subroutine addc(ifail,n) include 'dualcom.inc' integer ifail,n c*** local integer i,j,k double precision cc,ss,h,c1,s1,t1,t2,xny,dsq1 character*40 line1,line2 external dsq1 intrinsic max,dabs data line1/' matrix r '/ data line2/' matrix j '/ if (proto) then write(10,*) ' add new restriction ' endif ifail=0 do j=n,iq+2,-1 cc=d(j-1) ss=d(j) h=dsq1(cc,ss) if(h .eq. 0.d0) goto 20 d(j)=0.d0 s1=ss/h c1=cc/h if ( c1 .lt. 0.d0 ) then c1=-c1 s1=-s1 d(j-1)=-h else d(j-1)=h endif xny=s1/(1.d0+c1) do k=1,n t1=xj(k,j-1) t2=xj(k,j) xj(k,j-1)=t1*c1+t2*s1 xj(k,j)=xny*(t1+xj(k,j-1))-t2 enddo 20 continue enddo iq=iq+1 do i=1,iq r(i,iq)=d(i) enddo if ( dabs(d(iq)) .le. epsmac*rnorm ) then if (proto) then write(10,*) ' problem degenerate ' endif ifail=1 return endif rnorm=max(rnorm,dabs(d(iq))) if (proto) then if ( iq .ne. 0 ) f call matdru(r,mx,mx,iq,iq,line1,10,.false.,.true.) call matdru(xj,nv,nv,n,iq,line2,10,.true.,.false.) endif return end c***************************************************************** c subprogram for structured output of a matrix c***************************************************************** subroutine matdru(a,me,ne,ma,na,head,lognum,fix,triang) implicit none integer anz,i,j,k,me,ne,ma,na,lognum,ju,jo double precision a(me,ne) logical fix,triang character*40 head character*131 prline character*15 blank,aij(4) character*11 itext data blank/' '/ if ( ma .gt. 99 .or. na .gt. 99 ) then write(lognum,fmt='(/1x,a40)') head write(lognum,*) 'matdru mat too large' return endif write(lognum,fmt='(/1x,a40)') head anz=4 jo=0 do while ( jo .lt. na ) ju=jo+1 jo=min(ju+anz-1,na) write(lognum,fmt='(''row/column'',4(6x,i3,6x))') (j,j=ju,jo) do i=1,ma if ( .not. triang ) then if ( fix ) then write(lognum,fmt='(1x,3x,i4,3x,4(g14.7,1x))') f i,(a(i,j),j=ju,jo) else write(lognum,fmt='(1x,3x,i4,3x,4(d14.7,1x))') f i,(a(i,j),j=ju,jo) endif else write(itext,fmt='(1x,3x,i4,3x)') i do k=1,4 aij(k)=blank enddo k=0 do j=ju,jo k=k+1 if ( j .ge. i ) then if ( fix ) then write(aij(k),fmt='(g14.7,1x)') a(i,j) else write(aij(k),fmt='(d14.7,1x)') a(i,j) endif endif enddo prline=itext//aij(1)//aij(2)//aij(3)//aij(4) write(lognum,fmt='(a80)') prline endif enddo enddo return end c******************************************************************* c scalar product of two vectors or parts of vectors c******************************************************************* double precision function skal1(i,j,a,b) implicit none integer i,j,k double precision a(*),b(*) double precision s if ( i .gt. j ) then skal1=0.d0 return else s=0.d0 do k=i,j s=s+a(k)*b(k) enddo skal1=s return endif end c******************************************************************* subroutine leftel(a,b,y,yl,n) implicit none include 'dualpar.inc' integer n,i,j double precision a(n,n),b(nv),y(nv),yl,h c leftel assumes that the cholesky-factor of a c a=r(transpose)*r is stored in the upper half of a. c b is a right hand side. leftel solves c y(transpose)*r = b(transpose) c yl=norm(y)**2 yl=0.d0 do i=1,n h=b(i) do j=1,i-1 h = h - a(j,i)*y(j) enddo h=h/a(i,i) y(i)=h yl = h**2 + yl enddo return end c********************************************************** double precision function dsq1(a,b) c c computes sqrt(a**2+b**2) numerically safe c implicit none double precision a,b,a1,b1 intrinsic dabs,dsqrt a1=dabs(a) b1=dabs(b) if ( a1 .gt. b1 ) then dsq1=a1*dsqrt(1.d0+(b1/a1)**2) else if ( b1 .gt. a1 ) then dsq1=b1*dsqrt(1.d0+(a1/b1)**2) else dsq1=a1*dsqrt(2.d0) endif endif return end c******************************************************************* c cholesky decomposition of a symmetric positive matrix a c the strict lower triangle remains unaffected c the upper triangle including the diagonal holds the cholesky- c factor. initially the strict upper triangle may be undefined c******************************************************************* subroutine chol(a,n,r,ifail) implicit none include 'dualpar.inc' integer n,ifail double precision a(n,n),r(n,n),h,s c computes the cholesky-factor r of a=r(transp)*r c and stores it in the upper triangle of r c a and r may be identical. the strict lower triangle c of a remains unchanged anyway. c ifail .ne. 0 if the decomposition does'nt exist, c =0 otherwise integer i,j,k ifail=0 do i=1,n h=a(i,i) do j=1,i-1 c h=h-r(j,i)**2 h=h-(r(j,i)*r(j,i)) enddo if ( h .le. 0.d0 ) then write(14,*)'hhhhh',h ifail=1 return endif h=dsqrt(h) r(i,i)=h h=1.d0/h do k=i+1,n s=0.d0 do j=1,i-1 s=s+r(j,i)*r(j,k) enddo s=(a(k,i)-s)*h c s=(a(k,i)-s)/h r(i,k)=s enddo enddo return end c******************************************************** subroutine cholso(a,n,b,x) c copyright = spellucci implicit none include 'dualpar.inc' integer n,i,j double precision a(n,n),b(n),c(nv),x(n) double precision s c solves the linear system a*x=b, a symmetric positive definite c it is assumed that the cholesky-factor r of c a = r(transpose)*r c is stored in the upper triangle of a by chol. do i=1,n s=0.d0 do j=1,i-1 s=s+a(j,i)*c(j) enddo c(i)=(b(i)-s)/a(i,i) enddo do i=n,1,-1 s=0.d0 do j=i+1,n s=s+a(i,j)*x(j) enddo x(i)=(c(i)-s)/a(i,i) enddo return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Solving the quadtatic problem by goldfarb idnani dual method with elleminig the c equality constraint ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine resoudre2(M,C,U,dvaf,nu,nm,bor) c$debug include"stokcom.fd" integer nu,nm,bor,nu1,arr real*8 m(nx1,nx1),c(nx1),u(nx1),dmx,BQ(nu-1,nu-1),AQ(nu-1,nu), + CQ(nu-1),SOL(nu-1),CEO(0),CIO(nu),MCE(nu-1,0),f(nx1),dvaf, + dcpvv external dcpvv c external DLSLSF c call RSEL(M,d,d1,nx1,nm) bor=1 nu1=nu-1 do i=1,nu1 dmx=M(i,nu)-M(nu,nu) do j=1,nu1 BQ(i,j)=M(i,j)-dmx-M(nu,j) AQ(i,j)=0.d0 enddo dmx=M(i,nu)-(M(nu,nu)+c(nu)) CQ(i)=C(i)+dmx CIO(i)=0.d0 c MCE(i,1)=1.d0 AQ(i,i)=1.d0 enddo do i=1,nu1 AQ(i,nu)=-1.d0 enddo CIO(nu)=1.d0 c CEO(1)=1.d0 call driver(BQ,AQ, MCE, CQ, CIO, CEO, SOL,nu1,nu,0) c CALL DQPROG (nu,nu+1,1,AQ,nu+1, BSQ,CQ,BQ,nu, DIAG, c + SOL, NACT, IACT, ALAMDA) U(nu)=1.d0 do i=1,nu1 U(i)=sol(i) U(nu)=U(nu)-U(i) write( 14,*)'sol(',i,')=',sol(i) enddo write(14,*)'Sol(',nu,')=',U(nu) c call securute(M,C,U,nu,arr) c do i=1,nu c write( 14,*)'u(',i,')=',u(i) c enddo call DCPMV(m,u,f,nm,nx1,nx1) dvaf=(DCPVV(f,u,nm,nx1)/2.D0)+DCPVV(c,u,nm,nx1) c write(*,*)'arrrr',dvaf return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c stok of active constraints cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc Subroutine CACTIVE (GRADC,B,SD,BXL,BXU,VAC) include"stokcom.fd" Parameter(deps=1.d-8) integer VAC(m1+nx),k real*8 GRADC(mc,nx), B(mc),zs,SD(nx),BXL(nx+m1),BXU(nx+m1) zs=0.d0 k=0 do i=m2+1,mc do j=1,nx zs=zs+GRADC(i,j)*SD(j) enddo if ( DABS(zs-B(i)).le.deps) then VAC(i+nx-m2)=1 k=k+1 else VAC(i+nx-m2)=0 endif enddo do i=1,nx if ((DABS(sd(i)-BXL(i)).LE.deps).or.(DABS(sd(i)-BXU(i)).LE.deps + )) then VAC(i)=1 k=k+1 else VAC(i)=0 endif enddo write(14,*)'nombre des contraintes actives est',k return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc Subroutine simdual(A,As,MB,MBs,Cc,SM,SMM,Iposv,L1,lp1,L3,LL1,nll1, + nl1,nl3,LC,Lp2,Nvar,nlc,Icase,itn,m1,m2,m3,m,n,npe,isdp) $debug Integer m,n,Iposv(m),As(m+1,n),MBs(n+1,n), + L1(m),nl1,L2(m),L3(n+m),nl3,nl2,icase,Xp(n),Cp(m), + md1,md2,md3,Mbsd(n+1,n),Asd(n+1,m),Iposvd(n),Ld1(n), + nld1,nld3,itn,lp1(m),LL1(m),nll1,LC(m),nlc,Lp2(m),Lpp(m), + NVar(n),nlpp,npe Real*8 A(m,n),MB(m+1,m),Cc(n+m),dobj,SM(m+1),Ad(n,m),Ccd(m+n), +SMM(n+1), MBd(m+1,m),rrrs,deps deps=1.d-12 dobj=0.d0 ire=0 if(itn.eq.1) then if (isdp.eq.1) then call simdp(A,SM,SMM,Cc,MB,MBs,As,Iposv,L1,lp1,L3,LL1,nll1, + nl1,nl3,LC,Lp2,Nvar,nlc,Icase,itn,m1,m2,m3,m,n) else c stop call simprev(A,SM,Cc,MB,MBs,As,Iposv,L1,L3,nl1,nl3,Icase,itn, + m1,m2,m3,m,n,npe) endif else Do i=1,n k=1 do j=1,m if(DABS(A(j,i)).GE.deps) then As(k+1,i)=j k=k+1 endif enddo As(1,i)=k-1 enddo do i=1,m c SM(i)=0.d0 Iposv(i)=n+i enddo do i=1,n Xp(i)=1 Nvar(i)=0 enddo if(isdp.eq.1) then call cvd(Xp,Cp,A,Ad,Sm,SMM,Cc,Ccd,m1,m2,m3,m,md1,md2,md3,n) c nl1=m do i=1,n do j=1,n MB(i,j)=0.d0 enddo MBs(1,i)=1 MB(i,i)=1.d0 L1(i)=i Lp1(i)=0 Lp2(i)=0 enddo nl1=n nll1=0 do i=n+1,m nll1=nll1+1 LL1(nll1)=i Lp1(i)=0 Lp2(i)=0 enddo do i=1,m1+m2 Lpp(i)=i enddo nlpp=m1+m2 c do i=1,m c write(6,*)'Smd(',i,')=',Ccd(i) c enddo c stop call simprev(Ad,SMM,Ccd,MBd,MBsd,Asd,Iposvd,Ld1,L3,nld1, +nld3,Icase,itn,md1,md2,md3,n,m,npe) c stop nl3=0 nlc=0 c nll1=0 do i=1,nld1 kkk=ld1(i) c kkkk=Cp(kkk) ii=iposvd(kkk) if(ii.GT.m) then if(ii.LE.m+m1+m2) then nll1=nll1+1 LL1(nll1)=ii-m endif nl3=nl3+1 c L3(nl3)=kkkk jj=ii-m kkkk=Cp(jj) L3(nl3)=kkkk do jjj=1,n if(L1(jjj).eq.jj) goto 87 enddo 87 continue do j=jjj,nl1-1 L1(j)=L1(j+1) enddo nl1=nl1-1 else nlc=nlc+1 Lc(nlc)=ii Lp2(ii)=nlc if(ii.LE.m1+m2) then do jjj=1,m1+m2 if(Lpp(jjj).eq.ii) goto 78 enddo 78 continue do j=jjj,nlpp-1 Lpp(j)=Lpp(j+1) enddo nlpp=nlpp-1 endif endif enddo do i=1,nl1 kkk=L1(i) kkkkk=nll1+i LL1(kkkkk)=L1(i) c write(6,*)'kkkkkk',i,kkk Lp1(kkk)=i rrrs=0.d0 do j=1,Mbsd(1,kkk) k=Mbsd(j+1,kkk) kk=iposvd(k) if (kk.LE.m) then rrrs=rrrs+Ccd(kk)*MBd(k,kkk) endif enddo if(((kkk.GT.md1).and.(kkk.Le.md1+md2)).or.((kkk.GT.md1+md2).and. + (Cc(kkk).Le.0.d0))) then rrrs=-rrrs endif kkkk=Cp(kkk) SMM(i)=-rrrs iposv(kkk)=kkkk nl3=nl3+1 L3(nl3)=Lc(i)+n Nvar(kkkk)=1 enddo do i=1,nll1 k=LL1(i) kk=Lpp(i) iposv(k)=kk+n enddo c if((i.GT.md1+md2).and.(Cc(i).Le.0.d0)) then c SM(i)=-SM(i) c endif c iposv(ii)=kkkk do i=1,nld1 k=Ld1(i) c if((k.GT.md1+md2).and.(Cc(k).Le.0.d0)) then if(((k.GT.md1).and.(k.Le.md1+md2)).or.((k.GT.md1+md2).and. + (Cc(k).Le.0.d0))) then do j=1,nld1 kk=Ld1(j) MBd(kk,k)=-Mbd(kk,k) enddo endif enddo do i=1,nld1 k=Ld1(i) kk=iposvd(k) if(kk.le.m1) then do j=1,nld1 kkk=Ld1(j) MBd(k,kkk)=-MBd(k,kkk) enddo endif enddo jj=0 do j=1,nld1 k=Ld1(j) kk=iposvd(k) if(kk.LE.m) then jj=jj+1 ii=1 do i=1,nl1 kkk=L1(i) MB(i,jj)=MBd(k,kkk) if(DABS(MB(i,jj)).GE.deps) then MBs(ii+1,jj)=i ii=ii+1 endif enddo MBs(1,jj)=ii-1 endif enddo c stop dobj=0.d0 do i=1,nl1 ii=L1(i) k=iposv(ii) if(k.LE.n) then dobj=dobj+SMM(i)*Cc(k) endif enddo write(6,*)'nld1',nld1,nll1,nl1 c endif c enddo c stop else call simprev(A,SM,Cc,MB,MBs,As,Iposv,L1,L3,nl1,nl3,Icase,itn, + m1,m2,m3,m,n,npe) endif endif return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine simprev(A,SM,Cc,MB,MBs,As,Iposv,L1,L3,nl1,nl3, +Icase,itn,m1,m2,m3,m,n,npe) $debug Integer m,n,Iposv(m),As(m+1,n),MBs(m+1,m), + L1(m),nl1,L2(m),L3(n+m),nl2,k,kk,kkk,kkkk,icase,born, + Vaps(m+1),pp,p,q,qq,ich,MBsq(m+1),LL1(m2+m3),nll1,itn,nl3,npe Real*8 A(m,n), MB(m+1,m),Cr(n+m2),Cc(n+m),rrrs,dobj,SM(m+1),Binf, + Vap(m),deps c Phase une ntest=0 ich=0 deps=1.d-10 born=0 if(itn.eq.1) goto 30 Do i=1,m do j=1,m MB(i,j)=0.d0 enddo MB(i,i)=1.d0 MBs(1,i)=1 MBs(2,i)=i enddo do i=1,m1 Iposv(i)=n+i L2(i)=1 L1(i)=i enddo nl1=m nll1=m2+m3 nl2=m nl3=n+m2 do i=m1+1,m Iposv(i)=n+i LL1(i-m1)=i L2(i)=1 L1(i)=i enddo Do i=1,n k=1 do j=1,m if(DABS(A(j,i)).GE.deps) then As(k+1,i)=j k=k+1 endif enddo As(1,i)=k-1 L3(i)=i enddo Do i=1,m2 L3(n+i)=n+m+i enddo 44 continue c write(6,*)'debut1' ntest=ntest+1 call Calculcr(A,As,MB,Mbs,iposv,Cc,Cr,L1,LL1,L2,L3, + nl1,nll1,nl2,nl3,m1,m2,m3,m,n,ich) Call SIMP1(Cr,L3,nl3,Binf,pp,p,m1,m2,m3,m,n) dobj=0.d0 do i=1,m k=iposv(i) if((k.GT.n+m1).and.(k.LE.n+m)) then dobj=dobj+SM(i) endif enddo c if (ntest.eq.4) stop c if (ntest.eq.3) stop if ((Binf.GE.-deps).and.(dobj.GE.deps)) then Write(6,*)'le probleme est non admissible' stop return elseif(Binf.GE.-deps) then c Write(6,*)'fin phase une',nll1 c stop do i=1,nll1 k=LL1(i) q=k if(k.LE.m1+m2) then p=n+m+k-m1 do j=1,nl3 if(L3(j).eq.p) goto 16 enddo 16 continue rrrs=1.d0 pp=j else do j=1,nl3 kk=L3(j) if (kk.GT.n) then if (kk.GT.n+m) then kkk=kk-n-m+m1 else kkk=kk-n endif rrrs=MB(k,kkk) else kkk=As(1,kk) rrrs=0.d0 do kkkk=1,kkk kkkkk=As(kkkk+1,kk) rrrs=rrrs+MB(k,kkkkk)*A(kkkkk,kk) enddo endif if(DABS(rrrs).GT.deps) then p=kk pp=j if (p.LE.n) then do ii=1,m Vap(ii)=A(ii,p) enddo do ii=1,As(1,p) Vaps(ii+1)=As(ii+1,p) enddo Vaps(1)=As(1,p) endif goto 17 endif enddo 17 continue endif if(DABS(rrrs).LE.deps) then do jj=1,nl1 if(L1(jj).eq.k) goto 22 enddo 22 continue nl1=nl1-1 do j=jj,nl1 L1(j)=L1(j+1) enddo do j=1,n if (DABS(A(k,j)).GE.Deps) then do jj=1,As(1,j) if(As(jj+1,j).eq.k) goto 9 enddo 9 continue As(1,j)=As(1,j)-1 do jjj=jj,As(1,j) As(jjj+1,j)=As(jjj+2,j) enddo endif enddo do j=1,m if (DABS(MB(k,j)).GE.Deps) then do jj=1,MBs(1,j) if(MBs(jj+1,j).eq.k) goto 10 enddo 10 continue MBs(1,j)=MBs(1,j)-1 do jjj=jj,MBs(1,j) MBs(jjj+1,j)=MBs(jjj+2,j) enddo endif enddo else call CalculCA(MB,MBS,Vap,Vaps,p,nl1,L1,ich,m,m1,m2,m3,n) kk=1 do ii=1,m if (DABS(MB(q,ii)).GE.deps) then MBsq(kk+1)=ii kk=kk+1 endif enddo MBsq(1)=kk-1 L2(q)=0 call MJMB(Vap,Vaps,q,MB,MBs,MBsq,Sm,L1,m,m1,m2,m3,n, + nl1,ich) npe=npe+1 nl3=nl3-1 do ii=pp,nl3 L3(ii)=L3(ii+1) enddo endif enddo do i=1,nl3 if (L3(i).GT.n+m) then L3(i)=L3(i)+m1-m endif enddo do i=1,m if (iposv(i).GT.n+m) then iposv(i)=iposv(i)+m1-m endif enddo goto 30 endif if (p.LE.n) then do i=1,As(1,p) k=As(i+1,p) Vap(k)=A(k,p) Vaps(i+1)=k enddo Vaps(1)=As(1,p) endif call CalculCA(MB,MBS,Vap,Vaps,p,nl1,L1,ich,m,m1,m2,m3,n) call SIMP2(A,As,MB,MBs,Sm,Vap,qq,q,L1,nl1,born,m1,m2,m3,m,n) k=1 do i=1,m if (DABS(MB(q,i)).GE.deps) then MBsq(k+1)=i k=k+1 endif enddo MBsq(1)=k-1 L2(q)=0 call MJMB(Vap,Vaps,q,MB,MBs,MBsq,Sm,L1,m,m1,m2,m3,n,nl1,ich) npe=npe+1 do i=1,m do j=1,MBS(1,i) k=MBs(j+1,i) c write(6,*)'MBs',k,MB(k,i) enddo enddo do i=1,vaps(1) k=vaps(i+1) c write(6,*)'Vap',k,vap(k) enddo c write(6,*)'p,pp,q,qq',p,pp,q,qq if ((iposv(q).GT.n+m1).and.(iposv(q).le.n+m)) then do i=1,nll1 if(LL1(i).eq.qq) goto 71 enddo 71 continue nll1=nll1-1 do j=i,nll1 LL1(j)=LL1(j+1) enddo nl3=nl3-1 do i=pp,nl3 L3(i)=L3(i+1) enddo iposv(q)=p else c Write(14,*)'p,q,pp,qq,iposv(q),L3(pp)',p,q,pp,qq,iposv(q),L3(pp) k=iposv(q) iposv(q)=p L3(pp)=k endif c do i=1,m c write(6,*)'iposv(',i,')=',iposv(i) c enddo goto 44 30 continue ntest=ntest+1 c if (ntest.eq.6) stop ich=1 c write(6,*)'nnnnnl1',nl1 call Calculcr(A,As,MB,Mbs,iposv,Cc,Cr,L1,LL1,L2,L3, + nl1,nll1,nl2,nl3,m1,m2,m3,m,n,ich) c if (ntest.eq.6) stop c write(6,*)'matrice inverse' c do i=1,m c do j=1, mbs(1,i) c write(6,*)'MB(',MBS(j+1,i),',',i,')=',MB(MBs(j+1,i),i) c enddo c enddo Call SIMP1(Cr,L3,nl3,Binf,pp,p,m1,m2,m3,m,n) if (Binf.GE.-deps) then icase=0 write(6,*)' probleme borne' return endif c write(6,*)'asssss',p,As(1,p) if (p.LE.n) then do i=1,As(1,p) k=As(i+1,p) Vap(k)=A(k,p) c write(14,*)'p,vap',p,k,nl1,vap(k) Vaps(i+1)=k enddo Vaps(1)=As(1,p) endif c write(14,*)'mbs',MBs(2,2),MBs(2,3) c write(14,*)'lllllllll' c write(6,*)'MBBB',MB(1,1),MB(2,1),MB(3,1),MBS(1,1) call CalculCA(MB,MBS,Vap,Vaps,p,nl1,L1,ich,m,m1,m2,m3,n) c write(14,*)'vvvvvvvvvvvvv',vap c write(6,*)'vvvap',vaps(1),vap(1),vap(2),vap(3) call SIMP2(A,As,MB,MBs,Sm,Vap,qq,q,L1,nl1,born,m1,m2,m3,m,n) if (born.eq.1) then write(6,*)'the linear problem is unbounded, add for example + the constraint xi<100000, i=1,nx ' write(14,*)'the linear problem is unbounded, add for example + the constraint xi<100000, i=1,nx ' stop return endif k=1 do i=1,m if (DABS(MB(q,i)).GE.deps) then MBsq(k+1)=i k=k+1 endif enddo MBsq(1)=k-1 L2(q)=0 call MJMB(Vap,Vaps,q,MB,MBs,MBsq,Sm,L1,m,m1,m2,m3,n,nl1,ich) npe=npe+1 c Write(14,*)'p,q,pp,qq,iposv(q),L3(pp)',p,q,pp,qq,iposv(q),L3(pp) k=iposv(q) iposv(q)=p L3(pp)=k dobj=0.d0 do i=1,nl1 k=L1(i) kk=iposv(k) if(kk.Le.n) then dobj=dobj+cc(kk)*SM(k) endif enddo goto 30 end Subroutine Calculcr(A,As,MB,Mbs,iposv,Cc,Cr,L1,LL1,L2, + L3,nl1,nll1,nl2,nl3,m1,m2,m3,m,n,ich) Integer m1,m2,m3,m,n,nl2,iposv(m),L1(m), + L2(m),L3(n+m),nl3,nl1,Lams(m+1),As(m+1,n),Mbs(m+1,m), + k,kk,LL1(m2+m3),nll1 Real*8 A(m,n),Mb(m+1,m),Cr(n+m2),rrrs,Cc(n+m),Lam(m),deps c probleme artificiel c L3 represente les variables hors base c nl3 nombre de variables hors base c L1 represente les variables artificielle de base c nl1 represente le nbe de variables artificielles courant do i=1,m lam(i)=0.d0 enddo deps=1.d-10 if (ich.eq.0) then k=1 do i=1,m if (L2(i).eq.1) then c Write(6,*)'lllllllll',i if ((iposv(i).GT.n+m1 ).and.(iposv(i).LE.n+m)) then Lam(i)=MB(i,i) if(DABS(Lam(i)).GE.deps) then lams(k+1)=i k=k+1 endif endif else rrrs=0.d0 do j=1,nll1 kk=LL1(j) rrrs=rrrs+MB(kk,i) enddo Lam(i)=rrrs if(DABS(Lam(i)).GE.deps) then lams(k+1)=i k=k+1 endif endif enddo lams(1)=k-1 do i=1,nl3 k=L3(i) rrrs=0.d0 if(k.GT.n) then if (k.Le.n+m) then rrrs=-lam(k-n) else rrrs=lam(k-n-m+m1) endif Cr(i)=rrrs else c write(6,*)'kk',k if(lams(1).LE.As(1,k)) then kk=Lams(1) do j=1,kk kkk=lams(j+1) rrrs=rrrs+lam(kkk)*A(kkk,k) enddo else kk=As(1,k) do j=1,kk kkk=As(j+1,k) rrrs=rrrs+lam(kkk)*A(kkk,k) enddo endif Cr(i)=-rrrs endif enddo c Phase deux else kkkk=1 do i=1,nl1 rrrs=0.d0 k=L1(i) kk=MBs(1,k) do j=1,kk kkk=MBs(j+1,k) c write(6,*)'kkkkkkkk',j,kk,kkk c if(iposv(kkk).Le.n) then c write(6,*)'kkk,k,iposv(kkk)', kkk,k,iposv(kkk) rrrs=rrrs+Cc(iposv(kkk))*MB(kkk,k) c endif enddo Lam(k)=rrrs if(DABS(Lam(k)).GE.deps) then lams(kkkk+1)=k kkkk=kkkk+1 endif enddo lams(1)=kkkk-1 do i=1,nl3 k=L3(i) rrrs=0.d0 if(k.GT.n) then if (k.Le.n+m1) then rrrs=-lam(k-n) else rrrs=lam(k-n) endif Cr(i)=Cc(k)+rrrs else if(lams(1).LE.As(1,k)) then kk=Lams(1) do j=1,kk kkk=lams(j+1) rrrs=rrrs+lam(kkk)*A(kkk,k) enddo else kk=As(1,k) do j=1,kk kkk=As(j+1,k) rrrs=rrrs+lam(kkk)*A(kkk,k) enddo endif Cr(i)=Cc(k)-rrrs endif enddo endif return end Subroutine CalculCA(MB,MBS,Vap,Vaps,p,nl1,L1,ich,m,m1,m2,m3,n) Integer MBs(m+1,m),Vaps(m+1),L1(m),nl1,m,m1,m2,m3,n,p,k,kk Real*8 MB(m+1,m),Vap(m),Vap1(m),deps do i=1,m vap1(i)=0.d0 enddo deps=1.d-10 If (ich.eq.0) then if (p.GT.n) then if(p.LE.n+m1) then k=p-n kk=MBs(1,k) do i=1,kk kkk=MBs(i+1,k) Vap1(kkk)=MB(kkk,k) Vaps(i+1)=kkk enddo Vaps(1)=kk else c write(6,*)'p,n,m,m1,m2',p,n,m,m1,m2 k=p-n-m+m1 c write(6,*)'kkkkk',k kk=MBs(1,k) do i=1,kk kkk=MBs(i+1,k) Vap1(kkk)=-MB(kkk,k) Vaps(i+1)=kkk enddo Vaps(1)=kk endif else kk=Vaps(1) do i=1,kk k=Vaps(i+1) kkk=MBs(1,k) do j=1,kkk kkkk=MBs(j+1,k) Vap1(kkkk)=Vap1(kkkk)+Vap(k)*MB(kkkk,k) enddo enddo kk=1 do i=1,m Vap(i)=Vap1(i) c write(6,*)'***p,vap',p,vap(i) if(Dabs(Vap(i)).Ge.deps) then VAps(kk+1)=i kk=kk+1 endif enddo Vaps(1)=kk-1 endif else if (p.GT.n) then do i=1,m vap(i)=0 enddo if(p.LE.n+m1) then k=p-n kk=MBs(1,k) do i=1,kk kkk=MBs(i+1,k) Vap1(kkk)=MB(kkk,k) Vaps(i+1)=kkk enddo Vaps(1)=kk c write(14,*)'kkkkkkkkkkkk',p,k,MB(58,k),Vap(58) else k=p-n kk=MBs(1,k) do i=1,kk kkk=MBs(i+1,k) Vap1(kkk)=-MB(kkk,k) Vaps(i+1)=kkk enddo Vaps(1)=kk endif else k=Vaps(1) c write(14,*)'7777777777k',k do i=1,k kk=Vaps(i+1) c write(6,*)'888888kk',kk kkk=MBs(1,kk) c write(6,*)'999999999kkk',kkk do j=1,kkk kkkk=MBs(j+1,kk) Vap1(kkkk)=Vap1(kkkk)+Vap(kk)*MB(kkkk,kk) enddo enddo kk=1 do i=1,nl1 k=L1(i) Vap(k)=Vap1(k) if(Dabs(Vap(k)).Ge.deps) then VAps(kk+1)=k kk=kk+1 endif enddo Vaps(1)=kk-1 endif endif do i=1,m Vap(i)=vap1(i) enddo return end Subroutine MJMB(Vap,Vaps,q,MB,MBs,MBsq,Sm,L1,m,m1,m2,m3,n,nl1,ich) Integer L1(m),MBs(m+1,m),q,Vaps(m+1),m1,m2,m3,m,n,nl1,MBsq(m+1), + kkk,k,kk Real*8 Vap(m),MB(m+1,m),Sm(m+1),deps,dpiv,rrrs deps=1.d-10 k=MBsq(1) do i=1,Vaps(1) kk=Vaps(i+1) c write(6,*)'vap',q dpiv=1.d0/Vap(q) if((kk.LT.q).or.(kk.GT.q)) then rrrs=-Vap(kk)*dpiv do j=1,k kkk=MBsq(j+1) MB(kk,kkk)=MB(kk,kkk)+rrrs*MB(q,kkk) enddo Sm(kk)=Sm(kk)+rrrs*Sm(q) endif enddo rrrs=dpiv do j=1,k kkk=MBsq(j+1) MB(q,kkk)=rrrs*MB(q,kkk) enddo Sm(q)=rrrs*Sm(q) if(ich.eq.0) then do i=1,k kkk=1 kk=MBsq(i+1) do j=1,m if(DABS(MB(j,kk)).GE.deps) then MBs(kkk+1,kk)=j kkk=kkk+1 endif enddo MBs(1,kk)=kkk-1 enddo else do i=1,k kkk=1 kk=MBsq(i+1) do j=1,nl1 kkkk=L1(j) if(DABS(MB(kkkk,kk)).GE.deps) then MBs(kkk+1,kk)=kkkk kkk=kkk+1 endif enddo MBs(1,kk)=kkk-1 enddo endif return end c ****************************************************** Subroutine SIMP1(Cr,L3,nl3,Binf,pp,p,m1,m2,m3,m,n) Integer p,pp,m1,m2,m3,n,L3(n+m),nl3,k Real*8 Cr(n+m2),Binf,test,deps Binf=Cr(1) p=L3(1) pp=1 do 121 k=2,nl3 TEST=Cr(k)-Binf if(test.Le.0.d0)then Binf=Cr(k) p=L3(k) pp=k endif 121 continue c if(p.GT.n+m+m2) stop return end c **************************************************************** c Locate a pivot element taking degeneracy into account Subroutine SIMP2(A,As,MB,MBs,Sm,Vap,qq,q,L1,nl1,born,m1,m2,m3,m,n) Integer L1(m),nl1,q,m1,m2,m3,n,born,MBs(m+1,m),As(m+1,n),qq real*8 Sm(m+1),Vap(m),Q1,Q0,Q2,QP,MB(m+1,m),A(m,n) deps=1.d-10 do 11 i=1,Nl1 k=L1(i) c write(6,*)'Vap',k,vap(k) c write(24,*)' A(',L2(i)+1,',',KP+1,')=',A(L2(i)+1,KP+1) if(Vap(k).GT.deps) GOTO 2 11 continue born=1 return 2 Q1=Sm(k)/Vap(k) q=k qq=i if(i+1.GT.NL1)return Do 13 ii=i+1,NL1 k=L1(ii) if(Vap(k).GT.deps) then Q2=Sm(k)/Vap(k) if (Q2.LE.Q1) then q=k qq=ii Q1=Q2 c We have a degenery elseif (Q2.EQ.Q1) then do 12 j=1,N rrrs1=0.d0 rrrs2=0.d0 do jj=1,As(1,j) jjj=As(jj+1,j) rrrs1= rrrs1+MB(k,jjj)*A(jjj,j) rrrs2=rrrs2+MB(q,jjj)*A(jjj,j) enddo QP=rrrs1/Vap(k) Q0=rrrs2/Vap(q) if(Q0.NE.QP) GOTO 6 12 continue 6 if(Q0.GT.QP) then q=k qq=ii endif endif endif 13 continue Return end c **************************************************************** ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine simdp(A,SM,SMM,Cc,MB,MBs,As,Iposv,L1,lp1,L3,LL1,nll1, + nl1,nl3,LC,Lp2,Nvar,nlc,Icase,itn,m1,m2,m3,m,n) $debug Integer m,n,Iposv(m),As(m+1,n),MBs(n+1,n),lp1(m), + L1(m),nl1,L2(n),L3(n+m),nl2,k,kk,kkk,icase,born, + Vaps(n+1),pp,p,q,qq,ich,MBsq(n+1),LL1(m),nll1,itn, + Avecs(n+1),LC(m),nlc,Lp2(m),nit,Nvar(n) Real*8 A(m,n), MB(n+1,n),Cr(n+m2),Cc(n+m),rrrs,SM(m+1),Binf, + Vap(n+1),deps,Avec(n),Vapp(m),SMM(n+1),SMMM(m),rrrs1 c Phase une ntest=0 deps=1.d-14 born=0 nit=0 30 continue ich=1 c write(6,*)'nnnnnl1',nl1 call Calculdcr(A,As,MB,Mbs,iposv,Cc,Cr,L1,L2, + L3,Lc,LL1,nll1,nlc,nl1,nl2,nl3,m1,m2,m3,m,n) c write(6,*)'coefficient cout reduit' c do i=1,n+m2 c write(6,*)'Cr(',L3(i),')=',Cr(i) c enddo nit=nit+1 Call SIMPd1(Cr,L3,nl3,Binf,pp,p,m1,m2,m3,m,n) rrrs=0.d0 do i=1,nl1 k=L1(i) rrrs=rrrs+Cc(iposv(k))*SMM(i) enddo write(14,*)'p,pp,,qq,Bi,rrr',p,pp,q,qq,Binf,rrrs if (Binf.GE.-1.d-12) then icase=0 write(6,*)' probleme borne' c do i=1,m c k=iposv(i) c if(k.le.n) then c do ii=1,nl1 c rrrs=0.d0 c rrrs1=0.d0 c do j=1,nl1 c kk=Lc(j) c rrrs=rrrs+MB(ii,j)*A(kk,k) c rrrs1=rrrs1+MB(ii,j)*SM(kk) c enddo c enddo c endif c enddo return endif if(p.le.n) then kk=1 do i=1,nlc k=LC(i) Vap(i)=A(k,p) if(DABS(Vap(i)).GE.deps) then Vaps(kk+1)=i kk=kk+1 endif enddo Vaps(1)=kk-1 endif c write(6,*)'asssss',p,As(1,p) c write(6,*)'MBBB',MB(1,1),MB(2,1),MB(3,1),MBS(1,1) call CalculdCA(A,MB,MBS,Vap,Vapp,Vaps,SM,SMM,SMMM,Iposv, + p,pp,nl1,L1,L3,Lp1,Lp2,LL1,nll1,nl3,LC,nlc,m,m1,m2,m3,n) c write(6,*)'llllllllll1',nl1,nll1 call SIMPd2(A,As,MB,MBs,SMMM,Vapp,Nvar,qq,q,LL1,lp2,nll1,born,m1, + Iposv,L1,Lc,m2,nl1,m3,m,n) c write(14,*)'q,qq',q,qq c write(6,*)'llllllllll1',nl1,nll1 if (born.eq.1) then write(6,*) 'the linear problem is unbounded, add for example + the constraint xi<100000, i=1,nx ' write(14,*) 'the linear problem is unbounded, add for example + the constraint xi<100000, i=1,nx ' return endif c L2(q)=0 c do i=1,n c write(6,*)'vap(',i,')=',vap(i) c enddo kkk=1 if (qq.le.nll1) then ii=iposv(q)-n else ii=iposv(q) endif if(ii.GT.m1) then do i=1,nl1 k=L1(i) kk=iposv(k) Avec(i)=-A(ii,kk) if (DABS(Avec(i)).GE.deps) then Avecs(kkk+1)=i kkk=kkk+1 endif enddo Avecs(1)=kkk-1 kkk=nll1 else do i=1,nl1 k=L1(i) kk=iposv(k) c write(6,*)'kkkkk',kk,ii Avec(i)=A(ii,kk) if (DABS(Avec(i)).GE.deps) then Avecs(kkk+1)=i kkk=kkk+1 endif enddo Avecs(1)=kkk-1 kkk=nll1 endif call MJMBd(Avec,Avecs,Vap,Vapp,Vaps,q,qq,p,pp,MB,MBs, + MBsq,SMM,SMMM,L1,LL1,L3,Lp1,m,m1,m2,m3,n,nl1,nll1,nl3, + iposv,LC,Lp2,nlc) c if( (p.LE.n).and.(qq.LE.kkk)) then c Iposv(q)=p c NVAR(p)=1 c else kk=iposv(q) if (kk.LE.n) NVAR(kk)=0 Iposv(q)=p if (p.LE.n) NVAR(p)=1 L3(pp)=kk c endif goto 30 return end Subroutine Calculdcr(A,As,MB,Mbs,iposv,Cc,Cr,L1,L2, + L3,Lc,LL1,nll1,nlc,nl1,nl2,nl3,m1,m2,m3,m,n) Integer m1,m2,m3,m,n,nl2,iposv(m),L1(m), + L2(m),L3(n+m),nl3,nl1,Lams(m+1),As(m+1,n),Mbs(n+1,n), + k,kk,kkk,kkkk,kkkkk,LC(m),nlc,LL1(m),nll1 Real*8 A(m,n),Mb(n+1,n),Cr(n+m2),rrrs,Cc(n+m),Lam(m),deps, + Ccc(m) c probleme artificiel c L3 represente les variables hors base c nl3 nombre de variables hors base c L1 represente les variables artificielle de base c nl1 represente le nbe de variables artificielles courant do i=1,m lam(i)=0.d0 enddo do i=1,nl1 ii=iposv(l1(i)) Ccc(i)=Cc(ii) enddo deps=1.d-14 c Phase deux do ii=1,nll1 jj=iposv(LL1(ii))-n if(jj.Le.m1+m2) then if(jj.le.m1) then do iii=1,nl1 jjj=iposv(L1(iii)) c write(6,*)'lllllll',iii,jjj,jj Ccc(iii)=Ccc(iii)-Cc(jj+n)*A(jj,jjj) enddo else do iii=1,nl1 jjj=iposv(L1(iii)) Ccc(iii)=Ccc(iii)+Cc(jj+n)*A(jj,jjj) enddo endif endif enddo kkkk=1 do i=1,nlc rrrs=0.d0 k=LC(i) kk=MBs(1,i) do j=1,kk kkk=MBs(j+1,i) rrrs=rrrs+Ccc(kkk)*MB(kkk,i) enddo c write(6,*)'kkkk',k Lam(k)=rrrs if(DABS(Lam(k)).GE.deps) then lams(kkkk+1)=k kkkk=kkkk+1 endif enddo do i=1,nll1 k=LL1(i) kk=iposv(k) kkk=kk-n if (Cc(kk).GT.deps) then if(kkk.Le.m1) then Lam(kkk)=Cc(kk) else Lam(kkk)=-Cc(kk) endif lams(kkkk+1)=kkk kkkk=kkkk+1 endif enddo lams(1)=kkkk-1 do i=1,nl3 k=L3(i) rrrs=0.d0 if(k.GT.n) then if (k.Le.n+m1) then rrrs=-lam(k-n) else rrrs=lam(k-n) endif Cr(i)=rrrs else if(lams(1).LE.As(1,k)) then kk=Lams(1) do j=1,kk kkk=lams(j+1) rrrs=rrrs+lam(kkk)*A(kkk,k) enddo else kk=As(1,k) do j=1,kk kkk=As(j+1,k) rrrs=rrrs+lam(kkk)*A(kkk,k) enddo endif Cr(i)=Cc(k)-rrrs endif enddo return end Subroutine CalculdCA(A,MB,MBS,Vap,Vapp,Vaps,SM,SMM,SMMM,Iposv, + p,pp,nl1,L1,L3,Lp1,Lp2,LL1,nll1,nl3,LC,nlc,m,m1,m2,m3,n) Integer MBs(n+1,n),Vaps(n+1),L1(m),nl1,m,m1,m2,m3,n,p,pp,k,kk, + Iposv(m),L3(n+m),nl3,LP1(m),LL1(m),nll1,LC(m),nlc, + kkk,kkkk,kkkkk,Lp2(m) Real*8 MB(n+1,n),Vap(n+1),Vap1(n),deps,Vapp(m),SMM(n+1),Sm(m+1), + rrrs ,rrrs1,A(m,n),SMMM(m) do i=1,n vap1(i)=0.d0 enddo do i=1,m vapp(i)=0.d0 SMMM(i)=0.d0 enddo deps=1.d-14 if (p.GT.n) then if(p.LE.n+m1) then k=p-n kk=Lp2(k) kkk=MBs(1,kk) do i=1,kkk kkkk=MBs(i+1,kk) Vap(kkkk)=MB(kkkk,kk) kkkkk=L1(kkkk) Vapp(kkkkk)=Vap(kkkk) SMMM(kkkkk)=SMM(kkkk) Vaps(i+1)=kkkk enddo Vaps(1)=kkk do i=1,nll1 kkk=LL1(i) ii=iposv(kkk)-n rrrs=0.d0 do j=1,MBs(1,kk) jj=MBs(j+1,kk) jjj=L1(jj) jjjj=iposv(jjj) rrrs=rrrs+A(ii,jjjj)*MB(jj,kk) enddo if(iposv(kkk).GT.n+m1) then Vapp(kkk)=rrrs else Vapp(kkk)=-rrrs endif rrrs1=0.d0 do j=1,nl1 jj=L1(j) jjj=iposv(jj) rrrs1=rrrs1+A(ii,jjj)*SMM(j) enddo kkkk=iposv(kkk) kkkkk=kkkk-n if(kkkkk.GT.m1) then SMMM(kkk)=rrrs1-SM(kkkkk) else SMMM(kkk)=-rrrs1+SM(kkkkk) endif enddo else k=p-n kk=lp2(k) kkk=MBs(1,kk) do i=1,kkk kkkk=MBs(i+1,kk) Vap(kkkk)=-MB(kkkk,kk) kkkkk=L1(kkkk) Vapp(kkkkk)=Vap(kkkk) SMMM(kkkkk)=SMM(kkkk) Vaps(i+1)=kkkk enddo Vaps(1)=kkk do i=1,nll1 kkk=LL1(i) ii=iposv(kkk)-n rrrs=0.d0 do j=1,MBs(1,kk) jj=MBs(j+1,kk) jjj=L1(jj) jjjj=iposv(jjj) rrrs=rrrs+A(ii,jjjj)*MB(jj,kk) enddo if(iposv(kkk).GT.n+m1) then Vapp(kkk)=-rrrs else Vapp(kkk)=rrrs endif rrrs1=0.d0 do j=1,nl1 jj=L1(j) jjj=iposv(jj) rrrs1=rrrs1+A(ii,jjj)*SMM(j) enddo kkkk=iposv(kkk) kkkkk=kkkk-n if(kkkkk.GT.m1) then SMMM(kkk)=rrrs1-SM(kkkkk) else SMMM(kkk)=-rrrs1+SM(kkkkk) endif enddo endif else k=Vaps(1) do i=1,k kk=Vaps(i+1) kkk=MBs(1,kk) do j=1,kkk kkkk=MBs(j+1,kk) Vap1(kkkk)=Vap1(kkkk)+Vap(kk)*MB(kkkk,kk) enddo enddo kk=1 do i=1,nl1 k=L1(i) Vap(i)=Vap1(i) Vapp(k)=Vap(i) if(Dabs(Vap(i)).Ge.deps) then Vaps(kk+1)=i kk=kk+1 endif enddo Vaps(1)=kk-1 do i=1,nll1 kkk=LL1(i) ii=iposv(kkk)-n rrrs=0.d0 do j=1,Vaps(1) jj=Vaps(j+1) jjj=L1(jj) jjjj=iposv(jjj) c Write(6,*)'jjjjjj',jjj,jjjj,kkk,jj rrrs=rrrs+A(ii,jjjj)*Vap(jj) enddo rrrs1=0.d0 do j=1,nl1 jj=L1(j) SMMM(jj)=SMM(j) jjj=iposv(jj) rrrs1=rrrs1+A(ii,jjj)*SMM(j) enddo kkkk=iposv(kkk)-n if(kkkk.GT.m1) then Vapp(kkk)=rrrs-A(kkkk,p) SMMM(kkk)=rrrs1-SM(kkkk) else Vapp(kkk)=-rrrs+A(kkkk,p) SMMM(kkk)=-rrrs1+SM(kkkk) endif enddo endif c do i=1,m c write(14,*)'SMMM,Vapp',SMMM(i),Vapp(i) c enddo return end Subroutine MJMBd(Avec,Avecs,Vap,Vapp,Vaps,q,qq,p,pp,MB,MBs, + MBsq,SMM,SMMM,L1,LL1,L3,Lp1,m,m1,m2,m3,n,nl1,nll1,nl3, + iposv,LC,Lp2,nlc ) Integer L1(m),MBs(n+1,n),q,Vaps(n+1),m1,m2,m3,m,n,nl1, + MBsq(n+1),kk,Avecs(n+1),Lp1(m),qq,p,nnl1,nt,qp,LL1(m),nll1, + L3(n+m),nl3,Lp2(m),LC(m),nlc,pp,qqq,iposv(m),k,kkk,kkkk, + kkkkk,ittt Real*8 Vap(n+1),MB(n+1,n),SMMM(m),SMM(n+1),Avec(n),deps,dpiv, + rrrs,Vapp(m) deps=1.d-14 nt=0 ntt=0 if (p.LE.n) then if( qq.GT.nll1) then qp=qq-nll1 c write(14,*)'11111111111111111111111' goto 111 else c write(14,*)'22222222222222222222222' nnl1=nl1+1 qp=nnl1 do i=1,nl1 rrrs=0.d0 do j=1,Avecs(1) k=Avecs(j+1) c kk=L1(k) c kkk=Lp1(kk) rrrs=rrrs+Avec(k)*MB(k,i) enddo MB(nnl1,i)=-rrrs MB(i,nnl1)=0.d0 if(Dabs(rrrs).GE.deps) then MBs(1,i)=MBs(1,i)+1 kkkkk=MBs(1,i) MBs(kkkkk+1,i)=nnl1 endif enddo SMM(nnl1)=SMMM(q) if(iposv(q).GT.n+m1) then MB(nnl1,nnl1)=-1.d0 else MB(nnl1,nnl1)=1.d0 endif MBs(1,nnl1)=1 MBs(2,nnl1)=nnl1 Vap(nnl1)=Vapp(q) if(Dabs(Vap(nnl1)).GE.deps) then Vaps(1)=Vaps(1)+1 kkkkk=Vaps(1) Vaps(kkkkk+1)=nnl1 endif c do i=pp,nl3-1 c L3(i)=L3(i+1) c enddo do i=qq,nl1+nll1-1 LL1(i)=LL1(i+1) enddo nll1=nll1-1 nl1=nl1+1 nlc=nlc+1 L1(nl1)=q kk=iposv(q)-n LC(nl1)=kk LL1(nl1+nll1)=q Lp1(q)=nl1 Lp2(kk)=nl1 c nl3=nl3-1 endif else if (qq.GT.nll1) then nt=1 c write(14,*)'3333333333333333' qp=qq-nll1 goto 111 else c write(14,*)'4444444444444444444' qqq=Lp2(p-n) qp=nl1+1 ii=iposv(q)-n LC(qqq)=ii Lp2(ii)=qqq Lp2(p-n)=0 do i=1,nl1 rrrs=0.d0 do j=1,Avecs(1) k=Avecs(j+1) kk=L1(k) kkk=Lp1(kk) c write(14,*)'uuuuuuuuuuuu',k,kkk rrrs=rrrs+Avec(k)*MB(k,i) enddo c Write(6,*)'44444',n,nl1,qp MB(qp,i)=-rrrs MB(i,qqq)=0.d0 enddo SMM(qp)=SMMM(q) c Vap(qp)=MB(qp,qqq) Vap(qp)=Vapp(q) if(iposv(q).GT.n+m1) then MB(qp,qqq)=-1.d0 else MB(qp,qqq)=1.d0 endif MBs(1,qqq)=1 MBs(2,qqq)=qp Vaps(1)=Vaps(1)+1 kkkkk=Vaps(1) c write(*,*)'kkkkkkk',kkkkk Vaps(kkkkk+1)=qp endif endif 111 continue k=1 do i=1,nl1 if (DABS(MB(qp,i)).GE.deps) then MBsq(k+1)=i k=k+1 endif enddo MBsq(1)=k-1 DPiv=1.d0/Vap(qp) c write(14,*)'dddddddddddddpiv',dpiv k=MBsq(1) do i=1,Vaps(1) kk=Vaps(i+1) c k=MBsq(1) c write(14,*)'vap',vap(kk) if((kk.LT.qp).or.(kk.GT.qp)) then rrrs=-Vap(kk)*DPiv do j=1,k kkk=MBsq(j+1) MB(kk,kkk)=MB(kk,kkk)+rrrs*MB(qp,kkk) enddo Smm(kk)=Smm(kk)+rrrs*Smm(qp) else ntt=1 endif enddo if(ntt.eq.1) then rrrs=DPiv do j=1,k kkk=MBsq(j+1) MB(qp,kkk)=rrrs*MB(qp,kkk) enddo Smm(qp)=rrrs*Smm(qp) c write(14,*)'ooooooooooooooooooooui' endif do i=1,k kkk=1 kk=MBsq(i+1) do j=1,nl1 if(DABS(MB(j,kk)).GE.deps) then MBs(kkk+1,kk)=j kkk=kkk+1 endif enddo MBs(1,kk)=kkk-1 enddo if(nt.eq.1) then qqq=qq-nll1 qqqq=Lp2(p-n) Lp2(p-n)=0 do i=1,nl1 do j=qqqq,nl1-1 MB(i,j)=MB(i,j+1) MBs(i+1,j)=MBs(i+1,j+1) enddo enddo c Lp2(LC(qqqq))=0 do i=qqqq,nl1-1 MBs(1,i)=MBs(1,i+1) LC(i)=LC(i+1) Lp2(LC(i))=i enddo do i=qqq,nl1-1 SMM(i)=SMM(i+1) do j=1,nl1 MB(i,j)=MB(i+1,j) enddo L1(i)=L1(i+1) Lp1(L1(i))=i enddo do i=qq,nll1+2,-1 LL1(i)=LL1(i-1) enddo do i=1,nl1-1 ittt=0 do j=1,MBs(1,i) if (MBs(j+1,i).EQ.qqq) then ittt=1 goto 13 endif if (MBs(j+1,i).GT.qqq) then ittt=2 goto 13 endif enddo 13 continue if(ittt.eq.1) then do jj=j,MBs(1,i)-1 MBs(jj+1,i)=MBs(jj+2,i)-1 enddo MBs(1,i)=MBs(1,i)-1 endif if(ittt.eq.2) then do jj=j,MBs(1,i) MBs(jj+1,i)=MBs(jj+1,i)-1 enddo endif enddo nll1=nll1+1 LL1(nll1)=q nl1=nl1-1 nlc=nlc-1 c write(6,*)'nnnnnnlll1',nl1,nll1 endif c do i=1,nl1 c write(14,*)'SMMM',SMM(i) c enddo return end c ****************************************************** Subroutine SIMPd1(Cr,L3,nl3,Binf,pp,p,m1,m2,m3,m,n) Integer p,pp,m1,m2,m3,n,L3(n+m),nl3,k Real*8 Cr(n+m2),Binf,test Binf=Cr(1) p=L3(1) pp=1 do 121 k=2,nl3 TEST=Cr(k)-Binf if(test.LT.0.d0)then Binf=Cr(k) p=L3(k) pp=k endif 121 continue return end c **************************************************************** c Locate a pivot element taking degeneracy into account Subroutine SIMPd2(A,As,MB,MBs,SM,Vap,Nvar,qq,q,L1,lp2,nl1, + born,m1,Iposv,LL1,Lc,m2,nll1,m3,m,n) Integer L1(m),nl1,q,m1,m2,m3,n,born,MBs(n+1,n),As(m+1,n),qq, + Lp2(m),nll1,Nvar(n),Iposv(m),LL1(n),k,kk,kkk,Lc(n) real*8 SM(m),Vap(m),Q1,Q0,Q2,QP,MB(n+1,n),A(m,n),Va1(m),Va2(m), + rrrs,rrrs1,rrrs2 deps=1.d-14 itt=0 ittt=0 c write(*,*)'nl1',nll1,nl1 iii=nl1+nll1 do 11 i=1,iii k=L1(i) c write(6,*)'Vap',k,vap(k) c write(24,*)' A(',L2(i)+1,',',KP+1,')=',A(L2(i)+1,KP+1) if(Vap(k).GT.1.d-12) GOTO 2 11 continue born=1 write(6,*)'bornnn',born return 2 Q1=Sm(k)/Vap(k) q=k qq=i if(i+1.GT.iii)return Do 13 ii=i+1,iii k=L1(ii) if(Vap(k).GT.1.d-12) then Q2=Sm(k)/Vap(k) if (Q2.LT.Q1) then q=k qq=ii Q1=Q2 c We have a degenery elseif (Q2.EQ.Q1) then c write(14,*)'deeeeeeeeeeeeeeee' do 12 j=1,N rrrs1=0.d0 rrrs2=0.d0 If(Nvar(j).eq.1) then If(Iposv(k).eq.j) then rrrs1=1.d0 rrrs2=0.d0 elseif(Iposv(q).eq.j) then rrrs1=0.d0 rrrs2=1.d0 endif else if(qq.GT.nl1) then kk=qq-nl1 do jj=1,As(1,j) jjj=As(jj+1,j) jjjj=lp2(jjj) if (jjjj.GE.1) then rrrs1= rrrs1+MB(kk,jjjj)*A(jjj,j) endif enddo else kkk=L1(qq) kkk=iposv(kkk)-n if(itt.eq.0) then do kk=1,m Va1(kk)=0.d0 enddo do kk=1,nll1 rrrs=0.d0 kkkk=LC(kk) c write(6,*)'kkkkk',kk do jj=1,MBs(1,kk) jjj=MBs(jj+1,kk) jjjj=LL1(jjj) jjjjj=iposv(jjjj) rrrs=rrrs+A(kkk,jjjjj)*MB(jjj,kk) enddo Va1(kkkk)=rrrs enddo itt=1 endif do jj=1,AS(1,j) jjj=As(jj+1,j) rrrs1= rrrs1+Va1(jjj)*A(jjj,j) enddo rrrs1=rrrs1-A(kkk,j) if (kkk.Le.m1) rrrs1=-rrrs1 endif if(ii.GT.nl1) then kk=ii-nl1 do jj=1,As(1,j) jjj=As(jj+1,j) jjjj=lp2(jjj) if (jjjj.GE.1) then rrrs2=rrrs2+MB(kk,jjjj)*A(jjj,j) endif enddo else if(ittt.eq.0) then kkk=L1(ii) kkk=iposv(kkk)-n do kk=1,m Va2(kk)=0.d0 enddo do kk=1,nll1 rrrs=0.d0 kkkk=Lc(kk) do jj=1,MBs(1,kk) jjj=MBs(jj+1,kk) jjjj=LL1(jjj) jjjjj=iposv(jjjj) rrrs=rrrs+A(kkk,jjjjj)*MB(jjj,kk) enddo Va2(kkkk)=rrrs enddo ittt=1 endif do jj=1,AS(1,j) jjj=As(jj+1,j) rrrs2= rrrs2+Va2(jjj)*A(jjj,j) enddo rrrs2=rrrs2-A(kkk,j) if (kkk.Le.m1) rrrs2=-rrrs2 endif endif QP=rrrs1/Vap(q) Q0=rrrs2/Vap(k) if(Q0.NE.QP) GOTO 6 12 continue 6 if(Q0.GT.QP) then q=k qq=ii endif endif endif 13 continue Return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine cvd(Xp,Cp,A,Ad,Sm,Smd,Cc,Ccd,m1,m2,m3,m,md1,md2,md3,n) $debug Integer Xp(n),Cp(m),md1,md2,md3 Real*8 A(m,n),Ad(n,m),Cc(n+m),Ccd(m+n),Sm(m+1),Smd(n+1),rrrs md1=0 md2=0 md3=0 do i=1,n+m Ccd(i)=0.d0 enddo do i=1,n if(Xp(i).eq.0) then md3=md3+1 else if(Cc(i).Le.0.d0) then md2=md2+1 else md1=md1+1 endif endif enddo nd1=0 nd2=0 nd3=0 do i=1,n if(Xp(i).eq.0) then nd3=nd3+1 if (Cc(i).le.0.d0) then k=md1+md2+nd3 do j=1,m Ad(k,j)=-A(j,i) enddo Cp(k)=i Smd(k)=-Cc(i) else k=md1+md2+nd3 do j=1,m Ad(k,j)=A(j,i) enddo Cp(k)=i Smd(k)=Cc(i) endif else if(Cc(i).le.0.d0) then nd2=nd2+1 k=md1+nd2 do j=1,m Ad(k,j)=-A(j,i) enddo Cp(k)=i Smd(k)=-Cc(i) else nd1=nd1+1 k=nd1 do j=1,m Ad(k,j)=A(j,i) enddo Cp(k)=i Smd(k)=Cc(i) endif endif enddo do i=1,m Ccd(i)=-Sm(i) enddo do i=1,m1 do j=1,n Ad(j,i)=-Ad(j,i) enddo Ccd(i)=-Ccd(i) enddo do i=1,n rrrs=0.d0 do j=m1+m2+1,m rrrs=rrrs+Ad(i,j) enddo c write(6,*)'iiiii',i Smd(i)=Smd(i)+rrrs*1.d2 enddo c write(6,*)'md1,md2,md3',md1,md2,md3 Return end c **************************************************************** c ************************************************************* c fonction objective c ************************************************************ subroutine OBJE(nx,X,FO) integer nx real*8 X(nx),FO FO=0 do 25 i=1,nx-20 FO=FO+DEXP((X(i)*X(i)-4))*(X(i)-4) 25 continue do i=nx-20,nx-1 FO=FO+((X(i)*X(i)-1))*(X(i)-1) enddo fo=fo+(X(nx)+1)*(X(nx)+1) return end c c subroutine OBJE(nx,X,FO) c integer nx c real*8 X(nx),FO,sq c sq=1.D0 c do 47 i=1,nx c sq=sq*X(i) c 47 continue c FO=DEXP(sq) c return c end c ******************************************************** c fonction de contraintes c ******************************************************** c subroutine CONTR(mc,nx,X,FC) c integer mc,nx c real*8 X(nx),FC(mc),sq c sq=0.D0 c do 42 i=1,nx c sq=sq+(X(i)*X(i)) c 42 continue c FC(1)=sq-10.D0 c FC(2)=X(2)*X(3)-(5.D0*X(4)*X(5)) c FC(3)=(X(1)**3.D0)+(X(2)**3.D0)+1.D0 c return c end c ********************************************* subroutine CONTR(mc,nx,X,FC) integer mc,nx real*8 X(nx),FC(mc) do 11 i=1,mc-1 FC(i)= (X(i)*X(i))+(X(i+1)*X(i+1))-5 11 continue FC(mc)= (X(nx)*X(nx))+(X(nx-1)*X(nx-1))-5 return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccc Last Line