Fixed SBGMap preImage

sbg/lexp.cpp

@@ -125,6 +125,7 @@ LExp LExp::inverse() const
   return LExp(new_slope, new_offset);
 }
 
+
 LExp LExp::mod(const LExp &other) const
 {
   Util::ERROR_UNLESS(other.slope_ == 0

sbg/map.cpp

@@ -27,12 +27,16 @@ namespace LIB {
 
 void compatible(Interval i, LExp le)
 {
+  Util::RATIONAL rat_inf(Util::INT_Inf, 1);
+  if (le.slope() == rat_inf || le.slope() > rat_inf)
+    return;
+
   Util::RATIONAL st_rat(i.step()), min_rat(i.minElem());
   Util::RATIONAL im_step = st_rat * le.slope()
                  , im_begin = min_rat * le.slope() + le.offset();
   if (im_step.denominator() != 1 || im_begin.denominator() != 1) 
     if (im_step.numerator() != 0 || im_begin.numerator() != 0) {
-      Util::ERROR("LIB::SBGMap::compatible: uncompatible i/le");
+      Util::ERROR("LIB::SBGMap::compatible: incompatible i/le");
       return;
     }
 
@@ -63,6 +67,9 @@ Interval image(Interval i, LExp le) {
   }
 
   Util::RATIONAL rat_inf(Util::INT_Inf, 1);
+  // TODO: is it right? if (m == rat_inf && h == -rat_inf)
+  //  return i;
+
   if (m == rat_inf || m > rat_inf)
     return Interval(0, 1, Util::Inf);
 
@@ -186,10 +193,22 @@ Set SBGMap<Set>::preImage() const { return dom_; }
 template<typename Set>
 Set SBGMap<Set>::preImage(const Set &subcodom) const
 {
-  SBGMap inv(image(), exp_.inverse());
   Set im = image();
   Set cap_subcodom = im.intersection(subcodom);
-  Set inv_im = inv.image(cap_subcodom);
+  if (cap_subcodom.isEmpty()) 
+    return Set();
+
+  Util::MD_NAT min = cap_subcodom.minElem();
+  Exp inv_exp = exp_.inverse();
+  for (unsigned int j = 0; j < inv_exp.size(); ++j) {
+    Util::RATIONAL m = exp_[j].slope(), h = exp_[j].offset();
+    Util::RATIONAL rat_inf(Util::INT_Inf, 1);
+    if (m == rat_inf && h == -rat_inf) {
+      inv_exp[j] = LExp(0, min[j]);
+    }
+  }
+  SBGMap inv(cap_subcodom, inv_exp);
+  Set inv_im = inv.image();
 
   return dom_.intersection(inv_im);
 }

sbg/sbg.cpp

@@ -17,7 +17,6 @@
 
  ******************************************************************************/
 
-#include <iostream>
 #include "sbg/sbg.hpp"
 
 namespace SBG {

test/implicit_indices.test

@@ -0,0 +1,25 @@
+// model implicit_indices
+//   Real a[100][100], b[100];
+// equation
+//   for i in 1:100 loop
+//     a[i] = b;
+//   end loop;
+// end implicit_indices
+
+dims = 2
+
+V %= {[0:1:99]x[0:1:99], [101:1:200]x[200:1:299], [100:1:100]x[100:1:199]};
+Vmap %= <<{[100:1:100]x[100:1:199]} -> 0*x+2|0*x+2, {[0:1:199]x[0:1:199]} -> 0*x+1|0*x+1
+          , {[101:1:200]x[200:1:299]} -> 0*x+3|0*x+3>>;
+map1 %= <<{[0:1:99]x[0:1:99]} -> 1*x+101|1*x+200, {[100:1:199]x[100:1:199]} -> 1*x+1|1*x+100>>;
+map2 %= <<{[100:1:199]x[100:1:199]} -> 0*x+100|1*x+0, {[0:1:99]x[0:1:99]} -> 1*x+0|1*x+0>>;
+Emap %= <<{[0:1:99]x[0:1:99]} -> 0*x+1|0*x+1, {[100:1:199]x[100:1:199]} -> 0*x+2|0*x+2>>;
+
+matching(
+  V %= {[0:1:99]x[0:1:99], [101:1:200]x[200:1:299], [100:1:100]x[100:1:199]};
+  Vmap %= <<{[100:1:100]x[100:1:199]} -> 0*x+2|0*x+2, {[0:1:199]x[0:1:199]} -> 0*x+1|0*x+1
+            , {[101:1:200]x[200:1:299]} -> 0*x+3|0*x+3>>;
+  map1 %= <<{[0:1:99]x[0:1:99]} -> 1*x+101|1*x+200, {[100:1:199]x[100:1:199]} -> 1*x+1|1*x+100>>;
+  map2 %= <<{[100:1:199]x[100:1:199]} -> 0*x+100|1*x+0, {[0:1:99]x[0:1:99]} -> 1*x+0|1*x+0>>;
+  Emap %= <<{[0:1:99]x[0:1:99]} -> 0*x+1|0*x+1, {[100:1:199]x[100:1:199]} -> 0*x+2|0*x+2>>;
+)

test/sbg9.test

@@ -1,36 +1,20 @@
-// Matching Test4 of the paper, multiple copies
+// Matching Test1 of the paper, multiple copies
 
 N = 10000000
 
-F1 = N-2
-F2 = N-2+F1
-F3 = N-2+F2
-F4 = 1+F3
-F5 = 1+F4
-F6 = 1+F5
-F7 = 1+F6
-F8 = 1+F7
-F9 = 1+F8
-U1 = N+F9 // a
+F1 = N
+F2 = N+F1
+F3 = N+F2
+U1 = N+F3 // a
 U2 = N+U1 // b
 U3 = N+U2 // der(x)
 
-E1 = N-2 // F1 - a[i]
-E2 = N-2+E1 // F1 - der(x[i])
-E3 = N-2+E2 // F2 - a[i]
-E4 = N-2+E3 // F2 - a[i+1]
-E5 = N-2+E4 // F2 - b[i]
-E6 = N-2+E5 // F3 - b[i]
-E7 = 1+E6 // F4 - b[1]
-E8 = 1+E7 // F5 - b[N]
-E9 = 1+E8 // F6 - a[1]
-E10 = 1+E9 // F6 - der(x[1])
-E11 = 1+E10 // F7 - a[N]
-E12 = 1+E11 // F7 - der(x[N])
-E13 = 1+E12 // F8 - a[1]
-E14 = 1+E13 // F9 - a[N-1]
-E15 = 1+E14 // F9 - a[N]
-E16 = 1+E15 // F9 - b[N-1]
+E1 = N // F1 - a[i]
+E2 = N+E1 // F1 - b[i]
+E3 = N+E2 // F2 - a[i]
+E4 = N+E3 // F2 - b[i]
+E5 = N+E4 // F3 - a[i]
+E6 = N+E5 // F3 - der(x[i])
 
 off1d = 0
 off2d = r(F1, 1)-r(E2, 1)
@@ -66,9 +50,8 @@ off14b = r(U1, 1)-E14
 off15b = r(U1, 1)-E15-1
 off16b = r(U2, 1)-E16-1
 
-V %= {[1:1:F1], [F1+1:1:F2], [F2+1:1:F3], [F3+1:1:F4], [F4+1:1:F5], [F5+1:1:F6]
-  , [F6+1:1:F7], [F7+1:1:F8], [F8+1:1:F9]
-  , [F9+1:1:U1], [U1+1:1:U2], [U2+1:1:U3]};
+V %= {[1:1:F1], [F1+1:1:F2], [F2+1:1:F3]
+  , [F3+1:1:U1], [U1+1:1:U2], [U2+1:1:U3]};
 Vmap %= <<{[1:1:F1]} -> 0*x+1, {[F1+1:1:F2]} -> 0*x+2, {[F2+1:1:F3]} -> 0*x+3
   , {[F3+1:1:F4]} -> 0*x+4, {[F4+1:1:F5]} -> 0*x+5, {[F5+1:1:F6]} -> 0*x+6
   , {[F6+1:1:F7]}-> 0*x+7, {[F7+1:1:F8]} -> 0*x+8, {[F8+1:1:F9]} -> 0*x+9
@@ -118,4 +101,4 @@ Emap %= <<{[1:1:E1]} -> 0*x+1, {[E1+1:1:E2]} -> 0*x+2, {[E2+1:1:E3], [E3+1:1:E4]
   , {[E9+1:1:E10]} -> 0*x+9, {[E10+1:1:E11]} -> 0*x+10, {[E11+1:1:E12]} -> 0*x+11
   , {[E12+1:1:E13]} -> 0*x+12, {[E13+1:1:E14], [E14+1:1:E15]} -> 0*x+13
   , {[E15+1:1:E16]} -> 0*x+14>>;
-  )
+  )