[Boost-commit] svn:boost r75674 - sandbox/variadic_templates/sandbox/stepper/libs/array_stepper/examples

Subject: [Boost-commit] svn:boost r75674 - sandbox/variadic_templates/sandbox/stepper/libs/array_stepper/examples
From: cppljevans_at_[hidden]
Date: 2011-11-26 13:58:15

Author: cppljevans
Date: 2011-11-26 13:58:14 EST (Sat, 26 Nov 2011)
New Revision: 75674
URL: http://svn.boost.org/trac/boost/changeset/75674

Log:
Added loop over all possible axis rotations.
Output=
===>Enter:{{rot_solve
  :axis_rot=0
  lengths_soln={ 4, 2, 2} 
  :axis_soln=0:rank=3:rot_lengths=1:axis_inner=2
  :istk_unrot=
  { offset= 0
  , space=16
  , rotation=0
  , luvs={ { lower=0, upper=0, value=0}, { lower=0, upper=0, value=0}, { lower=0, upper=0, value=0}}
  , rotated indices={ 0, 0, 0} 
  }
  :istk_rot=
  { offset= 0
  , space=16
  , rotation=1
  , luvs={ { lower=0, upper=0, value=0}, { lower=0, upper=0, value=0}, { lower=0, upper=0, value=0}}
  , rotated indices={ 0, 0, 0} 
  }
  r=
  { { { 7, 14} 
    , { 21, 28} 
    } 
  , { { 6, 12} 
    , { 18, 24} 
    } 
  , { { -4, -8} 
    , { -12, -16} 
    } 
  , { { -15, -30} 
    , { -45, -60} 
    } 
  } 
  .
  :istk_rot(axis_inner_at_min)=
  { offset= 0
  , space= 4
  , rotation=1
  , luvs={ { lower=0, upper=3, value=0}, { lower=0, upper=0, value=0}, { lower=0, upper=0, value=0}}
  , rotated indices={ 0, 0, 0} 
  }
  :r_w=
  { { { 7} 
    , { 14} 
    } 
  , { { 21} 
    , { 28} 
    } 
  } 
  .
  :istk_borders_rot=
  { offset= 0
  , space=16
  , rotation=1
  , luvs={ { lower=0, upper=0, value=0}, { lower=0, upper=0, value=0}, { lower=0, upper=0, value=0}}
  , rotated indices={ 0, 0, 0} 
  }
  :rb_w=
  { { { 7, 6, -4, -15} 
    , { 14, 12, -8, -30} 
    } 
  , { { 21, 18, -12, -45} 
    , { 28, 24, -16, -60} 
    } 
  } 
  .
  ***solution***
  x=
  { { { 2, 4} 
    , { 6, 8} 
    } 
  , { { 1, 2} 
    , { 3, 4} 
    } 
  , { { 0, 0} 
    , { 0, 0} 
    } 
  , { { -3, -6} 
    , { -9, -12} 
    } 
  } 
  .
  x_w=
  { { { 2, 1, 0, -3} 
    , { 4, 2, 0, -6} 
    } 
  , { { 6, 3, 0, -9} 
    , { 8, 4, 0, -12} 
    } 
  } 
  .
===>Exit:}}rot_solve
===>Enter:{{rot_solve
  :axis_rot=1
  lengths_soln={ 2, 4, 2} 
  :axis_soln=1:rank=3:rot_lengths=2:axis_inner=2
  :istk_unrot=
  { offset= 0
  , space=16
  , rotation=0
  , luvs={ { lower=0, upper=0, value=0}, { lower=0, upper=0, value=0}, { lower=0, upper=0, value=0}}
  , rotated indices={ 0, 0, 0} 
  }
  :istk_rot=
  { offset= 0
  , space=16
  , rotation=2
  , luvs={ { lower=0, upper=0, value=0}, { lower=0, upper=0, value=0}, { lower=0, upper=0, value=0}}
  , rotated indices={ 0, 0, 0} 
  }
  r=
  { { { 7, 21} 
    , { 6, 18} 
    , { -4, -12} 
    , { -15, -45} 
    } 
  , { { 14, 28} 
    , { 12, 24} 
    , { -8, -16} 
    , { -30, -60} 
    } 
  } 
  .
  :istk_rot(axis_inner_at_min)=
  { offset= 0
  , space= 4
  , rotation=2
  , luvs={ { lower=0, upper=0, value=0}, { lower=0, upper=3, value=0}, { lower=0, upper=0, value=0}}
  , rotated indices={ 0, 0, 0} 
  }
  :r_w=
  { { { 7} 
    , { 14} 
    } 
  , { { 21} 
    , { 28} 
    } 
  } 
  .
  :istk_borders_rot=
  { offset= 0
  , space=16
  , rotation=2
  , luvs={ { lower=0, upper=0, value=0}, { lower=0, upper=0, value=0}, { lower=0, upper=0, value=0}}
  , rotated indices={ 0, 0, 0} 
  }
  :rb_w=
  { { { 7, 6, -4, -15} 
    , { 14, 12, -8, -30} 
    } 
  , { { 21, 18, -12, -45} 
    , { 28, 24, -16, -60} 
    } 
  } 
  .
  ***solution***
  x=
  { { { 2, 6} 
    , { 1, 3} 
    , { 0, 0} 
    , { -3, -9} 
    } 
  , { { 4, 8} 
    , { 2, 4} 
    , { 0, 0} 
    , { -6, -12} 
    } 
  } 
  .
  x_w=
  { { { 2, 1, 0, -3} 
    , { 4, 2, 0, -6} 
    } 
  , { { 6, 3, 0, -9} 
    , { 8, 4, 0, -12} 
    } 
  } 
  .
===>Exit:}}rot_solve
===>Enter:{{rot_solve
  :axis_rot=2
  lengths_soln={ 2, 2, 4} 
  :axis_soln=2:rank=3:rot_lengths=0:axis_inner=2
  :istk_unrot=
  { offset= 0
  , space=16
  , rotation=0
  , luvs={ { lower=0, upper=0, value=0}, { lower=0, upper=0, value=0}, { lower=0, upper=0, value=0}}
  , rotated indices={ 0, 0, 0} 
  }
  :istk_rot=
  { offset= 0
  , space=16
  , rotation=0
  , luvs={ { lower=0, upper=0, value=0}, { lower=0, upper=0, value=0}, { lower=0, upper=0, value=0}}
  , rotated indices={ 0, 0, 0} 
  }
  r=
  { { { 7, 6, -4, -15} 
    , { 14, 12, -8, -30} 
    } 
  , { { 21, 18, -12, -45} 
    , { 28, 24, -16, -60} 
    } 
  } 
  .
  :istk_rot(axis_inner_at_min)=
  { offset= 0
  , space= 4
  , rotation=0
  , luvs={ { lower=0, upper=0, value=0}, { lower=0, upper=0, value=0}, { lower=0, upper=3, value=0}}
  , rotated indices={ 0, 0, 0} 
  }
  :r_w=
  { { { 7} 
    , { 14} 
    } 
  , { { 21} 
    , { 28} 
    } 
  } 
  .
  :istk_borders_rot=
  { offset= 0
  , space=16
  , rotation=0
  , luvs={ { lower=0, upper=0, value=0}, { lower=0, upper=0, value=0}, { lower=0, upper=0, value=0}}
  , rotated indices={ 0, 0, 0} 
  }
  :rb_w=
  { { { 7, 6, -4, -15} 
    , { 14, 12, -8, -30} 
    } 
  , { { 21, 18, -12, -45} 
    , { 28, 24, -16, -60} 
    } 
  } 
  .
  ***solution***
  x=
  { { { 2, 1, 0, -3} 
    , { 4, 2, 0, -6} 
    } 
  , { { 6, 3, 0, -9} 
    , { 8, 4, 0, -12} 
    } 
  } 
  .
  x_w=
  { { { 2, 1, 0, -3} 
    , { 4, 2, 0, -6} 
    } 
  , { { 6, 3, 0, -9} 
    , { 8, 4, 0, -12} 
    } 
  } 
  .
===>Exit:}}rot_solve

