OPTIMIZATION PART
2008/01/16 | J. SCHOENMAEKERS | NEW
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(optgra), | intent(inout) | :: | me | |||
| real(kind=wp), | intent(inout) | :: | Varacc |
ITERATION SCALED DISTANCE ACCUMULATED |
||
| integer(kind=ip), | intent(out) | :: | Finish |
0=LIMIT 1=OPTIM |
||
| real(kind=wp) | :: | Desnor | ||||
| logical, | intent(out) | :: | error |
if there was a fatal error |
subroutine ogopti(me,Varacc,Finish,Desnor,error) !! OPTIMIZATION PART !! !! 2008/01/16 | J. SCHOENMAEKERS | NEW class(optgra),intent(inout) :: me real(wp),intent(inout) :: Varacc !! ITERATION SCALED DISTANCE ACCUMULATED integer(ip),intent(out) :: Finish !! 0=LIMIT 1=OPTIM logical,intent(out) :: error !! if there was a fatal error integer(ip) :: staflg , faccnt , numcor integer(ip) :: con , var , cos , act , ind , len , inc integer(ip) :: nnn , typ , des , prv , met real(wp) :: Desnor , foldis , cosimp , cornor , quacor , refdis real(wp) :: co0 , co1 , co2 , nor real(wp) :: cosco2 , cosco1 real(wp) :: maxdis , norprv real(wp) :: val , max , det , ccc , dis real(wp) :: fac , del , exc , eps , imp real(wp) :: bet , tht character(len=str_len) :: str character(len=name_len) :: nam real(wp) , dimension(:) , allocatable :: cosact real(wp) , dimension(:) , allocatable :: varvec real(wp) , dimension(:) , allocatable :: varwrk real(wp) , dimension(:) , allocatable :: corvec real(wp) , dimension(:) , allocatable :: desder real(wp) , dimension(:) , allocatable :: desprv real(wp) , dimension(:) , allocatable :: varprv real(wp) , dimension(:) , allocatable :: convec real(wp) , dimension(:) , allocatable :: conqua integer(ip) , dimension(:) , allocatable :: concor error = .false. ! initialize: - JW : these could also be in the class. ! that would save the allocate/reallocate ! step every time this routine is called. allocate (cosact(me%Numcon)) allocate (varvec(me%Numvar)) allocate (varwrk(me%Numvar)) allocate (corvec(me%Numcon)) allocate (desder(me%Numvar)) allocate (desprv(me%Numvar)) allocate (varprv(me%Numvar)) allocate (convec(me%Numcon+1)) allocate (conqua(me%Numcon+1)) allocate (concor(me%Numcon+1)) ! ====================================================================== ! OPTIMIZATION PART ! ---------------------------------------------------------------------- cos = me%Numcon + 1 des = me%Numcon + 2 prv = me%Numcon + 3 faccnt = 0 me%Conred(1:cos,:) = me%Conder(1:cos,:) me%Conred(des,:) = me%Vardes desder = me%Vardes me%Conred(prv,:) = me%Funvar me%Conder(prv,:) = me%Funvar call me%ogwrit(3,"") call me%ogwrit(3,"OPTIMIZATION PART") ! ---------------------------------------------------------------------- ! WRITE (STR,'("NUMACT = ",I4)') NUMACT ! CALL me%ogwrit (3,STR) numcor = me%Numact concor = me%Actcon ! ACTCON = list of indices of active constraints me%Numact = 0 ! number of active constraints me%Conact(1:des) = 0 main: do ! DO COR = 1, NUMCOR ! CON = CONCOR(COR) ! NAM = CONSTR(CON) ! LEN = CONLEN(CON) ! WRITE (STR,'("ACT = ",I4,5X,1X,A)') CON, NAM(1:LEN) ! CALL me%ogwrit (3,STR) ! CALL me%OGINCL (CON) ! end do ! ====================================================================== ! VECTOR OF STEEPEST ASCENT ! ---------------------------------------------------------------------- call me%ogwrit(3,"") call me%ogwrit(3,"VECTOR OF STEEPEST ASCENT") call me%ogwrit(3,"") call me%ogwrit(3,"REMOVE AND INCLUDE CONSTRAINTS:") call me%ogwrit(3,"") call me%ogwrit(3," REM INC CONSTRAINT") ! ---------------------------------------------------------------------- ! REMOVE PASSIVE INEQUALITY CONSTRAINTS ! ---------------------------------------------------------------------- do act = me%Numact , 1 , -1 con = me%Actcon(act) ! CON=index of active constraint in CONVAL if ( me%Conval(con)<=1.0_wp ) cycle nam = me%Constr(con) len = me%Conlen(con) write (str,'(I4,5X,1X,A)') con , nam(1:len) call me%ogwrit(2,str) call me%ogexcl(act,error) ! Exclude constraint ACT from active set if (error) return me%Conact(con) = -1 ! Mark CON as blocked from inclusion in active set end do ! ---------------------------------------------------------------------- ! INCLUDE VIOLATED INEQUALITY CONSTRAINTS AND SELECT PASSIVE ONES ! SELECT PASSIVE INEQUALITY CONSTRAINTS ! ---------------------------------------------------------------------- do con = 1 , me%Numcon if ( me%Contyp(con)==-2 ) cycle if ( me%Conact(con)>0 ) then ! CON is active, CONACT(CON) is the position of CON in active set elseif ( me%Contyp(con)==0 ) then ! CON in inactive me%Conact(con) = 0 elseif ( me%Conval(con)<-1.0_wp ) then me%Conact(con) = 0 elseif ( me%Conval(con)<=+1.0_wp ) then me%Conact(con) = 0 else me%Conact(con) = -1 end if end do ! ====================================================================== nnn = 1 inner: do ! Active set refinement loop nnn = nnn + 1 ! JW : some kind of max iter here? should this be an input? if ( nnn>999 ) then ! saveguard active set refinement loop to continue forever Finish = 0 write (str,*) "NNN=" , nnn call me%ogwrit(2,str) exit main end if ! ====================================================================== ! DERIVATIVES OF MERIT W.R.T. ACTIVE CONSTRAINTS ! ---------------------------------------------------------------------- cosact = me%ogrigt(-me%Conred(cos,1:me%Numact)) ! COSACT=Lagrange multipliers??? Desnor = sqrt(sum(me%Conred(cos,me%Numact+1:me%Numvar)**2)) ! norm of cost function gradient ! ---------------------------------------------------------------------- ! CONSTRAINT REMOVAL ! ---------------------------------------------------------------------- ind = 0 exc = -1.0e-12_wp max = exc do act = 1 , me%Numact con = me%Actcon(act) ! ACT=index in active set, CON=index of constraint if ( me%Contyp(con)==0 ) cycle val = cosact(act) fac = dot_product(me%Conred(con,1:me%Numvar),me%Conred(con,1:me%Numvar)) fac = sqrt(fac) val = val*fac if ( val>=exc ) cycle if ( val>max ) cycle max = val ind = act end do ! ---------------------------------------------------------------------- if ( ind/=0 ) then con = me%Actcon(ind) nam = me%Constr(con) len = me%Conlen(con) write (str,'(I4,5X,3(1X,D10.3),1X,A)') con , Desnor , max , me%Varmax , nam(1:len) call me%ogwrit(3,str) call me%ogexcl(ind,error) ! Exclude constraint IND from active set if (error) return cycle end if ! ---------------------------------------------------------------------- ! CONSTRAINT INCLUSION ! ---------------------------------------------------------------------- if ( Desnor/=0.0_wp ) then ! ---------------------------------------------------------------------- inc = 0 eps = 1.0e-03_wp max = -1.0e10_wp ! JW : what is going on here ? max = 0.0_wp ! ---------------------------------------------------------------------- do con = 1 , me%Numcon if ( me%Contyp(con)==-2 ) cycle if ( me%Conact(con)/=0 ) cycle del = dot_product(me%Conred(con,me%Numact+1:me%Numvar),& me%Conred(cos,me%Numact+1:me%Numvar))/Desnor val = abs(del)*me%Varmax if ( val<eps ) cycle fac = dot_product(me%Conred(con,1:me%Numvar),me%Conred(con,1:me%Numvar)) fac = sqrt(fac) del = del/fac if ( del<0.0_wp .and. del<max ) then max = del inc = con end if if ( me%Contyp(con)/=0 ) cycle del = -del if ( del<0.0_wp .and. del<max ) then max = del inc = con end if end do ! ---------------------------------------------------------------------- if ( inc/=0 ) then con = inc nam = me%Constr(con) len = me%Conlen(con) write (str,'(5X,I4,3(1X,D10.3),1X,A)') con , Desnor , max , me%Varmax , nam(1:len) call me%ogwrit(3,str) call me%ogincl(inc) cycle end if end if ! ---------------------------------------------------------------------- ! MATCHED INEQUALITY CONSTRAINTS + STEEPEST ASCENT VECTOR NORM ! ---------------------------------------------------------------------- call me%ogwrit(3,"") call me%ogwrit(3,"STATUS OF MATCHED INEQUALITY CONSTRAINTS:") call me%ogwrit(3,"") call me%ogwrit(3," ACT PAS CONSTRAINT") do con = 1 , me%Numcon if ( me%Contyp(con)==-2 ) cycle nam = me%Constr(con) len = me%Conlen(con) if ( me%Contyp(con)==0 ) then elseif ( me%Conact(con)>0 ) then write (str,'(I4,5X,1X,A)') con , nam(1:len) call me%ogwrit(3,str) elseif ( me%Conact(con)==0 ) then write (str,'(5X,I4,1X,A)') con , nam(1:len) call me%ogwrit(3,str) end if end do call me%ogwrit(3,"") write (str,'("STEEPEST ASCENT NORM: ",D13.6)') Desnor call me%ogwrit(3,str) write (str,'("MAXIMUM DISTANCE....: ",D13.6)') me%Varmax call me%ogwrit(3,str) ! ====================================================================== Finish = 0 ! ====================================================================== ! IF CONVERGENCE ! ---------------------------------------------------------------------- cosimp = Desnor*me%Varmax if ( abs(cosimp)<=1.0_wp ) then foldis = 0.0_wp Finish = 1 write (str,'("FINAL...............:",1X,D13.6,'//'11X,1(1X,D10.3),1X,D16.9)') & foldis , cosimp , me%Conval(cos) + cosimp call me%ogwrit(2,str) exit main end if ! ====================================================================== ! IF CONSTRAINT IS HIT ! ---------------------------------------------------------------------- do var = 1 , me%Numvar val = me%Conder(cos,var) do act = 1 , me%Numact con = me%Actcon(act) val = val + me%Conder(con,var)*cosact(act) end do me%Vardes(var) = val end do ! ---------------------------------------------------------------------- ind = 0 dis = 1.0e10_wp do con = 1 , me%Numcon if ( me%Contyp(con)==-2 ) cycle if ( me%Conact(con)/=-1 ) cycle val = dot_product(me%Conred(con,me%Numact+1:me%Numvar),me%Conred(cos,me%Numact+1:me%Numvar)) if ( val==0.0_wp ) cycle val = -me%Conval(con)/val*Desnor if ( val<=0.0_wp ) cycle if ( val>=dis ) cycle dis = val ind = con end do ! ---------------------------------------------------------------------- if ( ind/=0 ) then val = sqrt(sum((me%Varval-me%Varref+me%Vardes*dis/Desnor)**2)) if ( val>me%Varmax ) ind = 0 end if ! ---------------------------------------------------------------------- if ( ind/=0 ) then if ( me%Confix(ind)<=0 ) then if ( val>me%Varmax*1.0e-1_wp ) ind = 0 end if end if ! ---------------------------------------------------------------------- if ( ind/=0 ) then imp = dis*Desnor con = ind tht = 1.0_wp bet = 0.0_wp nam = me%Constr(con) len = me%Conlen(con) write (str,'( "CONSTRAINT REACHED..:",'//'1X,D13.6,11X,1(1X,D10.3),1X,D16.9,22X,1X,I4,1X,A)') & dis , imp , me%Conval(cos) + cosimp , con , nam(1:len) call me%ogwrit(2,str) Varacc = Varacc + dis me%Varval = me%Varval + dis*me%Vardes/Desnor do con = 1 , me%Numcon + 1 val = dot_product(me%Conred(con,me%Numact+1:me%Numvar),me%Conred(cos,me%Numact+1:me%Numvar)) me%Conval(con) = me%Conval(con) + val*dis/Desnor end do cycle main end if exit inner end do inner inner2: do ! ====================================================================== cosact = me%ogrigt(-me%Conred(cos,1:me%Numact)) do var = 1 , me%Numvar val = me%Conder(cos,var) do act = 1 ,me%Numact con = me%Actcon(act) val = val + me%Conder(con,var)*cosact(act) end do desprv(var) = val end do Desnor = sqrt(sum(desprv**2)) write (str,'("DESNOR=",D13.6)') Desnor ! CALL me%ogwrit (2,STR) ! ---------------------------------------------------------------------- cosact = me%ogrigt(-me%Conred(prv,1:me%Numact)) do var = 1 , me%Numvar val = me%Conder(prv,var) do act = 1 , me%Numact con = me%Actcon(act) val = val + me%Conder(con,var)*cosact(act) end do varprv(var) = val end do norprv = sqrt(sum(varprv**2)) write (str,'("NORPRV=",D13.6)') norprv ! CALL me%ogwrit (2,STR) ! ---------------------------------------------------------------------- cosact = me%ogrigt(-me%Conred(des,1:me%Numact)) do var = 1 , me%Numvar val = desder(var) do act = 1 , me%Numact con = me%Actcon(act) val = val + me%Conder(con,var)*cosact(act) end do me%Vardir(var) = val end do nor = sqrt(sum(me%Vardir**2)) write (str,'("NOR=",D13.6)') nor ! CALL me%ogwrit (2,STR) met = me%Optmet tht = 1.0_wp bet = 0.0_wp select case (met) case ( 0 ) ! STEEPEST DESCENT METHOD case ( 1 ) ! MODIFIED SPECTRAL CONJUGATE GRADIENT METHOD varvec = desprv - varprv if ( norprv**2>1.0e-12_wp ) then tht = -dot_product(me%Vardir,varvec)/norprv**2 bet = Desnor**2/norprv**2 end if case ( 2 ) ! SPECTRAL CONJUGATE GRADIENT METHOD varvec = desprv - varprv varwrk = me%Varref - me%Vargrd val = dot_product(varwrk,varvec) fac = dot_product(me%Vardir,varvec) if ( abs(val)>1.0e-12_wp .and. abs(fac)>1.0e-12_wp ) then tht = -dot_product(varwrk,varwrk)/val varwrk = -varvec*tht - varwrk bet = dot_product(varwrk,desprv)/fac end if case ( 3 ) ! CONJUGATE GRADIENT METHOD if ( norprv/=0.0_wp ) then tht = 1.0_wp bet = Desnor**2/norprv**2 end if end select ! WRITE (STR,'("THT=",D13.6)') THT ! CALL me%ogwrit (3,STR) ! WRITE (STR,'("BET=",D13.6)') BET ! CALL me%ogwrit (3,STR) ! ---------------------------------------------------------------------- eps = 1.0e-03_wp ! WRITE (STR,*) "THT/BET=",THT,BET ! CALL me%ogwrit (2,STR) me%Vardes = tht*desprv + bet*me%Vardir Desnor = sqrt(sum(me%Vardes**2)) nor = Desnor do con = 1 , me%Numcon if ( me%Contyp(con)==-2 ) cycle if ( me%Conact(con)/=0 ) cycle del = dot_product(me%Conder(con,1:me%Numvar),me%Vardes(1:me%Numvar))/nor val = abs(del)*me%Varmax if ( val<eps ) cycle fac = dot_product(me%Conder(con,1:me%Numvar),me%Conder(con,1:me%Numvar)) del = del/sqrt(fac) nam = me%Constr(con) len = me%Conlen(con) typ = me%Contyp(con) if ( del<0.0_wp ) then ! CALL OGINCL (CON) act = me%Conact(con) ! WRITE (STR,'(5X,2I4,3(1X,D10.3),1X,A)') & ! CON,ACT,+NOR,VAL,DEL,NAM(1:LEN) ! CALL me%ogwrit (2,STR) if ( act/=0 ) cycle inner2 bet = 0.0_wp tht = 1.0_wp end if if ( me%Contyp(con)/=0 ) cycle del = -del if ( del<0.0_wp ) then ! CALL OGINCL (CON) act = me%Conact(con) ! WRITE (STR,'(5X,2I4,3(1X,D10.3),1X,A)') & ! CON,ACT,-NOR,VAL,DEL,NAM(1:LEN) ! CALL me%ogwrit (2,STR) if ( act/=0 ) cycle inner2 bet = 0.0_wp tht = 1.0_wp del = dot_product(me%Conder(con,1:me%Numvar),me%Conder(cos,1:me%Numvar)) val = abs(del)*me%Varmax/Desnor ! WRITE (STR,'(5X,2I4,3(1X,D10.3),1X,A)') & ! CON,TYP,-DESNOR,VAL,DEL,NAM(1:LEN) ! CALL me%ogwrit (2,STR) end if end do cosco1 = dot_product(me%Vardes,me%Conder(cos,1:me%Numvar))/Desnor if ( cosco1<0.0_wp ) then write (str,*) "COSCO1=" , cosco1 call me%ogwrit(2,str) bet = 0.0_wp tht = 1.0_wp end if me%Vardes = tht*desprv + bet*me%Vardir Desnor = sqrt(sum(me%Vardes**2)) exit inner2 end do inner2 inner3: do ! WRITE (STR,*) "THT/BET=",THT,BET ! CALL me%ogwrit (2,STR) ! ====================================================================== ! SECOND ORDER DERIVATIVES BY NUMERICAL DIFFERENCING ! ---------------------------------------------------------------------- call me%ogwrit(3,"") write (str,'("SECOND ORDER CORRECTION")') call me%ogwrit(3,str) me%Vardes = tht*desprv + bet*me%Vardir Desnor = sqrt(sum(me%Vardes**2)) ! ---------------------------------------------------------------------- ! SAFEGUARD: ZERO STEEPEST-ASCENT NORM => TREAT AS CONVERGENCE ! ---------------------------------------------------------------------- ! NOTE: otherwise, DESNOR used in denominator later on may cause ! NaN issues ! see: https://github.com/esa/pyoptgra/pull/18/files IF (Desnor == 0.0_wp) THEN write (me%Loglup,'("OGOPTI WARN: Zero second order correction")') write (me%Loglup,'(" Treating as convergence")') foldis = 0.0_wp quacor = 0.0_wp cosimp = 0.0_wp Finish = 1 exit main END IF ! ---------------------------------------------------------------------- ! MAXIMUM TRAVEL DISTANCE ! ---------------------------------------------------------------------- dis = me%Varmax varvec = me%Varval - me%Varref co0 = dot_product(varvec,varvec) - dis**2 if ( co0>=0.0_wp ) then dis = 0.0_wp else co1 = dot_product(me%Vardes,varvec) co2 = Desnor**2 det = co1**2 - co0*co2 dis = (sqrt(det)-co1)/co2 end if dis = dis*Desnor maxdis = dis ! ====================================================================== ! COMPUTE SECOND ORDER EFFECTS ! ---------------------------------------------------------------------- nnn = 0 if ( me%Senopt>=+1 ) then do con = 1 , me%Numcon if ( me%Contyp(con)==-2 ) cycle act = me%Conact(con) ind = me%Senact(con) if ( act==-1 .and. ind==-1 ) cycle if ( act==0 .and. ind==0 ) cycle if ( act>0 .and. ind>0 ) cycle nam = me%Constr(con) len = me%Conlen(con) write (str,'(I4,1X,I4,1X,I4,1X,A)') con , act , ind , nam(1:len) call me%ogwrit(2,str) nnn = 1 end do end if if ( me%Senopt<=0 .or. nnn==1 ) then fac = me%Varstp/Desnor varvec = me%Varref + me%Vardes*me%Varstp/Desnor call me%ogeval(varvec,convec,0,me%Conder) conqua = matmul(me%Conder(1:cos,1:me%Numvar),me%Vardes(1:me%Numvar)) conqua = 2.0_wp*(convec-me%Conref-conqua*fac)/fac**2 end if if ( me%Senopt==-1 ) then me%Senqua = conqua elseif ( me%Senopt>=+1 .and. nnn==0 ) then conqua = me%Senqua end if ! ====================================================================== ! COMPUTE CORRECTION VECTOR ! ---------------------------------------------------------------------- do act = 1 , me%Numact con = me%Actcon(act) corvec(act) = conqua(con) end do corvec = me%ogleft(corvec) ! ---------------------------------------------------------------------- cornor = sqrt(sum(corvec(1:me%Numact)**2))*0.5_wp/Desnor/Desnor call me%ogwrit(3,"") write (str,'("STEEPEST ASCENT NORM: ",D13.6)') Desnor call me%ogwrit(3,str) write (str,'("ACCUMULATED DISTANCE: ",D13.6)') Varacc call me%ogwrit(3,str) ! ====================================================================== ! GET EXTREMUM ! ---------------------------------------------------------------------- cosco1 = dot_product(me%Vardes,me%Conder(cos,1:me%Numvar))/Desnor cosco2 = conqua(cos) - dot_product(me%Conred(cos,1:me%Numact),corvec(1:me%Numact)) cosco2 = cosco2*0.5_wp/Desnor/Desnor write (str,*) "COSCO2/COSCO1=" , cosco2 , cosco1 call me%ogwrit(3,str) if ( cosco1<0.0_wp ) call me%ogwrit(2,str) ! ---------------------------------------------------------------------- foldis = 0.0_wp quacor = cornor*foldis*foldis cosimp = foldis*(cosco1+foldis*cosco2) call me%ogwrit(3,"") write (str,'( "STEEPEST ASCENT FOLLOW",'//' 5X,"DISTANCE",'//& ' 1X,"CORRECTION",'//' 2X,"MERIT_DEL",'//' 6X,"MERIT_VALUE")') call me%ogwrit(3,str) write (str,'("INITIAL.............:",1X,D13.6,'//' 2(1X,D10.3),1X,D16.9)') & foldis , quacor , cosimp ,me%Conval(cos) + cosimp call me%ogwrit(3,str) ! ====================================================================== if ( cosco2<0.0_wp ) then foldis = -0.5_wp*cosco1/cosco2 quacor = cornor*foldis*foldis cosimp = foldis*(cosco1+foldis*cosco2) write (str,'("MERIT MAXIMUM.......:",1X,D13.6,'//'2(1X,D10.3),1X,D16.9)') & foldis , quacor , cosimp , me%Conval(cos) + cosimp call me%ogwrit(3,str) staflg = 1 elseif ( cosco2>0.0_wp ) then foldis = me%Varmax quacor = cornor*foldis*foldis cosimp = foldis*(cosco1+foldis*cosco2) write (str,'("MERIT CONVEX........:",1X,D13.6,'//'2(1X,D10.3),1X,D16.9)') & foldis , quacor , cosimp , me%Conval(cos) + cosimp call me%ogwrit(3,str) staflg = 2 else foldis = me%Varmax quacor = cornor*foldis*foldis cosimp = foldis*(cosco1+foldis*cosco2) write (str,'("MERIT LINEAR........:",1X,D13.6,'//'2(1X,D10.3),1X,D16.9)') & foldis , quacor , cosimp , me%Conval(cos) + cosimp call me%ogwrit(3,str) staflg = 2 end if ! ====================================================================== ! IF MAXIMUM DISTANCE IS HIT ! ---------------------------------------------------------------------- if ( foldis>me%Varmax ) then foldis = me%Varmax quacor = cornor*foldis*foldis cosimp = foldis*(cosco1+foldis*cosco2) write (str,'("MAXIMUM DISTANCE....:",1X,D13.6,'//'2(1X,D10.3),1X,D16.9)') & foldis , quacor , cosimp , me%Conval(cos) + cosimp call me%ogwrit(3,str) staflg = 2 end if ! ====================================================================== ! IF CONVERGENCE ! ---------------------------------------------------------------------- if ( abs(cosimp)<=1.0_wp ) then write (str,'("FINAL...............:",1X,D13.6,'//' 2(1X,D10.3),1X,D16.9,2D11.3)') & foldis , quacor , cosimp , me%Conval(cos) + cosimp , tht , bet call me%ogwrit(2,str) if ( tht/=1.0_wp .or. bet/=0.0_wp ) then tht = 1.0_wp bet = 0.0_wp cycle end if foldis = 0.0_wp Finish = 1 exit main end if ! ====================================================================== ! IF REMAINING DISTANCE IS HIT ! ---------------------------------------------------------------------- if ( foldis>maxdis ) then foldis = maxdis quacor = cornor*foldis*foldis cosimp = foldis*(cosco1+foldis*cosco2) write (str,'("REMAINING DISTANCE..:",1X,D13.6,'//'2(1X,D10.3),1X,D16.9)') & foldis , quacor , cosimp , me%Conval(cos) + cosimp call me%ogwrit(3,str) staflg = 2 end if ! ====================================================================== ! IF CONSTRAINT IS HIT ! ---------------------------------------------------------------------- ind = 0 do con = 1 , me%Numcon if ( me%Contyp(con)==-2 ) cycle if ( me%Conact(con)/=-1 ) cycle co2 = conqua(con) - dot_product(me%Conred(con,1:me%Numact),corvec(1:me%Numact)) co1 = dot_product(me%Conder(con,1:me%Numvar),me%Vardes(1:me%Numvar)) co0 = me%Conval(con)*2.0_wp if ( co2/=0.0_wp ) then det = co1**2 - co2*co0 if ( det<0.0_wp ) cycle det = sqrt(det) val = 1.0e10_wp fac = (-co1+det)/co2 if ( fac>0.0_wp .and. fac<val ) val = fac fac = (-co1-det)/co2 if ( fac>0.0_wp .and. fac<val ) val = fac elseif ( co1/=0.0_wp ) then val = -co0/co1*0.5_wp else cycle end if val = val*Desnor if ( val>0.0_wp .and. val<foldis ) then foldis = val ind = con end if end do ! ---------------------------------------------------------------------- if ( ind/=0 ) then quacor = cornor*foldis*foldis cosimp = foldis*(cosco1+foldis*cosco2) con = ind nam = me%Constr(con) len = me%Conlen(con) write (str,'( "CONSTRAINT REACHED..:",'//'1X,D13.6,2(1X,D10.3),1X,D16.9,1X,I4,1X,A)') & foldis , quacor , cosimp , me%Conval(cos) + cosimp , con , nam(1:len) call me%ogwrit(3,str) staflg = 3 end if ! ====================================================================== ! UPDATE ! ---------------------------------------------------------------------- refdis = foldis exit inner3 end do inner3 inner4: do write (str,'("FINAL...............:",1X,D13.6,'//'2(1X,D10.3),1X,D16.9)') & foldis , quacor , cosimp , me%Conval(cos) + cosimp call me%ogwrit(3,str) ! ---------------------------------------------------------------------- fac = foldis/Desnor ! ---------------------------------------------------------------------- ! VARIABLE DELTA ! ---------------------------------------------------------------------- cosact = me%ogrigt(corvec) dis = 0.0_wp do var = 1 , me%Numvar val = 0.0_wp do act = 1 , me%Numact ind = me%Actcon(act) val = val - cosact(act)*me%Conder(ind,var) end do varvec(var) = fac*(me%Vardes(var)+(val*fac*0.5_wp)) dis = dis + varvec(var)*varvec(var) end do dis = sqrt(dis) ! ---------------------------------------------------------------------- write (str,*) "REFDIS=" , refdis call me%ogwrit(3,str) write (str,*) "FOLDIS=" , foldis call me%ogwrit(3,str) write (str,*) "DIS=" , dis call me%ogwrit(3,str) if ( dis>refdis*1.2_wp .and. me%Senopt>0 ) then faccnt = faccnt + 1 if ( faccnt>=10 ) exit inner4 foldis = foldis*0.5_wp quacor = cornor*foldis*foldis cosimp = foldis*(cosco1+foldis*cosco2) cycle end if ! ---------------------------------------------------------------------- ! UPDATE VARIABLES ! ---------------------------------------------------------------------- Varacc = Varacc + foldis me%Varval = me%Varval + varvec ccc = sqrt(sum((me%Varval-me%Varref)**2)) - me%Varmax**2 if ( ccc>=0.0_wp ) then write (str,*) "CCC > 0" , ccc call me%ogwrit(3,str) staflg = 2 end if select case (staflg) case ( 1 ) ! MAXIMUM REACHED: NEXT ITERATION write (str,'("MERIT MAXIMUM.......:",1X,D13.6,'//'2(1X,D10.3),1X,D16.9,2D11.3)') & foldis , quacor , cosimp , me%Conval(cos) + cosimp , tht , bet call me%ogwrit(2,str) if ( me%Senopt>0 ) Finish = 1 exit inner4 case ( 2 ) ! MAXIMUM TRAVEL DISTANCE REACHED: NEXT ITERATION write (str,'("REMAINING DISTANCE..:",1X,D13.6,'//'2(1X,D10.3),1X,D16.9,2D11.3)') & foldis , quacor , cosimp , me%Conval(cos) + cosimp , tht , bet call me%ogwrit(2,str) exit inner4 case ( 3 ) ! CONSTRAINT HIT: UPDATE CONSTRAINT + CORRECT nam = me%Constr(con) len = me%Conlen(con) write (str,'( "CONSTRAINT REACHED..:",'//'1X,D13.6,2(1X,D10.3),1X,D16.9,2D11.3,1X,I4,1X,A)') & foldis , quacor , cosimp , me%Conval(cos) + cosimp , tht , bet , con , nam(1:len) call me%ogwrit(2,str) convec = conqua - matmul(me%Conred(1:cos,1:me%Numact),corvec(1:me%Numact)) convec = convec*fac*0.5_wp convec = convec + matmul(me%Conder(1:cos,1:me%Numvar),me%Vardes(1:me%Numvar)) me%Conval = me%Conval + convec*fac cycle main end select exit inner4 end do inner4 exit main end do main me%Funvar = desprv me%Confix = me%Conact(1:me%Numcon) if ( me%Senopt==-1 ) me%Senact = me%Conact(1:me%Numcon) deallocate (cosact) deallocate (varvec) deallocate (varwrk) deallocate (corvec) deallocate (desder) deallocate (desprv) deallocate (varprv) deallocate (convec) deallocate (conqua) deallocate (concor) end subroutine ogopti