?? asa_cg.cpp
字號:
if ( phi <= fquad ) { if ( Parm->PrintLevel >= 4 ) { printf ("alpha: %14.6e phi: %14.6e fquad: %14.6e\n", alpha, phi, fquad) ; } if ( asa_Wolfe (alpha, phi, dphi, Com) ) { status = 0 ; goto Exit ; } } } if ( phi <= Com->fpert ) { a = alpha ; dphia = dphi ; } else { /* contraction phase, only break at termination or Secant step */ b = alpha ; while ( TRUE ) { alpha = .5*(a+b) ; nshrink++ ; if ( nshrink > Parm->nexpand ) { status = 6 ; goto Exit ; } asa_step (xtemp, x, d, alpha, nfree) ; asa_g (gtemp, xtemp, Com) ; dphi = asa_dot (gtemp, d, nfree) ; if ( dphi >= ZERO ) goto Secant ; phi = asa_f (xtemp, Com) ; if ( Parm->PrintLevel >= 4 ) { printf ("contract, a: %14.6e b: %14.6e alpha: %14.6e phi: " "%14.6e dphi: %14.6e\n", a, b, alpha, phi, dphi) ; } if ( Com->QuadOK && (phi <= fquad) ) { if ( asa_Wolfe (alpha, phi, dphi, Com) ) { status = 0 ; goto Exit ; } } if ( phi <= Com->fpert ) { a = alpha ; dphia = dphi ; } else { b = alpha ; } } }/* expansion phase */ ngrow++ ; if ( ngrow > Parm->nexpand ) { status = 3 ; goto Exit ; } alphaold = alpha ; alpha = MIN (rho*alpha, minstep) ; if ( alpha != alphaold ) { asa_step (xtemp, x, d, alpha, nfree) ; asa_g (gtemp, xtemp, Com) ; dphi = asa_dot (gtemp, d, nfree) ; if ( Parm->PrintLevel >= 4 ) { printf ("expand, a: %14.6e alpha: %14.6e phi: " "%14.6e dphi: %14.6e\n", a, alpha, phi, dphi) ; } } else /* a new constraint is active */ { do /* while statement */ { alphaold = alpha ; phiold = phi ; if ( alpha < maxstep ) { alpha = rho*alphaold ; asa_project (xtemp, x, d, alpha, Com) ; phi = asa_f (xtemp, Com) ; } } while ( phi < phiold ) ; if ( alphaold == minstep ) { asa_step (xtemp, x, d, minstep, nfree) ; status = -3 ; } else { asa_project (xtemp, x, d, alphaold, Com) ; asa_g (gtemp, xtemp, Com) ; status = -4 ; } phi = phiold ; goto Exit ; } }Secant: b = alpha ; dphib = dphi ; if ( Com->QuadOK ) { phi = asa_f (xtemp, Com) ; if ( ngrow + nshrink == 0 ) fquad = MIN (phi, Com->f0) ; if ( phi <= fquad ) { if ( asa_Wolfe (alpha, phi, dphi, Com) ) { status = 0 ; goto Exit ; } } } nsecant = Parm->nsecant ; for (iter = 1; iter <= nsecant; iter++) { if ( Parm->PrintLevel >= 4 ) { printf ("secant, a: %14.6e b: %14.6e da: %14.6e db: %14.6e\n", a, b, dphia, dphib) ; } width = Parm->gamma*(b - a) ; if ( -dphia <= dphib ) alpha = a - (a-b)*(dphia/(dphia-dphib)) ; else alpha = b - (a-b)*(dphib/(dphia-dphib)) ; c = alpha ; a0 = a ; b0 = b ; da0 = dphia ; db0 = dphib ; status = asa_update (&a, &dphia, &b, &dphib, &alpha, &phi, &dphi, Com) ; if ( status >= 0 ) goto Exit ; else if ( status == -2 ) { if ( c == a ) { if ( dphi > da0 ) alpha = c - (c-a0)*(dphi/(dphi-da0)) ; else alpha = a ; } else { if ( dphi < db0 ) alpha = c - (c-b0)*(dphi/(dphi-db0)) ; else alpha = b ; } if ( (alpha > a) && (alpha < b) ) { if ( Parm->PrintLevel >= 4 ) printf ("2nd secant\n") ; status = asa_update (&a, &dphia, &b, &dphib, &alpha, &phi, &dphi, Com) ; if ( status >= 0 ) goto Exit ; } }/* bisection iteration */ if ( b-a >= width ) { alpha = .5*(b+a) ; if ( Parm->PrintLevel >= 4 ) printf ("bisection\n") ; status = asa_update (&a, &dphia, &b, &dphib, &alpha, &phi, &dphi, Com) ; if ( status >= 0 ) goto Exit ; } else if ( b <= a ) { status = 7 ; goto Exit ; } } status = 4 ;Exit: Com->alpha = alpha ; Com->f = phi ; Com->df = dphi ; return (status) ;}/* ========================================================================= === asa_lineW =========================================================== ========================================================================= Ordinary Wolfe line search routine. This routine is identical to asa_line except that the function psi [a] = phi [a] - phi [0] - a*delta*dphi [0] is minimized instead of the function phi ========================================================================= */int asa_lineW( double dphi0 , /* function derivative at starting point (alpha = 0) */ asa_com *Com /* cg com structure */){ int i, iter, nfree, nsecant, nshrink, ngrow, status ; double a, dpsia, b, dpsib, c, alpha, phi, dphi, alphaold, phiold, a0, da0, b0, db0, width, fquad, rho, psi, dpsi, minstep, maxstep, *x, *d, *xtemp, *gtemp ; asacg_parm *Parm ; nfree = Com->nfree ; x = Com->x ; /* current iterate */ d = Com->d ; /* current search direction */ xtemp = Com->xtemp ; /* x + alpha*d */ gtemp = Com->gtemp ; /* gradient at x + alpha*d */ minstep = Com->minstep ; maxstep = Com->maxstep ; alpha = Com->alpha ; if ( alpha > minstep ) { alpha = minstep ; Com->QuadOK = FALSE ; } phi = Com->f ; Parm = Com->cgParm ; rho = Parm->rho ; asa_step (xtemp, x, d, alpha, nfree) ; asa_g (gtemp, xtemp, Com) ; dphi = asa_dot (gtemp, d, nfree) ; /* check if gradient is nan; if so, reduce stepsize */ if ( dphi != dphi ) { for (i = 0; i < Parm->nexpand; i++) { alpha /= rho ; asa_step (xtemp, x, d, alpha, nfree) ; asa_g (gtemp, xtemp, Com) ; dphi = asa_dot (gtemp, d, nfree) ; if ( dphi == dphi ) break ; } if ( i == Parm->nexpand ) { status = -2 ; goto Exit ; } Com->QuadOK = FALSE ; rho = Parm->nan_rho ; } dpsi = dphi - Com->wolfe_hi ; /*Find initial interval [a,b] such that dphia < 0, dphib >= 0, and phia <= phi0 + feps*Ck */ a = ZERO ; dpsia = dphi0 - Com->wolfe_hi ; ngrow = 0 ; nshrink = 0 ; while ( dpsi < ZERO ) { phi = asa_f (xtemp, Com) ; psi = phi - alpha*Com->wolfe_hi ; /* if quadstep in effect and quadratic conditions hold, check Wolfe condition*/ if ( Com->QuadOK ) { if ( ngrow == 0 ) fquad = MIN (phi, Com->f0) ; if ( phi <= fquad ) { if ( Parm->PrintLevel >= 4 ) { printf ("alpha: %14.6e phi: %14.6e fquad: %14.6e\n", alpha, phi, fquad) ; } if ( asa_Wolfe (alpha, phi, dphi, Com) ) { status = 0 ; goto Exit ; } } } if ( psi <= Com->fpert ) { a = alpha ; dpsia = dphi ; } else { /* contraction phase, only break at termination or Secant step */ b = alpha ; while ( TRUE ) { alpha = .5*(a+b) ; nshrink++ ; if ( nshrink > Parm->nexpand ) { status = 6 ; goto Exit ; } asa_step (xtemp, x, d, alpha, nfree) ; asa_g (gtemp, xtemp, Com) ; dphi = asa_dot (gtemp, d, nfree) ; dpsi = dphi - Com->wolfe_hi ; if ( dpsi >= ZERO ) goto Secant ; phi = asa_f (xtemp, Com) ; psi = phi - alpha*Com->wolfe_hi ; if ( Parm->PrintLevel >= 4 ) { printf ("contract, a: %14.6e b: %14.6e alpha: %14.6e phi: " "%14.6e dphi: %14.6e\n", a, b, alpha, phi, dphi) ; } if ( Com->QuadOK && (phi <= fquad) ) { if ( asa_Wolfe (alpha, phi, dphi, Com) ) { status = 0 ; goto Exit ; } } if ( psi <= Com->fpert ) { a = alpha ; dpsia = dpsi ; } else { b = alpha ; } } }/* expansion phase */ ngrow++ ; if ( ngrow > Parm->nexpand ) { status = 3 ; goto Exit ; } alphaold = alpha ; alpha = MIN (rho*alpha, minstep) ; if ( alpha != alphaold ) { asa_step (xtemp, x, d, alpha, nfree) ; asa_g (gtemp, xtemp, Com) ; dphi = asa_dot (gtemp, d, nfree) ; dpsi = dphi - Com->wolfe_hi ; if ( Parm->PrintLevel >= 4 ) { printf ("expand, a: %14.6e alpha: %14.6e phi: " "%14.6e dphi: %14.6e\n", a, alpha, phi, dphi) ; } } else /* a new constraint is active */ { do /* while statement */ { alphaold = alpha ; phiold = phi ; if ( alpha < maxstep ) { alpha = rho*alphaold ; asa_project (xtemp, x, d, alpha, Com) ; phi = asa_f (xtemp, Com) ; } } while ( phi < phiold ) ; if ( alphaold == minstep ) { asa_step (xtemp, x, d, minstep, nfree) ; status = -3 ; } else { asa_project (xtemp, x, d, alphaold, Com) ; asa_g (gtemp, xtemp, Com) ; status = -4 ; } phi = phiold ; goto Exit ; ?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -