' ******************************************************************
' Minimization of a sum of squared functions along a line
' (Used internally by equation solvers)
' ******************************************************************

' ------------------------------------------------------------------
' Error codes
' ------------------------------------------------------------------

#DEFINE OptOk 0
' No error

#DEFINE OptNonConv -1
' Non-convergence

#DEFINE MatErrDim -3
' Non-compatible dimensions

' ------------------------------------------------------------------
' Mathematical constants
' ------------------------------------------------------------------

#DEFINE Gold 1.618033988749895#   
' Golden Mean (1 + Sqrt(5))/2

#DEFINE CGold 0.3819660112501052#  
' 2 - Gold

' ------------------------------------------------------------------
' Global variable
' ------------------------------------------------------------------

COMMON SHARED ErrCode AS INTEGER
' Error code from the latest function evaluation

' ******************************************************************

FUNCTION SumSqrFn(Equations AS SUB(X() AS DOUBLE, F() AS DOUBLE),_
                  X()       AS DOUBLE,_ 
                  F()       AS DOUBLE) AS DOUBLE
' ----------------------------------------------------------------------
' Computes the sum of squared functions F(i)^2 at point X
' ----------------------------------------------------------------------

  DIM Sum AS DOUBLE, I AS INTEGER

  Equations(X(), F())

  Sum = 0
  FOR I = LBOUND(F) TO UBOUND(F)
    Sum = Sum + F(I) * F(I)
  NEXT I

  SumSqrFn = Sum
END FUNCTION

SUB LinMinEq(      Equations AS SUB(X() AS DOUBLE, F() AS DOUBLE),_
                   X()       AS DOUBLE,_ 
                   DeltaX()  AS DOUBLE,_ 
             BYREF R         AS DOUBLE,_ 
                   MaxIter   AS INTEGER,_ 
                   Tol       AS DOUBLE,_ 
                   F()       AS DOUBLE)
' ----------------------------------------------------------------------
' Minimizes a sum of squared functions from point X in the direction 
' specified by DeltaX, using golden search as the minimization algorithm
' ----------------------------------------------------------------------
' Input parameters  : Equations = system of equations
'                     X         = starting point
'                     DeltaX    = search direction
'                     R         = initial step, in fraction of |DeltaX|
'                     MaxIter   = maximum number of iterations
'                     Tol       = required precision
' ----------------------------------------------------------------------
' Output parameters : X = refined minimum coordinates
'                     F = function values at minimum
'                     R = step corresponding to the minimum
' ----------------------------------------------------------------------
' Possible results  : OptOk      = no error
'                     OptNonConv = non-convergence
' ----------------------------------------------------------------------

  CONST MinTol = 1.5D-8  ' Minimal tolerance (~ Sqrt(MachEp))
 
  DIM AS INTEGER I, L, N, Iter
  DIM AS DOUBLE  A, B, C, Fa, Fb, Fc, R0, R1, R2, R3, F1, F2, Norm

  L = LBOUND(X)
  N = UBOUND(X)

  IF LBOUND(DeltaX) <> L OR LBOUND(F) <> L OR _
     UBOUND(DeltaX) <> N OR UBOUND(F) <> N THEN
       ErrCode = MatErrDim
       EXIT SUB
  END IF

  DIM P(L TO N) AS DOUBLE
 
  IF R <= 0 THEN R = 1

  Norm = 0
  FOR I = L TO N
    Norm = Norm + DeltaX(I) * DeltaX(I)
  NEXT I
  Norm = SQR(Norm)

  ' Bracket the minimum (see procedure MinBrack in module GOLDSRCH)

  A = 0: B = R * Norm
 
  FOR I = L TO N
    P(I) = X(I)
  NEXT I
 
  Fa = SumSqrFn(Equations, P(), F())

  FOR I = L TO N
    P(I) = X(I) + B * DeltaX(I)
  NEXT I

  Fb = SumSqrFn(Equations, P(), F())

  IF Fb > Fa THEN
    SWAP A, B
    SWAP Fa, Fb
  END IF

  C = B + Gold * (B - A)

  FOR I = L TO N
    P(I) = X(I) + C * DeltaX(I)
  NEXT I

  Fc = SumSqrFn(Equations, P(), F())

  DO WHILE Fc < Fb
    A = B
    B = C
    Fa = Fb
    Fb = Fc
    C = B + Gold * (B - A)

    FOR I = L TO N
      P(I) = X(I) + C * DeltaX(I)
    NEXT I

    Fc = SumSqrFn(Equations, P(), F())
  LOOP

  IF A > C THEN
    SWAP A, C
    SWAP Fa, Fc
  END IF

  ' Refine the minimum (see procedure GoldSearch in module GOLDSRCH)

  IF Tol < MinTol THEN Tol = MinTol
 
  R0 = A: R3 = C

  IF (C - B) > (B - A) THEN
    R1 = B
    R2 = B + CGold * (C - B)
    F1 = Fb
  
    FOR I = L TO N
      P(I) = X(I) + R2 * DeltaX(I)
    NEXT I
  
    F2 = SumSqrFn(Equations, P(), F())
  ELSE
    R1 = B - CGold * (B - A)
    R2 = B

    FOR I = L TO N
      P(I) = X(I) + R1 * DeltaX(I)
    NEXT I

    F1 = SumSqrFn(Equations, P(), F())
    F2 = Fb
  END IF

  Iter = 0

  DO WHILE (Iter <= MaxIter) AND (ABS(R3 - R0) > Tol * (ABS(R1) + ABS(R2)))
    IF F2 < F1 THEN
      R0 = R1
      R1 = R2
      F1 = F2
      R2 = R1 + CGold * (R3 - R1)

      FOR I = L TO N
        P(I) = X(I) + R2 * DeltaX(I)
      NEXT I
    
      F2 = SumSqrFn(Equations, P(), F())
    ELSE
      R3 = R2
      R2 = R1
      F2 = F1
      R1 = R2 - CGold * (R2 - R0)
    
      FOR I = L TO N
        P(I) = X(I) + R1 * DeltaX(I)
      NEXT I
    
      F1 = SumSqrFn(Equations, P(), F())
    END IF

    Iter = Iter + 1
  LOOP

  IF F1 < F2 THEN R = R1 ELSE R = R2
  
  FOR I = L TO N
    X(I) = X(I) + R * DeltaX(I)
  NEXT I

  Equations(X(), F())
  
  IF Iter > MaxIter THEN
    ErrCode = OptNonConv
  ELSE
    ErrCode = OptOk
  END IF
END SUB