Text files modified: 
   sandbox/variadic_templates/sandbox/stepper/libs/array_stepper/examples/solve_tridiag_rotated.cpp |   320 +++++++++++++++++++++++++++------------ 
   1 files changed, 217 insertions(+), 103 deletions(-)

Modified: sandbox/variadic_templates/sandbox/stepper/libs/array_stepper/examples/solve_tridiag_rotated.cpp
==============================================================================
--- sandbox/variadic_templates/sandbox/stepper/libs/array_stepper/examples/solve_tridiag_rotated.cpp	(original)
+++ sandbox/variadic_templates/sandbox/stepper/libs/array_stepper/examples/solve_tridiag_rotated.cpp	2011-11-26 13:58:14 EST (Sat, 26 Nov 2011)
@@ -2,6 +2,8 @@
 #include <boost/array_stepper/numeric/solve_tridiag.hpp>
 #include <boost/array_stepper/array_dyn.hpp>
 #include <boost/array_stepper/array_store_print.hpp>
+#include <boost/assert.hpp>
+
 using namespace boost::array_stepper;
 void test_tridiag()
  /**@brief
@@ -21,7 +23,7 @@
       typedef numeric::tridiag<value_t> trid_t;
       typedef std::vector<value_t> vv_t;
         trid_t 
-      a( numeric::tridiag_seq_seq_tag()
+      a0( numeric::tridiag_seq_seq_tag()
        , std::vector< vv_t>
          ( { { -3.0, -2.0, 1.0} //upper
            , {  5.0,  4.0, 3.0, 5.0} //diagonal
@@ -29,7 +31,7 @@
            }
          )
        );
-      std::cout<<"a.size()="<<a.size()<<"\na=\n"<<a<<"\n";
+      //std::cout<<"a0.size()="<<a0.size()<<"\na=\n"<<a0<<"\n";
     #if 1 
        value_t const
       r0[n]
@@ -37,112 +39,224 @@
         ,   6.0
         , - 4.0
         , -15.0
-        };
-      std::size_t const m=2;
-      std::vector<unsigned> lengths_soln={n,m};
-      typedef ::boost::array_stepper::array_dyn<value_t> array_t;
-      array_t r(lengths_soln);
-      value_t const zero_tol=
-      #if 0
-        0.001
-      #else
-        std::numeric_limits<value_t>::epsilon()
-      #endif
+        }
+        //The rhs (from [MarLop1991])
+        ;
+        std::size_t const 
+      m=2
+        //lengths for axes besides the solution axis.
+        ;
+        std::size_t const 
+      rank=3
+        //number of axes(including solution axis).
         ;
+      auto rot_solve=[ =](std::size_t rot_soln)->void
       {
-              typedef
-            typename array_t::length_strides_t
-          length_strides_t
-            ;
-              typedef
-            typename length_strides_t::value_type::stride_t
-          stride_t
-            ;
-              typedef
-            index_stack_length_stride_crtp_indices<stride_t> 
+          std::vector<unsigned> 
+        lengths_soln(rank,m)
+          //shape of rhs and solution(the x's)
+          ;
+          std::size_t const 
+        axis_soln=(rank+rot_soln)%rank
+          //axis where solutions stored.
+          ;
+        lengths_soln[axis_soln]=n;
+        std::cout<<"lengths_soln="<<lengths_soln<<"\n";
+          std::size_t const 
+        rot_lengths=(axis_soln+1)%rank
+          //This is the rotation which rotates:
+          //  lengths_soln[axis_soln]
+          //to:
+          //  lengths_soln[rank-1]
+          ;
+          typedef ::boost::array_stepper::array_dyn<value_t> 
+        array_t;
+          array_t 
+        r(lengths_soln);
+            typedef
+          typename array_t::length_strides_t
+        length_strides_t;
+            typedef
+          typename length_strides_t::value_type::stride_t
+        stride_t;
+            typedef
+          index_stack_length_stride_crtp_indices<stride_t> 
+        istk_t;
+          length_strides_t const&
+        length_strides_v=r.length_strides();
           istk_t
-            ;
-            length_strides_t const&
-          length_strides_v=r.length_strides()
-            ;
-            istk_t
-          istk_outer
-            ( length_strides_v.begin()+1
-            , length_strides_v.end()
-            )
-            //All but the first axis.
-            ;
-          auto const n_outer=istk_outer.space();
-            istk_t
-          istk_inner
-            ( length_strides_v.begin()
-            , length_strides_v.begin()+1
-            )
-            //Only the first axis.
-            ;
-            value_t
-          ten_power=1.0
-            ;
-          auto const n_inner=istk_inner.space();
-          std::cout<<":n_outer="<<n_outer<<":n_inner="<<n_inner<<"\n";
-          for(int i_outer=0; i_outer<n_outer; ++i_outer, ++istk_outer, ten_power*=10.0)
-          {
-              std::cout<<"ten_power="<<ten_power<<"\n";
-              for(int i_inner=0; i_inner<n_inner; ++i_inner, ++istk_inner)
+        istk_unrot
+          ( length_strides_v.begin()
+          , length_strides_v.end()
+          );
+          istk_t
+        istk_rot(istk_unrot);
+        istk_rot.rotate
+          ( rot_lengths
+          );
+          auto const
+        axis_inner=(axis_soln+rank-rot_lengths)%rank;
+      #define TRACE_SOLVE_ROTATED_TRIDIAG 0
+      #if TRACE_SOLVE_ROTATED_TRIDIAG||1
+        std::cout
+          <<":axis_soln="<<axis_soln
+          <<":rank="<<rank
+          <<":rot_lengths="<<rot_lengths
+          <<":axis_inner="<<axis_inner
+          <<"\n:istk_unrot="<<istk_unrot
+          <<"\n:istk_rot="<<istk_rot
+          <<"\n";
+      #endif
+        BOOST_ASSERT_MSG(axis_inner==rank-1,"axis_inner!=rank-1");
+          auto const
+        ibs=istk_rot[axis_inner]
+          ;
+          auto const
+        stride_inner=ibs.stride_val()
+          //stride of axis_inner axis
+          ;
+          auto const
+        n_inner=ibs.length_val()
+          //length of axis_inner axis
+          ;
+        istk_rot.axis_offsets_put
+          ( axis_inner
+          , 0        //lower border
+          , n_inner-1//upper border
+          )
+          //Fix axis_inner index at it's min.
+          ;
+          auto const 
+        n_outer=istk_rot.space();
+      #if TRACE_SOLVE_ROTATED_TRIDIAG
+        std::cout
+          <<":n_outer="<<n_outer
+          <<":n_inner="<<n_inner
+          <<"\n";
+      #endif
+          auto
+        r_iter=r.data().begin();
+        {
+              value_t
+            factor=1.0;
+            for
+              ( int i_outer=0
+              ; i_outer<n_outer
+              ; ++i_outer
+                , ++istk_rot
+                , factor+=1.0
+              )
               {
-                    auto const 
-                  o = r.offset()
-                    + istk_outer.offset_space_val()
-                    + istk_inner.offset_space_val()
-                    ;
-                  r.data()[o]=ten_power*r0[i_inner];
+                  auto offset_inner=istk_rot.offset_space_val();
+                #if TRACE_SOLVE_ROTATED_TRIDIAG
+                  ::boost::trace_scope ts("init_outer");
+                  std::cout
+                    <<":i_outer="<<i_outer
+                    <<":offset_inner="<<offset_inner
+                    <<":istk_rot(axis_inner_at_min)="<<istk_rot
+                    <<"\n";
+                #endif
+                  for
+                    ( int i_inner=0
+                    ; i_inner<n_inner
+                    ; ++i_inner
+                      , offset_inner+=stride_inner
+                    )
+                  {
+                        value_t
+                      val=factor*r0[i_inner];
+                    #if TRACE_SOLVE_ROTATED_TRIDIAG
+                      ::boost::trace_scope ts("init_inner");
+                      std::cout
+                        <<":i_inner="<<i_inner
+                        <<":val="<<val
+                        <<":offset_inner="<<offset_inner
+                        <<":istk_rot(axis_inner_at_min)="<<istk_rot
+                        <<"\n";
+                    #endif
+                      r_iter[offset_inner]=val;
+                  }
               }
-          }
-          std::cout<<"r=\n"<<r<<".\n";
-      }
-      {
-        #if 1
-              typedef
-            typename array_t::length_strides_t
-          length_strides_t
-            ;
-              typedef
-            typename length_strides_t::value_type::stride_t
-          stride_t
-            ;
-              typedef
-            index_stack_length_stride_crtp_indices<stride_t> 
-          istk_t
-            ;
-            length_strides_t const&
-          length_strides_v=r.length_strides()
-            ;
-            istk_t 
-          istk_v
-            ( length_strides_v.begin()
-            , length_strides_v.end()
-            );
-          array_t x(lengths_soln);
-            auto
-          rhs_iter=r.data().begin();
-          istk_v.rotate(1);
-            array_iter_wrap<value_t> 
-          rhs_w(rhs_iter,istk_v);
-          std::cout<<"rhs_w=\n"
-            <<rhs_w<<".\n";
-          numeric::solve_tridiag<value_t> tds( istk_v, a, rhs_iter, x.data().begin());
-          if( tds.upper_triangulate_trid(zero_tol))
-          {
-              tds.upper_triangulate_rhs();
-              tds.back_substitute();
-              std::cout<<"x=\n"<<x<<".\n";
-              array_iter_wrap<value_t> x_w(x.data().begin(),istk_v);
-              std::cout<<"x_w=\n"
-                <<x_w<<".\n";
-          }
-        #endif
-      }
-    #endif      
+            #if 1
+              std::cout<<"r=\n"<<r<<".\n";
+            #endif
+            #if 1
+                array_iter_wrap<value_t> 
+              r_w(r_iter,istk_rot);
+              std::cout
+                <<":istk_rot(axis_inner_at_min)="<<istk_rot<<"\n"
+                <<":r_w=\n"
+                <<r_w<<".\n";
+            #endif
+            #if 1
+              istk_rot.axes_offsets_put(0,0);
+                array_iter_wrap<value_t> 
+              rb_w(r_iter,istk_rot);
+              std::cout
+                <<":istk_borders_rot="<<istk_rot<<"\n"
+                <<":rb_w=\n"
+                <<rb_w<<".\n";
+            #endif
+        }
+        {
+            trid_t
+          a(a0);
+          #if 1
+          //#define INPLACE
+          //ifdef INPLACE, then solution is done in-place.
+          //i.e. r contains the solution; otherwise,
+          //x (declared below) contains the solution.
+          //
+          #ifndef INPLACE
+              array_t 
+            x(lengths_soln);
+          #endif
+              numeric::solve_tridiag<value_t> 
+            tds
+              ( istk_rot
+              , a
+              , r_iter
+            #ifdef INPLACE
+              , r_iter
+            #else
+              , x.data().begin()
+            #endif
+              );
+              value_t const 
+            zero_tol=
+            #if 0
+              0.001
+            #else
+              std::numeric_limits<value_t>::epsilon()
+            #endif
+              ;
+            if( tds.upper_triangulate_trid(zero_tol))
+            {
+                tds.upper_triangulate_rhs();
+                tds.back_substitute();
+                std::cout<<"***solution***\n";
+              #ifdef INPLACE
+                std::cout<<"r=\n"<<r<<".\n";
+                array_iter_wrap<value_t> r_w(r_iter,istk_rot);
+                std::cout<<"r_w=\n"
+                  <<r_w<<".\n";
+              #else
+                std::cout<<"x=\n"<<x<<".\n";
+                array_iter_wrap<value_t> x_w(x.data().begin(),istk_rot);
+                std::cout<<"x_w=\n"
+                  <<x_w<<".\n";
+              #endif
+            }
+          #endif
+        }
+      #endif
+    };
+    for( std::size_t axis_rot=0; axis_rot<rank; ++axis_rot)
+    {
+        ::boost::trace_scope ts("rot_solve");
+        std::cout<<":axis_rot="<<axis_rot<<"\n";
+        rot_solve(axis_rot);
+    }
   }   
   
 int main()