Functions for C Code Generation

Differentiation by a Vector:  Dv

Dv[f_, x_] := 
      Module[{j, ans}, 
        ans = Table[
            (∂_(x[[j]]))f, {
              j, Length[x]}]; 
        If[Length[Dimensions[ans]] ≠ 1, ans = Transpose[ans]]; ans]

Constuction of an n by n Symmetric Matrix from a Vector

Sym[v_, n_] := 
      Module[{i, j, iv, ans}, ans = Table[0, {i, n}, {j, n}]; iv = 1; i = 0; 
        Do[Do[ans[[j, j + i]] = 
              v[[iv]]; 
            ans[[j + i, j]] = 
              v[[iv]]; 
            iv++, {j, n - i}], {i, 0, n - 1}]; ans]

Construction of a Vector from an n by n Symmetric Matrix

DeSym[S_, n_] := 
      Module[{i, j, ans = {}}, 
        Do[Do[AppendTo[ans, 
              S[[j, j + i]]], {j, 
              n - i}], {i, 0, n - 1}]; ans]

Initial Value

\(\(smpu = 3; 
optu = 1; 
lept = 0; 
leph = 0; 
leux = 0; 

    lejo1 = 0; 
leer = 0; 
lejo2 = 0; 
lejo3 = 0; 
leno = 0; 
\)\)

User Settings for RHC

Input Conditions

<< input.m

Setting Practice Conditions

If[optu != 0, optuse = "True", optuse = "False"];

Generate Vectors

\(xv = Array[x, dimx]; 
lmdv = Array[lmd, dimx]; 
uv = Array[u, dimu]; 
    
pv = Array[p, dimp]; \)

Generate RHC Algorithm

Euler-Lagrange Equations

H = L + lmdv . fxu; \(dlmddt = \(-Dv[H, xv]\); 

    If[smpu != 0, \(dlmddt = Simplify[dlmddt]; \)]; 
Hu = Dv[H, uv]; 

    \(usol = Solve[
        Table[Hu[[i]] == 0, {i, dimu}], 
        Table[uv[[i]], {i, dimu}]]; 
    \)\)uxl = 
      Array[u, dimu] /. usol[[1]];

Boudary Conditions for Backward Sweep

phix = Dv[phi, xv]; phix = Simplify[phix]; \(phix2 = Dv[phix, xv]; 
dsymphix2 = DeSym[phix2, dimx]; 

    If[smpu != 0, \(dsymphix2 = Simplify[dsymphix2]; \)]\)\(erropt = lmdv - phix; 
\(zetatemp = zeta; \)\)Clear[zeta]; \(\(If[dimp > 0, phixp = Dv[phix, pv]; 
	
      If[smu != 0, \(phixp = Simplify[phixp]; \)]; 
      cf = \((phix2 . fxu + phixp . ∂\_t pt - dlmddt)\)\ 
            onepdtfdt - zeta\ erropt, 
      \(cf = \((phix2 . fxu - dlmddt)\)\ onepdtfdt - zeta\ erropt; \)]; 
    \)\)

Differential Equations for Bacward Sweep

fx = Dv[fxu, xv]; fu = Dv[fxu, uv];

Hx = Dv[H, xv]; Hu = Dv[H, uv]; Hxx = Dv[Hx, xv]; Hxu = Dv[Hx, uv]; Hux = Transpose[Hxu]; Huu = Dv[Hu, uv]; InvHuu = Inverse[Huu];

\(As = fx - fu . InvHuu . Hux; 
Bs = fu . InvHuu . Transpose[fu]; 

    Cs = Hxx - Hxu . InvHuu . Hux; \)

\(\(If[smpu >= 2, 
As = Simplify[Expand[As]]; 

      Bs = Simplify[Expand[Bs]]; 
Cs = Simplify[Expand[Cs]]; ]; \)\)

S = Sym[Array[s, (1/2) dimx (dimx + 1)], dimx]; cv = Array[c, dimx];

Sdot = -Transpose[As] . S - S . As + S . Bs . S - Cs; cdot = -(Transpose[As] - S . Bs) . cv;

\(dsdt = DeSym[Sdot, dimx]; 

    If[smpu >= 3, 
dsdt = Simplify[dsdt]; 
cdot = Simplify[cdot]; ]; \)

dsdt >> "sdot1.tmp"cdot >> "cdot1.tmp"

\(dsdt = << "sdot1.tmp"; 
cdot = << "cdot1.tmp"; \)

Generate C Code

Load Packages

<< FormatFolder`Format2`

<< FormatFolder`Optimize2`

Setting Options

SetOptions[Optimize, OptimizeVariable -> {o, Array}, 
      OptimizeNull -> {List, x, u, lmd, p, s, c, o}, 
      OptimizeCoefficients -> True, OptimizeProcedure -> False];

\(\(SetOptions[CAssign, AssignOptimize -> ToExpression[optuse], 
      AssignReplace -> {\ "" -> ""}, AssignToArray -> {""}, 
      AssignPrecision -> 2, AssignBreak -> {75, "\
"
		}]; 
    \)\)

Optimize Dimensions

\(\(If[optu == 1, lept = Length[\(Optimize[pt]\)[\([1]\)]]; 

      leph = Length[\(Optimize[phix]\)[\([1]\)]]; 

      leux = Length[\(Optimize[uxl]\)[\([1]\)]]; 

      lejo1 = Length[\(Optimize[Join[fxu, dlmddt]]\)[\([1]\)]]; 

      leer = Length[\(Optimize[erropt]\)[\([1]\)]]; 

      lejo2 = Length[\(Optimize[Join[dsymphix2, cf]]\)[\([1]\)]]; 

      lejo3 = Length[\(Optimize[Join[dsdt, cdot]]\)[\([1]\)]]; 

      leno = Max[{lept, leph, leux, lejo1, leer, lejo2, lejo3}]; ]; \)\)

Conversion of Arrays

\(ArrNames = 
      Join[MyArrNames, {"x", "u", "lmd", "p", "s", 
          "c"}]; 

    ArrDims = 
      Join[MyArrDims, {dimx, dimu, dimx, dimp, 1\/2\ dimx\ \((dimx + 1)\), 
          dimx}]; 
ArrCnv = {}; 

    \(Do[ArrCnv = 
        Join[ArrCnv, 
          Table[ArrNames[[j]]<>"["<>
                ToString[i]<>"]" -> 
              ArrNames[[j]]<>"["<>
                ToString[i - 1]<>"]", {i, 
              ArrDims[[j]]}]], {j, 
        Length[ArrNames]}]; 
ArrCnv2 = ArrCnv; 

    Do[
	ArrCnv2 = 
        Join[ArrCnv2, 
          Table["o["<>ToString[i]<>"]" -> 
              "o["<>ToString[i - 1]<>"]", {i, leno}]]
]\)\)

Output CProgram

\(outfile = StringJoin[fnout, ".c"]; 

    Splice["AutoGen.mc", outfile, FormatType -> OutputForm]; \)