ghc/tests/programs/cholewo-eval/Main.lhs

   1 \begin{code}
   2 module Main(main) where
   3 import Arr
   4 \end{code}
   5
   6 \begin{code}
   7 type F a = Vector a -> a
   8 type DF a = Vector a -> Vector a
   9 \end{code}
  10
  11 \begin{code}
  12 data {-(Eval a) =>-} ScgData a = ScgData {k :: !Int, err :: !a,
  13                           w, p, r :: !(Vector a),
  14                           delta, pnorm2, lambda, lambdabar :: !a,
  15                           success :: !Bool}
  16 \end{code}
  17
  18 \begin{code}
  19 calculate2order :: Floating a => ScgData a -> a -> DF a -> ScgData a
  20 calculate2order d sigma1 df =
  21   let pnorm2' = vectorNorm2 (p d)
  22       sigma = sigma1 / (sqrt pnorm2')
  23       s = scaleVector (1/sigma) (df ((w d) + (scaleVector sigma (p d))) - df (w d))
  24   in d {pnorm2 = pnorm2', delta = inner (p d) s}
  25 \end{code}
  26
  27 \begin{code}
  28 hessPosDef :: (Floating a, Ord a) => ScgData a -> ScgData a
  29 hessPosDef d =
  30   let delta' = delta d + (lambda d - lambdabar d) * pnorm2 d {- step 3 -}
  31   in if delta' <= 0                             {- step 4 -}
  32      then let lambdabar' = 2.0 * (lambda d - delta' / pnorm2 d)
  33           in d {delta = -delta' + lambda d * pnorm2 d, lambda = lambdabar', lambdabar = lambdabar'}
  34      else d {delta = delta'}
  35 \end{code}
  36
  37 \begin{code}
  38 reduceError :: (Floating a, Ord a) => ScgData a -> DF a -> Bool -> a -> a -> ScgData a
  39 reduceError d df restart bdelta mu =
  40   let r' = negate (df (w d))
  41       p' = if restart
  42            then r'
  43            else let beta = (vectorNorm2 r' - inner r' (r d)) / mu
  44                 in r' + scaleVector beta (p d)
  45   in d {p = p', r = r', lambda = if bdelta >= 0.75 then lambda d / 4 else lambda d
  46     }
  47 \end{code}
  48
  49 \begin{code}
  50 data ScgInput a = ScgInput (F a) (DF a) (Vector a)
  51 \end{code}
  52
  53 \begin{code}
  54 scgIter :: (Floating a, Ord a) => ScgInput a -> [ScgData a]
  55 scgIter (ScgInput f df w1) =
  56     let p1 = negate (df w1)                     {- step 1 -}
  57         r1 = p1
  58         pnorm21 = vectorNorm2 p1
  59         n = vectorSize w1
  60         sigma1 = 1.0e-4
  61         lambda1 = 1.0e-6
  62         err1 = f w1
  63     in iterate (\d ->
  64            let d1 = if success d                {- step 2 -}
  65                     then calculate2order d sigma1 df
  66                     else d
  67                d2 = hessPosDef d1
  68                mu = inner (p d2) (r d2)         {- step 5 -}
  69                alpha = mu / (delta d2)
  70                w' = (w d2) + scaleVector alpha (p d2)
  71                err' = f w'
  72                bdelta = 2 * (delta d2) * ((err d2) - err') / (mu^2) {- step 6 -}
  73                success' = (bdelta >= 0)         {- step 7 -}
  74                restart = ((k d) `mod` n == 0)
  75                d3 = if success'
  76                     then (reduceError (d2 {w = w'}) df restart bdelta mu)
  77                             {err = err', lambdabar = 0}
  78                     else d2 {lambdabar = lambda d2}
  79                d4 = if bdelta < 0.25          {- step 8 -}
  80                     then d3 {lambda = (lambda d3) + (delta d3) * (1 - bdelta) / (pnorm2 d3)}
  81                     else d3
  82            in d4 {k = k d4 + 1, success = success'}
  83        )
  84        (ScgData 1 err1 w1 p1 r1 0.0 pnorm21 lambda1 0.0 True)
  85 \end{code}
  86
  87 \begin{code}
  88 rosenbrock = ScgInput
  89   (\ (Vector x) -> 100 * (x!2 - x!1^2)^2 + (1 - x!1)^2)
  90   (\ (Vector x) -> listVector [-2 * (1 - x!1) - 400 * x!1 * (x!2 - x!1^2),
  91                               200 * (x!2 -x!1^2)])
  92   (listVector [-1,-1.0])
  93 \end{code}
  94
  95
  96 \begin{code}
  97 main = print (w ((scgIter rosenbrock)!!1))
  98 \end{code}