SUBROUTINE MINOP(N,X,F,G,STEP,ACC,MAXFUN,NITYP,INIT) implicit none integer:: n,maxfun,nityp,i,j,k,nfcall,ngcall,iupd,iptest real:: step,acc,epssq,qsq,dss,f,fa,c,d,ggdiag,hdiag real:: gggg,hggg,dg,dga,dggd,wbsq,b0,c1,c2,c3 real:: a,b,e,cc,ca,cb,theta,dsq,flag,sum,zwc,zgv,vgv,vsq real:: bb,ctemp,zv,phi,usq,gsq,uv,sigma,bound,wasq,zgz,vwc logical:: init integer,parameter:: mq=100 real:: X(*), G(*),XA(mq), GA(mq),H(mq,mq), GG(mq,mq) real:: V(mq),U(mq), WA(mq), WB(mq), WC(mq), WT(mq),WS(mq) real:: Y(mq), Z(mq), GZ(mq), GV(mq), PZ(mq) real:: sum_zwc !another div/0 preventer real:: accmbb !yet another div/0 prot real:: ggtemp ! " ! ****************************************************************** ! * INITIATLIZE THE FOLLOWING PARAMETERS * ! * NFCALL & NGCALL = NO. OF FUNCTION & GRADIENT CALLS RESP.. * ! ****************************************************************** if(n.lt.1) then print *,' minop n too small ',n return endif QSQ=1. !FOR COMPILER NFCALL = 1 NGCALL = 1 IUPD = 1 EPSSQ = ACC * ACC IPTEST = 1 DSS = STEP*STEP CALL CALCF (N,X,F) CALL CALCG (N,X,G) IF(.not. INIT) GO TO 8 !SAVE HESSIAN IF NOT FIRST RUN GSQ = 0 DO 3 I = 1,N 3 GSQ = GSQ + G(I) * G(I) ! ! ****************************************************************** ! * PREPARE FOR THE FIRST ITERATION BY INITIALIZING MATRICES H&GG * ! ****************************************************************** ! c=0 ; if(gsq.gt.0) C = -SQRT( DSS / GSQ ) GGDIAG = 0.01 * SQRT( GSQ) / STEP hdiag=0 ; if(ggdiag.gt.0) HDIAG = 1. / GGDIAG DO 6 I= 1,N DO 5 J = 1,N H(I,J) = 0.0D0 GG(I,J) = 0.0D0 5 CONTINUE H(I,I) = HDIAG GG(I,I) = GGDIAG 6 Z(I) = C * G(I) 8 GSQ = 0 DO 9 I = 1, N 9 GSQ = GSQ + G(I) * G(I) IF(GSQ.LT.EPSSQ) RETURN ! ****************************************************************** ! * TEST WHETHER MAXFUN CALLS OF CALCFG HAVE BEEN MADE * ! ****************************************************************** ! 18 IF(NFCALL.EQ.1) continue !CALL OOO(0,NFCALL,F) IF(NFCALL.GT.MAXFUN) RETURN !FUNC CALLS EXHAUSTED BOUND = SQRT (DSS) !WAS FOR PRINTOUT ! ****************************************************************** ! * CALCUALTE NEWTON POINT * ! ****************************************************************** 26 WASQ = 0.0D0 DO 28 I = 1,N WA(I) = 0.0D0 DO 27 J= 1,N 27 WA(I) = WA(I) - H(I,J) * G(J) 28 WASQ = WASQ + WA(I) * WA(I) IF ( WASQ .LE. DSS) GO TO 39 GGGG = 0 HGGG = 0.0D0 DO 30 I = 1, N WB(I) = 0.0D0 DO 29 J = 1,N 29 WB(I) = WB(I) + GG(I,J) * G(J) HGGG = HGGG - WA(I) * G(I) 30 GGGG = GGGG + WB(I) * G(I) D = HGGG * GGGG if(abs(d).lt.1e-20) then d=gsq*gsq*1e10 !div 0 protect else D = GSQ * GSQ / D endif D = 0.2 + 0.8 * D if(abs(d).gt.1e6) goto 32 !protect overflow IF ( D * D * WASQ .GT. DSS ) GO TO 32 D = SQRT ( DSS / WASQ ) DO 31 I = 1,N 31 WA(I) = D * WA(I) GO TO 41 32 CONTINUE do i=1,20 if(abs(dss).lt.1e-12) dss=dss*3. if(abs(gggg).lt.1e-14) gggg=gggg*3. if(abs(gggg).gt.1e-12) gggg=gggg/3. enddo ! print *,'gggg dss qsq ',gggg,dss,qsq IF ((GGGG *(GGGG * DSS)).GT. GSQ ** 3 ) GO TO 34 ! ggtemp=gggg/gsq ! ggtemp=ggtemp*gggg/gsq ! ggtemp=ggtemp*dss/gsq ! IF (ggtemp .GT. 1.0) GO TO 34 ! Yet another patch to avoid overflow- sbp 11/08/2018 ! ! ****************************************************************** ! * THE CAUCHY POINT IS OUTSIDE CIRCLE, SO PICK THE POINT ON THE * ! * BOUNDARY * ! ****************************************************************** ! if((gsq*1.0e-25).lt.dss) then C = -SQRT(DSS)/sqrt(GSQ) ! else !avoid float underflow ! c=-1.0e-15 ! endif DO 33 I = 1,N 33 WA(I) = C * G(I) GO TO 41 ! ! ****************************************************************** ! * THE CAUCHY POINT IS INSIDE CIRCLE, USE DOG LEG STEP. * ! ****************************************************************** 34 DO 35 I = 1,N 35 WA(I) = D * WA(I) CA = 0.0 CB = 0 C = -GSQ / GGGG ! ****************************************************************** ! * SET THE STEEPEST DESCENT CORRECTION IN WB AND THE DIFFERENCE * ! * BETWEEN WA AND WB IN WC. * ! ****************************************************************** DO 36 I = 1, N WB(I) = C * G(I) WC(I)= WA(I) - WB(I) CA = CA + WB(I) * WC(I) 36 CB = CB + WC(I) * WC(I) ! ! * FIND CORRECTION VECTOR AT THE INTERSECTION OF WA-WB AND CIRCLE * C = DSS -C * C * GSQ if(( CA * CA + C * CB).gt.0) then !test sqrt neg val THETA = SIGN (C / (ABS(CA) + SQRT( CA * CA + C * CB)), CA) else goto 39 endif ! * TEST WHETHER TO USE THE GENERALIZED NEWTON STEP * IF (THETA-1.) 37, 39, 39 37 DO 38 I = 1, N 38 WA(I) = WB(I) + THETA * WC(I) GO TO 41 39 CONTINUE ! 41 DO 42 I = 1,N 42 XA(I) = X(I) + WA(I) NFCALL = NFCALL + 1 CALL CALCF(N,XA,FA) 43 DSQ = 0.0D0 DO 44 I = 1,N 44 DSQ = DSQ + WA(I) * WA(I) ! ****************************************************************** ! * TEST IF FUNCTION VALUE IS DECREASED * ! ****************************************************************** IF(MOD(NFCALL,NITYP).EQ.0) CALL OOO(0,NFCALL,FA) IF (FA .LT. F) GO TO 46 DSS = 0.25 * DSQ GO TO 18 !SBP MOD ! ! * CALCULATE SOME NUMBERS FOR REVISING STEP BOUND * 46 NGCALL = NGCALL + 1 CALL CALCG(N,XA,GA) DG = 0.0 DGA = 0.0 DGGD = 0.0 WBSQ = 0.0 DO 47 I = 1, N WC(I) = 0.0D0 DO 47 J = 1,N 47 WC(I) = WC(I) + GG(I,J) * WA(J) ! ! * SET (Y-GS) IN WB * DO 48 I = 1,N U(I) = WC(I) WB(I) = GA(I) -G(I) - WC(I) DG = DG + G(I) * WA(I) DGA = DGA + GA(I) * WA(I) DGGD = DGGD + WC(I) * WA(I) 48 WBSQ = WBSQ + WB(I) *WB(I) ! ****************************************************************** ! * TEST WHETHER TO DECREASE THE STEP-BOUND * ! ****************************************************************** IF (FA-F-0.1 * DG-0.05 * DGGD) 51, 51, 50 50 DSS = 0.25 * DSQ IF(DSS.LT.1E-13) RETURN !orig -15 GO TO 54 !NOT TOO SMALL ! ****************************************************************** ! * TEST WHETHER TO INCREASE THE STEP-BOUND * ! ****************************************************************** 51 DSS = DSQ IF (WBSQ - 0.25 * GSQ) 53, 53, 52 52 IF (DG -DGA -DGA) 54, 53, 53 53 DSS = 4.0 * DSQ ! * SET THE DIFFERENCE BETWEN GRADIENTS IN WC * 54 DO 55 I = 1,N 55 WC(I) = GA(I) - G(I) F = FA B0 = 0.0 DO 56 I = 1,N X(I) = XA(I) G(I) = GA(I) 56 B0 = B0 + WC(I) * WA(I) ! ****************************************************************** ! * TEST WHETHER WE NEED TO UPDATE H AND G * ! ****************************************************************** IF (B0 .GE. 1.0D-30) GO TO 58 ! WRITE(6,57) NFCALL 57 FORMAT('FUNCTION DECREASED ON',I4,'TH ITERATION, NO UPDATE IS MADE & ,SINCE Y*Z IS NEG. OR TOO SMALL') GO TO 8 58 IF (IUPD .EQ. 1) GO TO 63 IF (FLAG .EQ. 0.0) GO TO 60 DO 59 I = 1, N 59 Z(I) = WA(I) GO TO 63 60 C2=0.0 C3 = 0.0 SUM = 0.0 DO 61 I = 1, N C2 = C2 + WS(I) * WA(I) C3 = C3 + WT(I) * WA(I) 61 SUM = SUM +Y(I) * WA(I) DO I = 1, N sum_zwc=1. !more div/0 defense if(zwc.ne.0) sum_zwc=sum/zwc Z(I) = C1 *(C2 * Z(I) + C3 * V(I)) - sum_zwc * Z(I) z(i)=max(-1e15,min(1e15,z(i))) ENDDO 63 IUPD = IUPD+1 DO 64 I = 1, N GV(I) = U(I) - WC(I) V(I) = WA(I) GZ(I) = 0.0 DO 64 J = 1, N GZ(I) = GZ(I) + GG(I,J) * Z(J) 64 V(I) =V(I) - H(I,J) * WC(J) ZGV = 0.0 VGV = 0.0 ZWC = 0.0 ZGZ = 0.0 VWC = 0.0 DO 65 I = 1, N ZWC = ZWC + Z(I) * WC(I) ZGZ = ZGZ + Z(I) * GZ(I) VWC = VWC + V(I) * WC(I) ZGV = ZGV + Z(I) * GV(I) VGV = VGV + V(I) * GV(I) 65 CONTINUE zgz=min(1e10,max(-1e10,zgz)) vgv=min(1e10,max(-1e10,vgv)) ctemp=(zgz*vgv)-(zgv*zgv) if(ctemp.ne.0) then c1 = 1./ctemp else c1=.001 endif C2 = ZWC * VGV - VWC * ZGV C3 = -ZWC * ZGV + VWC * ZGZ c2=min(1e10,max(-1e10,c2)) c3=min(1e10,max(-1e10,c3)) VSQ =0.0 ZV =0.0 A = 0.0 E = 0.0 DO I = 1,N Y(I) = C1 * (C2 * GZ(I) + C3 * GV(I)) WS(I) = VGV * GZ(I) - ZGV * GV(I) WT(I) = -ZGV * GZ(I) + ZGZ * GV(I) ZV = ZV + Z(I) * V(I) VSQ = VSQ + V(I) * V(I) A = A + WS(I) * WA(I) E = E + WT(I) * WA(I) enddo B = 0.0 DO 67 I = 1,N PZ(I) = C1 * ( A * Z(I) + E * V(I)) 67 B = B +PZ(I) * Y(I) IF (B .GE. 1.0D-30) GO TO 70 B = B0 ! WRITE (6,68) 68 FORMAT(15X,'IMAGE PRODUCT NEG.OR TOO SMALL,ORIGINAL STEP USED') DO 69 I = 1, N PZ(I) = WA(I) 69 GA(I) = WC(I) GO TO 72 70 DO 71 I = 1, N 71 GA(I) = Y(I) 72 CC = 0.0 A = 0.0 DO 74 I = 1, N WB(I) = 0.0 GZ(I) = 0.0 DO 73 J = 1, N GZ(I) = GZ(I) + GG(I,J) *PZ(J) 73 WB(I) = WB(I) + H(I,J) * GA(J) A = A + GA(I) * WB(I) 74 CC = CC + PZ(I) * GZ(I) ! ****************************************************************** ! * OPTIMAL CONDITIONING * ! ****************************************************************** IF ( B .LE. 2.0D0 * A * CC / (A + CC)) GO TO 75 phi=1.0 if((b-a).ne.0) PHI = B / (B - A) GO TO 76 75 accmbb=a*cc-b*b !div/0 protect if(abs(accmbb).lt.1e-25) then if(accmbb.ge.0) then accmbb=1e-25 else accmbb=-1e-25 endif endif PHI = B * (CC - B)/accmbb 76 USQ = 0.0 UV = 0.0 ! ****************************************************************** ! * TEST WHETHER Z & V ARE LINEARLY INDEPENDENT * ! ****************************************************************** DO 77 I = 1, N if(vsq.ne.0) U(I) = Z(I) - ( ZV / VSQ ) * V(I) if(vsq.eq.0) U(I) = Z(I) USQ = USQ + U(I) * U(I) 77 UV = UV + U(I) * V(I) FLAG = 0.0 uv=min(1e10,max(-1e10,uv)) vsq=min(1e10,max(-1e10,vsq)) usq=min(1e10,max(-1e10,usq)) IF (1.0D6 * UV * UV .GE. VSQ * USQ ) FLAG =1.0 ! ****************************************************************** ! * IF FLAG = 1 THEN Z & V ARE LINEARLY DEPENDENT * ! ****************************************************************** bb=b*b if(abs(bb).lt.1e-20) bb=1e-20 !div by 0 protect if(abs(a).lt.1e-20) a=1e-20 !div by 0 protect SIGMA = 1.+ PHI * (A * CC/bb - 1.) IF ( abs(SIGMA) .LT. 1.0D-25) then WRITE (6,78),sigma 78 FORMAT(20X,' G IS ALMOST SINGULAR',e12.4) endif do I = 1, N if((a.ne.0).and.(b.ne.0)) then WA(I) = (1./B) * PZ(I) - (1./A) * WB(I) WC(I) = (1./B) * GZ(I) -CC/bb * GA(I) else wa(i)=0.0 ; wc(i)=0. endif U(I) = GA(I)- GZ(I) enddo ! ****************************************************************** ! * START TO UPDATE MATRICES H & GG * ! ****************************************************************** do I = 1, N do J = 1, N GG(I,J) = GG(I,J) + (1./B) * (U(I) * GA(J) + GA(I) * U(J)) & -(( B- CC ) / bb ) * GA(I) * GA(J) & -(PHI* A / SIGMA) * WC(I) * WC(J) H(I,J) = H(I,J) - (1./A) * WB(I) * WB(J) + (1./B) * PZ(I) * PZ(J) & + PHI * A * WA(I) * WA(J) enddo enddo GO TO 8 END !