Subject: [Boost-commit] svn:boost r75488 - sandbox/variadic_templates/sandbox/stepper/libs/array_stepper/examples
From: cppljevans_at_[hidden]
Date: 2011-11-14 07:54:50
Author: cppljevans
Date: 2011-11-14 07:54:49 EST (Mon, 14 Nov 2011)
New Revision: 75488
URL: http://svn.boost.org/trac/boost/changeset/75488
Log:
solve_tridiag demo with rotated axes
Added:
sandbox/variadic_templates/sandbox/stepper/libs/array_stepper/examples/solve_tridiag_rotated.cpp (contents, props changed)
Added: sandbox/variadic_templates/sandbox/stepper/libs/array_stepper/examples/solve_tridiag_rotated.cpp
==============================================================================
--- (empty file)
+++ sandbox/variadic_templates/sandbox/stepper/libs/array_stepper/examples/solve_tridiag_rotated.cpp 2011-11-14 07:54:49 EST (Mon, 14 Nov 2011)
@@ -0,0 +1,153 @@
+#include <boost/utility/trace_scope.hpp>
+#include <boost/array_stepper/numeric/solve_tridiag.hpp>
+#include <boost/array_stepper/array_dyn.hpp>
+#include <boost/array_stepper/array_store_print.hpp>
+using namespace boost::array_stepper;
+void test_tridiag()
+ /**@brief
+ * Test solve_tridiag with rotated rhs.
+ * The test problem is that
+ * from Example 3.2E on page 136 of [MarLop1991].
+ **@references
+ * [MarLop1991]
+ * Maron, Melvin J.
+ * Lopez, Robert J.
+ * _Numerical Analysis: A Practical Approach_
+ * Wadsworth(1991)
+ */
+ {
+ std::size_t const n=4;
+ typedef double value_t;
+ typedef numeric::tridiag<value_t> trid_t;
+ typedef std::vector<value_t> vv_t;
+ trid_t
+ a( 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
+ , { 1.0, -1.0, 2.0} //lower
+ }
+ )
+ );
+ std::cout<<"a.size()="<<a.size()<<"\na=\n"<<a<<"\n";
+ #if 1
+ value_t const
+ r0[n]
+ = { 7.0
+ , 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
+ ;
+ {
+ 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_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)
+ {
+ auto const
+ o = r.offset()
+ + istk_outer.offset_space_val()
+ + istk_inner.offset_space_val()
+ ;
+ r.data()[o]=ten_power*r0[i_inner];
+ }
+ }
+ 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
+ }
+
+int main()
+{
+ boost::iostreams::indent_scoped_ostreambuf<char> indent_outbuf(std::cout,2);
+ test_tridiag();
+ return 0;
+}