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 = 4; 
optu = 1; 
lept = 0; 
leph = 0; 
leet = 0; 
lewx = 0; 
    
ley = 0; 
lefx = 0; 
lejo1 = 0; 
ledx = 0; 
leer = 0; 
lejo2 = 0; 
    
lejo3 = 0; 
leno = 0; 
\)\)

User Settings for RHSE

Input Conditions

<< seinput.m

Setting Practice Conditions

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

Generate Vectors

\(lmdv = Array[lmd, dimx]; 
pv = Array[p, dimp]; 
wv = Array[w, dimw]; 
    
yv = Array[y, dimy]; 
dyv = Array[dy, dimy]; 
dhv = Array[dh, dimx]; 
    \)

Generate RHSE Algorithm

Euler-Lagrange Equations

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

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

    wsol = Solve[
        Table[Hw[[i]] == 0, {i, dimw}], 
        Table[wv[[i]], {i, dimw}]]; 

    wxl = Array[w, dimw] /. wsol[[1]]; 
    \)

Boudary Conditions for Backward Sweep

\(phix = Dv[phi, xv]; 
If[smpu != 0, \(phix = Simplify[phix]; \)]; 

    phix2 = Dv[phix, xv]; 
invphix2 = Inverse[phix2]; 

    dsyminvphix2 = DeSym[invphix2, dimx]; 

    If[smpu != 0, \(dsyminvphix2 = Simplify[dsyminvphix2]; \)]; 

    invphix2 = Sym[dsyminvphix2, dimx]; 
phixy = Dv[phix, yv]; 

    If[smpu != 0, \(phixy = Simplify[phixy]; \)]; 
phixhx = Dv[phix, hxv]; 

    If[smpu != 0, \(phixhx = Simplify[phixhx]; \)]; 
etax = Dv[eta, xv]; 

    If[smpu != 0, \(etax = Simplify[etax]; \)]; 
etax2 = Dv[etax, xv]; 

    dsymetax2 = DeSym[etax2, dimx]; 

    If[smpu != 0, \(dsymetax2 = Simplify[dsymetax2]; \)]; 

    etaxy = Dv[etax, yv]; 
If[smpu != 0, \(etaxy = Simplify[etaxy]; \)]; 

    erropt = lmdv + phix; 
zetatemp = zeta; 
\(Clear[zeta]; \)\)\(\(If[dimp > 0, phixp = Dv[phix, pv]; 
	etaxp = Dv[etax, pv]; 
	
      If[smu != 0, phixp = Simplify[phixp]; 
		etaxp = Simplify[etaxp]; ]; 
      cf = invphix2 . \ 
          \((\(-\((phix2 . fxu + phixp . ∂\_t pt + phixy . dyv + 
                      phixhx . dhv + dlmddt)\)\)\ onepdtfdt - zeta\ erropt)
            \), \(cf = 
        invphix2 . 
          \((\(-\((phix2 . fxu + dlmddt)\)\)\ onepdtfdt - zeta\ erropt)\); 
      \)]; \)\)

Differential Equations for Bacward Sweep

\(fx = Dv[fxu, xv]; 
fw = Dv[fxu, wv]; \)

\(Hx = Dv[H, xv]; 
Hw = Dv[H, wv]; 
Hxx = Dv[Hx, xv]; 

    Hxw = Dv[Hx, wv]; 
Hwx = Transpose[Hxw]; 
Hww = Dv[Hw, wv]; 

    InvHww = Inverse[Hww]; \)

\(As = fx - fw . InvHww . Hwx; 
Bs = fw . InvHww . Transpose[fw]; 

    Cs = Hxx - Hxw . InvHww . Hwx; \)

\(\(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 = As . S + S . Transpose[As] + S . Cs . S - Bs; 

    cdot = \((As + S . Cs)\) . cv; \)

\(dsdt = DeSym[Sdot, dimx]; 

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

dsdt >> "sdot.tmp"cdot >> "cdot.tmp"

\(dsdt = << "sdot.tmp"; 
cdot = << "cdot.tmp"; \)

Differential Equations for State Estimation

\(\(Id = IdentityMatrix[dimx]; 
invm = Inverse[Id - S . etax2]; 

    If[smpu >= 4, \(invm = Simplify[invm]; \)]; 

    dxdt = invm . 
        \((fxu + 
            S . \((Hx + etaxp . ∂\_t pt + etaxy . dyv + etax2 . dhv)
                \) + cv)\); 
If[smpu >= 4, \(dxdt = Simplify[dxdt]; \)]; 

    \)\)

Generate C Code

Load Packages

<< FormatFolder`Format2`

<< FormatFolder`Optimize2`

Setting Options

SetOptions[Optimize, OptimizeVariable -> {o, Array}, 
      OptimizeNull -> {List, x, u, y, dy, dh, w, lmd, p, s, c, o, hx}, 
      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]\)]]; 

      leet = Length[\(Optimize[etax]\)[\([1]\)]]; 

      lewx = Length[\(Optimize[wxl]\)[\([1]\)]]; 

      lefx = Length[\(Optimize[fxu]\)[\([1]\)]]; 

      ley = Length[\(Optimize[fy]\)[\([1]\)]]; 

      ledx = Length[\(Optimize[dxdt]\)[\([1]\)]]; 	

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

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

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

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

      leno = Max[{lept, leph, leet, lewx, lefx, ley, 
			\ \ ledx, lejo1, 
            leer, lejo2, lejo3}]; ]; \)\)

Conversion of Arrays

\(ArrNames = 
      Join[MyArrNames, {"x", "hx", "w", "lmd", "p", 
          "s", "c", "y", "dy", "dh"}]; 

    ArrDims = 
      Join[MyArrDims, {dimx, dimx, dimw, dimx, dimp, 
          1\/2\ dimx\ \((dimx + 1)\), dimx, dimy, dimy, 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 C Code

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

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