Merge branch 'dev'

Author: ldeniau

src/mad_tpsa_comp.c

@@ -233,7 +233,7 @@ FUN(compose) (ssz_t sa, const T *ma[sa], ssz_t sb, const T *mb[sb], T *mc[sa])
   printf("mb:\n"); print_damap(sb, mb, 0);
 #endif
 
-  if (hi_ord == 1) compose_ord1(sa,ma, sb,mb, mc);
+  if (hi_ord == 1) compose_ord1(sa,ma, sb,mb, mc_);
 
 #ifdef _OPENMP // TODO: find pcomp heuristic at desc init...
   else if (d->pcomp && hi_ord >= 6 && d->ord2idx[hi_ord+1] >= d->pcomp) {
@@ -249,6 +249,10 @@ FUN(compose) (ssz_t sa, const T *ma[sa], ssz_t sb, const T *mb[sb], T *mc[sa])
 
   else compose(sa,ma, sb,mb, mc_, hi_ord, mo_ord);
 
+#if DEBUG_COMPOSE
+  printf("mc:\n"); print_damap(sa, TC mc_, 0);
+#endif
+
   // copy back
   FOR(ia,sa) if (amc[ia]) {
     FUN(copy)(mc_[ia], mc[ia]);

src/madl_damap.mad

@@ -669,105 +669,6 @@ end
 
 -- map logarithm of exponential poisson bracket: r = exp(:log(x):) y
 
-local function fgrad (x, t)
-  local r = t:same()
-  for i=1,#x do r = r + x[i]*t:deriv(i) end
-  return r
-end
-
-local function liebra (x, y)
-  local r = x:same()
-  for i=1,#x do
-  for j=1,#x do
-    local t1 = y[j] * x[i]:deriv(j)
-    local t2 = x[j] * y[i]:deriv(j)
-    r[i] = r[i] + t2 - t1
-  end end
-  return r
-end
-
-local function mnrm (x)
-  local nrm = 0
-  for i=1,#x do nrm = nrm + x[i]:nrm() end
-  return nrm
-end
-
-local function logpb (x, r, eps) -- c_logf + c_flofacg
-  -- c_logf
-  local x0 = x:copy():set0(0)
-
-local fid = 1262 -- ..1289
-x0:write("fort/fort_n."..fid, "x0") ; fid = fid+1 ; -- 1263
-
-  -- c_flofacg
-  local v = x:same():eye() -- identity
-
-x0:write("fort/fort_n."..fid, "x") ; fid = fid+1 ; -- 1264
-
-  -- special case where eps = -no
-  if eps and eps <= 0 then
-    for k=1,round(eps)-1 do
-      local y = (-r):exppb(x0)
-      local t = y - v
-      r = r+t
-
-t:write("fort/fort_n."..fid, "t_0 @ "..k) ; fid = fid+1 ;
-r:write("fort/fort_n."..fid, "r_0 @ "..k) ; fid = fid+1 ;
-    end
-    return r
-  end
-
-  local nrm0 = mnrm(x0)
-  local epsone = eps or nrm0/1000
-  local xn, xn0 = 1e4, 1e36
-  local nrmax = 1000
-  local conv = false
-
-  for k=1,nrmax do
-    local y = (-r):exppb(x0)
-    local t = y - v
-
-t:write("fort/fort_n."..fid, "t_1 @ "..k) ; fid = fid+1 ; -- 1265
-
-    if xn < epsone then -- CBH (* == c_bra_v_ct)
-      local e2, e3 = x:same(), x:same()
-      for i=1,#x do e2[i] = -0.5*fgrad(t ,t[i]) end
-      for i=1,#x do e3[i] = -0.5*fgrad(e2,t[i])-(1/6)*fgrad(t,e2[i]) end
-      t = t+e2+e3
-
-t:write("fort/fort_n."..fid, "t_2 @ "..k) ; fid = fid+1 ; -- 1266
-
-      local z  = liebra(r,t) -- <r,t>
-      local e2 = liebra(r,z) -- <r,<r,t>>
-      local e3 = liebra(t,z) -- <t,<r,t>>
-      t = t + 0.5*z + (1/12)*(e2-e3)
-
-t:write("fort/fort_n."..fid, "t_3 @ "..k) ; fid = fid+1 ; -- 1267
-
-      local e4 = liebra(t,e2) -- <t,<r,<r,t>>>
-      t = t - (1/24)*e4
-
-t:write("fort/fort_n."..fid, "t_4 @ "..k) ; fid = fid+1 ; -- 1268
-
-    end
-    r = r+t
-
-r:write("fort/fort_n."..fid, "r_1 @ "..k) ; fid = fid+1 ;
-
-    local nrm = mnrm(t)
-    xn = nrm/nrm0
-
-print("xn=", xn, "epsone=", epsone, "nrm=", nrm, "nrm0=", nrm0)
-
-    -- same scheme as exppb
-    if xn < 1e-14 or conv and xn >= xn0 then break end
-    if xn < 1e-10 then conv = true end
-    xn0 = xn
-  end
-
-  return r
-end
-
 function MR.logpb (x, y0, r)
   if is_string(r) and r == 'in' then r = x end -- in place
   r = r or map_alloc(x)
@@ -1259,7 +1160,6 @@ function MR.__pow (x, n, r)
     if n == 0 then break end
     t:compose(t,t)                          -- t = t^2
   end
-
   return r
 end