自学教程：C++ A0类代码示例

    static type call(A0& a0, I const& indices)    {      Idx idx = nt2::function_index(indices, a0.extent(), meta::as_<typename A0::indexes_type>());      std::size_t b = nt2::run(idx, 0u, meta::as_<std::size_t>());      std::size_t e = nt2::run(idx, nt2::numel(idx.extent())-1u, meta::as_<std::size_t>())+1u;      typename A0::value_type* p = const_cast<typename A0::value_type*>(a0.raw());      return type(expr::make(container_ref(boost::make_shared<container>(idx.extent(), nt2::share(p+b, p+e)))));    }

 BOOST_FORCEINLINE void combine_eigens(const T& wr,  const T& wi, A0& w) {   typedef typename A0::value_type type_t;   int n = numel(wr);   w.resize(of_size(n, 1));   nt2::container::table<type_t, nt2::shared_> sw(of_size(n, 1), share(w.data(), w.data()+n));   sw = tocomplex(wr, wi); }

     BOOST_FORCEINLINE result_type operator()(A0& a0, A1 const& a1) const     {        result_type that;        nt2_la_int n   = nt2::width(a0);        nt2_la_int lda = a0.leading_size();        NT2_F77NAME(zpotrf) (&a1,&n,a0.raw(),&lda,&that);        return that;     }

     BOOST_FORCEINLINE result_type operator()(A0& a0, A1 const a1) const     {        result_type that;        nt2_la_int n   = nt2::width(a0);        nt2_la_int lda = a0.leading_size();        magma_zpotrf(a1,n,(cuDoubleComplex*)a0.raw(),lda,&that);        return that;     }

     BOOST_FORCEINLINE result_type operator()(A0& a0, A1& a1) const     {        result_type that;        nt2_la_int  m  = nt2::height(a0);        nt2_la_int  n  = nt2::width(a0);        nt2_la_int  ld = a0.leading_size();        a1.resize( nt2::of_size(std::min(n, m), 1) );        magma_cgetrf(m, n, (cuFloatComplex*)a0.data(), ld, a1.data(), &that);        return that;     }

 BOOST_FORCEINLINE result_type operator()( A0& a, A1& scale                                         , A2& ilo, A3& ihi, C0 job) const {    result_type info;    nt2_la_int  n   = nt2::width(a);    BOOST_ASSERT_MSG(n == nt2_la_int(nt2::height(a)), "input must be square");    nt2_la_int  lda = nt2::max(a.leading_size(), One<size_t>());    NT2_F77NAME(cgebal) (&job, &n                        , a.raw(), &lda                        , &ilo, &ihi                        , scale.raw(), &info);    return info; }

    BOOST_FORCEINLINE static void conf_bounds(const A0& a0, A1& a1,                                              const M2x2 pcov,                                              const value_type& normz,                                              Out0 & p)    {      typedef typename boost::proto::result_of::child_c<A1&,1>::type         Out1;      typedef typename boost::proto::result_of::child_c<A1&,2>::type         Out2;      //      NT2_DISPLAY(a0.extent());   NT2_DISPLAY(a1.extent());      p.resize(a0.extent());      const In0& x  = boost::proto::child_c<0>(a0);      //      NT2_DISPLAY(x);      const In1& a = boost::proto::child_c<1>(a0);      //      NT2_DISPLAY(a);      const In2& b = boost::proto::child_c<2>(a0);      //       NT2_DISPLAY(b);      BOOST_AUTO_TPL(z, x/b);      //       NT2_DISPLAY(z);      //       NT2_DISPLAY(nt2::gammainc(z, a));      //       NT2_DISPLAY(size(p));      p =  nt2::exp(z); //nt2::gammainc(z, a);      //       NT2_DISPLAY(p);      //      Out0& p = boost::proto::child_c<0>(a1);      //      NT2_DISPLAY(p);      BOOST_AUTO_TPL(itp, (p/nt2::oneminus(p)));      //       NT2_DISPLAY(itp);      BOOST_AUTO_TPL(dp, nt2::rec(p*oneminus(p))); // derivative of logit(p) w.r.t. p      //       NT2_DISPLAY(dp);      BOOST_AUTO_TPL(da, dgammainc(z,a)*dp);       // dlogitp/da = dp/da * dlogitp/dp      //       NT2_DISPLAY(da);      BOOST_AUTO_TPL(db, -nt2::exp(a*nt2::log(z)-z-nt2::gammaln(a)-nt2::log(b))* dp); // dlogitp/db = dp/db * dlogitp/dp      //      NT2_DISPLAY(db);      BOOST_AUTO_TPL(varLogitp, pcov(1,1)*sqr(da) + (Two<value_type>()*pcov(1,2)*da + pcov(2,2)*db)*db);      //       NT2_DISPLAY(varLogitp);      //       NT2_DISPLAY(normz);      BOOST_AUTO_TPL(exp_halfwidth, nt2::exp(normz*nt2::sqrt(varLogitp)));      //       NT2_DISPLAY(exp_halfwidth);      Out1 & plo = boost::proto::child_c<1>(a1);      Out2 & pup = boost::proto::child_c<2>(a1);      plo.resize(a0.extent());      pup.resize(a0.extent());      //      std::cout << "1" << std::endl;      plo = itp*exp_halfwidth;      //      std::cout << "2" << std::endl;      pup = itp/exp_halfwidth;      //      std::cout << "3" << std::endl;      plo /= nt2::oneplus(plo);      //      std::cout << "4" << std::endl;      pup /= nt2::oneplus(pup);      //      std::cout << "5" << std::endl;    }

     BOOST_FORCEINLINE result_type operator()(A0& a0, A1& a1) const     {        result_type that;        nt2_la_int  n  = nt2::width(a0);        nt2_la_int  ld = a0.leading_size();        nt2_la_int  ldb = a1.leading_size();        nt2_la_int  nhrs = nt2::width(a1);        char uplo = 'L';        NT2_F77NAME(zposv) ( &uplo, &n, &nhrs, a0.data(), &ld, a1.data(), &ldb                           , &that                            );        return that;     }

    BOOST_FORCEINLINE result_type operator()( A0& a0, A1& s, A2& u,A3 &vt, A4 const jobu            , A5 const jobvt) const    {        result_type that;        details::workspace<typename A0::value_type> w;        nt2_la_int  m  = nt2::height(a0);        nt2_la_int  n  = nt2::width(a0);        nt2_la_int  ld = a0.leading_size();        nt2_la_int ldu = u.leading_size() > 1 ? u.leading_size() : 1 ;        nt2_la_int ldvt= vt.leading_size() > 1 ? vt.leading_size() : 1 ;        nt2_la_int lwork_query = -1;        char c = std::toupper(jobu);        magma_vec_t magma_jobu  = c == 'A' ? MagmaAllVec : c == 'S' ? MagmaSomeVec : c == 'O' ? MagmaOverwriteVec : c == 'N' ? MagmaNoVec : (magma_vec_t)0;        c = std::toupper(jobvt);        magma_vec_t magma_jobvt = c == 'A' ? MagmaAllVec : c == 'S' ? MagmaSomeVec : c == 'O' ? MagmaOverwriteVec : c == 'N' ? MagmaNoVec : (magma_vec_t)0;        magma_sgesvd(magma_jobu,magma_jobvt,m, n, 0, ld, 0, 0, ldu                     , 0, ldvt, w.main()                     , lwork_query, &that                    );        w.prepare_main();        nt2::gesvd(a0,s,u,vt,jobu,jobvt,w);        return that;    }

 //========================================================================== /// INTERNAL ONLY - X = LINSOLVE(A,B) -- symmetric shape BOOST_FORCEINLINE void eval ( A0 const& a0, A1 const& a1 , A2 const& a2 ,boost::mpl::long_<1> const           , nt2::symmetric_ const&) const {   type_t rcond;   nt2::container::table<nt2_la_int> piv = nt2::zeros(a0.leading_size(), 1                                         , nt2::meta::as_<nt2_la_int>() );   boost::proto::child_c<0>(a2).resize(nt2::of_size(a0.leading_size(),1));   NT2_AS_TERMINAL_IN(desired_semantic1,a,a0);   NT2_AS_TERMINAL_IN(desired_semantic,b,a1);   nt2_la_int iter = nt2::sysvx( boost::proto::value(a),boost::proto::value(piv)                              , boost::proto::value(b)                              , boost::proto::value(boost::proto::child_c<0>(a2))                              , rcond);   boost::ignore_unused(iter); }

    result_type operator()(A0& yi, A1& inputs) const    {        yi.resize(inputs.extent());        const child0 & x   =  boost::proto::child_c<0>(inputs);        if (numel(x) <=  1)            BOOST_ASSERT_MSG(numel(x) >  1, "Interpolation requires at least two sample points in each dimension.");        else        {            BOOST_ASSERT_MSG(issorted(x, 'a'), "for 'nearest' interpolation x values must be sorted in ascending order");            const child1 & y   =  boost::proto::child_c<1>(inputs);            BOOST_ASSERT_MSG(numel(x) == numel(y), "The grid vectors do not define a grid of points that match the given values.");            const child2 & xi  =  boost::proto::child_c<2>(inputs);            bool extrap = false;            value_type extrapval = Nan<value_type>();            choices(inputs, extrap, extrapval, N1());            table<index_type>   index = bsearch (x, xi);            table<value_type>  dx    =  xi-x(index);            table<index_type> indexp1 =  oneplus(index);            yi = y(nt2::if_else(lt(nt2::abs(xi-x(index)), nt2::abs(xi-x(indexp1))), index,  indexp1));            value_type  b =  value_type(x(begin_));            value_type  e =  value_type(x(end_));            if (!extrap) yi = nt2::if_else(nt2::logical_or(boost::simd::is_nge(xi, b),                                               boost::simd::is_nle(xi, e)), extrapval, yi);        }        return yi;    }

 BOOST_FORCEINLINE result_type operator()( A0& a, A1& b, A2& alpha, A3& beta, A4& vr, A5& vl) const {   result_type info = 0;   char jobvr = 'V';   char jobvl = 'V';   details::workspace<typename A0::value_type> wk;   nt2_la_int n = nt2::width(a);   BOOST_ASSERT_MSG(  (n == nt2_la_int(nt2::height(a)))                      &&(n == nt2_la_int(nt2::height(b))) , "inputs must be square");   nt2_la_int lda  = nt2::max(a.leading_size(), One<size_t>());   nt2_la_int ldb  = nt2::max(b.leading_size(), One<size_t>());   nt2_la_int ldvl = n;   nt2_la_int ldvr = n;   NT2_F77NAME(zggev) ( &jobvl, &jobvr                      , &n                      , 0 /*a*/, &lda                      , 0 /*b*/, &ldb                      , 0/*alpha*/                      , 0/*beta*/                      , 0/*vl*/, &ldvl                      , 0/*vr*/, &ldvr                      , wk.main(), details::query()                      , wk.reals()                      , &info);   info = nt2::ggev_wvrvl(a, b, alpha, beta, vr, vl, wk);   return info; }

    BOOST_FORCEINLINE    result_type operator()(A0 const& a0, A1 ind, A2 dim) const    {      std::size_t along = dim-1;      typename A0::extent_type ex1 = a0.extent();      ex1[along] = 1;      return boost::proto::make_expr< nt2::tag::line_                                    , container::domain                                    > ( boost::cref(a0)                                      , along                                      , as_index( a0.extent()                                                , as_subscript(ex1, ind)                                                )                                      , nt2::linesstride(a0, dim)                                      );    }

 result_type operator()(A0& out, const A1& in) const {   size_t n =  boost::proto::child_c<0>(in);   out.resize(nt2::of_size(n, n));   nt2::table<value_type, _2D> l = nt2::abs(nt2::pascal(n, 1, nt2::meta::as_<value_type>()));   nt2::table<value_type, _2D> u = nt2::fliplr(nt2::flipud(l));   nt2::table<value_type, _2D> d = nt2::from_diag(nt2::pow(value_type(-2), nt2::_(value_type(0), value_type(n-1))));   out =  nt2::round2even(nt2::mtimes(l, nt2::mtimes(d, u)));   return out; }

    BOOST_FORCEINLINE result_type operator()(A0& a0, A1 a1, nt2::symmetric_ const&) const    {      result_type norm;      char uplo = 'L';      nt2_la_int n = nt2::width(a0);      nt2_la_int ld = n;     if(a1 =='I'|| a1 =='1'|| a1 =='O')     {        nt2::memory::container<tag::table_, result_type, nt2::_2D> work(nt2::of_size(n,1));        norm = NT2_F77NAME(zlansy)( &a1, &uplo, &n, a0.raw(), &ld, work.raw());      }      else      {        norm = NT2_F77NAME(zlansy)( &a1, &uplo, &n, a0.raw(), &ld, 0);      }      return norm;    }

 BOOST_FORCEINLINE void combine_vects(const T1& rv, const T2& wi, A0& v) {   typedef typename A0::value_type type_t;   int n = height(rv);   v.resize(of_size(n, n));   nt2::container::table<type_t, nt2::shared_> sv(of_size(n, n), share(v.data(), v.data()+numel(v)));   for(int j=1; j <= n; ++j)   {     if(wi(j))     {       sv(nt2::_, j  ) = tocomplex(rv(nt2::_, j), rv(nt2::_, j+1));       sv(nt2::_, j+1) = conj(sv(nt2::_, j));       ++j;     }     else       sv(nt2::_, j) = rv(nt2::_, j);   } }

void foo (int x)   {  // This should be a single variable declaration, but instead it is two seperate declarations.  // this is OK here, but a bug when it is handled this way inside of a for loop initialization.  // See test2012_106.C.// #if 1//   class A0 {public: int foo (int x) {return x;}} x2;// #else     for (class A0 {public: int foo (int x) {return x;}} x2; x2.foo (0); x) {}// #endif   }

 result_type operator()(A0& out, const A1& in) const {   BOOST_AUTO_TPL(x,boost::proto::child_c<0>(in));   _2D siz =  boost::proto::value(boost::proto::child_c<1>(in));   size_t n = siz[0];   size_t m = siz[1];   tab_t rnd = nt2::rand(n*m, 1, nt2::meta::as_<value_t>());   value_t p = nt2::numel(x);   out.resize(siz);   out(nt2::_) = x(nt2::iceil(rnd*p));   return out; }

    result_type operator()(A0& out, const A1& in) const    {      size_t n =  boost::proto::child_c<0>(in);      value_type fn =  value_type(n);      table<value_type> ns =  nt2::cons(of_size(1, 3), fn, fn/12, fn/20);      table<i_type> ee(nt2::of_size(1, 3));      table<value_type> ff(nt2::of_size(1, 3));//      nt2::frexp(ns, ff, ee);      nt2::frexp(ns(1), ff(1), ee(1));      nt2::frexp(ns(2), ff(2), ee(2));      nt2::frexp(ns(3), ff(3), ee(3));      BOOST_AUTO_TPL(kk, nt2::find(nt2::logical_and(eq(ff, nt2::Half<value_type>()), is_gtz(ee))));      BOOST_ASSERT_MSG(!nt2::isempty(kk), "n must be of form 2^m or  12*2^m or 20*2^m");      size_t k =  kk(1);      i_type e = nt2::minusone(ee(k));      out.resize(nt2::of_size(n, n));      table<value_type> h;      if (k == 1) //    n = 1 * 2^e;      {        h = nt2::One<value_type>();      }      else if ( k == 2) // N = 12 * 2^e;      {        h = nt2::vertcat(nt2::ones(1,12,meta::as_<value_type>()),                         nt2::horzcat(nt2::ones(11,1,meta::as_<value_type>()),                                      nt2::toeplitz(cons<value_type>(of_size(1, 11), -1, -1,  1, -1, -1, -1,  1,  1,  1, -1,  1),                                                    cons<value_type>(of_size(1, 11), -1,  1, -1,  1,  1,  1, -1, -1, -1,  1, -1)                                        )                           )          );      }      else if ( k == 3) // N = 20 * 2^e;      {          h = vertcat(nt2::ones(1,20,meta::as_<value_type>()),                      nt2::horzcat(nt2::ones(19,1,meta::as_<value_type>()),                              nt2::hankel(cons<value_type>(of_size(1, 19),-1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, 1),                                          cons<value_type>(of_size(1, 19),1, -1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1)                                )                        )            );      }      //  Kronecker product construction.      for (i_type i = 1; i <= e; ++i)      {        table<value_type> h1 =  nt2::vertcat(nt2::horzcat(h, h), nt2::horzcat(h, -h));        h =  h1;  //ALIASING      }      out =  h;      return out;    }

     BOOST_FORCEINLINE result_type operator()(A0& a0, A1& a1, A2& a2) const     {        result_type that;        nt2_la_int  n  = std::min(nt2::height(a0),nt2::width(a0));        nt2_la_int  lda = n;        nt2_la_int  nhrs = nt2::width(a2);        nt2_la_int  ldb = a2.leading_size();        a1.resize(nt2::of_size(n,1));        magma_zgesv(n,nhrs,(cuDoubleComplex*)a0.data(),lda,a1.data()                   ,(cuDoubleComplex*)a2.data(),ldb,&that);        return that;     }

    BOOST_FORCEINLINE result_type operator()(A0 const& a0, A1 const&a1, A2 const& a2) const    {      result_type rcond;      nt2_la_int n = nt2::height(a0);      nt2_la_int ld = n;      nt2_la_int info;      char uplo = 'L';      nt2::memory::container<tag::table_, v_t, nt2::_2D> work(nt2::of_size(2*n,1));      NT2_F77NAME(zsycon) ( &uplo, &n, a0.raw(), &ld, a1.raw(), &a2, &rcond                          , work.raw(), &info                          );      return rcond;    }

     BOOST_FORCEINLINE result_type operator()(A0& a0, A1& a1, A2& a2) const     {        result_type that;        details::workspace<typename A0::value_type> w;        nt2_la_int  m  = nt2::height(a0);        nt2_la_int  n  = nt2::width(a0);        nt2_la_int  ld = a0.leading_size();        nt2_la_int lwork_query = -1;        magma_sgeqp3 (m, n, 0, ld, 0, 0, w.main()                            , lwork_query, &that                            );          w.prepare_main();          nt2::geqp3(a0,a1,a2,w);        return that;     }

     BOOST_FORCEINLINE result_type operator()(A0& a0, A1& a1, A2& a2) const     {        result_type that;        details::workspace<typename A0::value_type> w;        nt2_la_int  m  = nt2::height(a0);        nt2_la_int  n  = nt2::width(a0);        nt2_la_int  ld = a0.leading_size();        NT2_F77NAME(sgeqp3) (&m, &n, 0, &ld, 0, 0, w.main()                            , details::query(), &that                            );          w.prepare_main();          nt2::geqp3(a0,a1,a2,w);        return that;     }

     BOOST_FORCEINLINE result_type operator()(A0& a0, A1& a1,A2& a2) const     {        result_type that;        details::workspace<typename A2::value_type> w;        nt2_la_int  m  = nt2::height(a0);        nt2_la_int  ld = a0.leading_size();        nt2_la_int  k  = a1.leading_size();        char side  = 'L';        char trans = 'N';        NT2_F77NAME(sormqr) (&side,&trans,&m, &m, &k, 0, &ld, 0, 0, &m, w.main()                            , details::query(), &that                            );        w.prepare_main();        nt2::mqr(a0,a1,a2,w);        return that;     }

 result_type operator()(A0& out, const A1& in) const {   size_t n = boost::proto::child_c<0>(in);   BOOST_ASSERT_MSG(n!= 2, "There is no 2x2 magic matrix");   out.resize(nt2::of_size(n, n));   if(n%2 == 1) //Odd order   {     oddOrderMagicSquare(out,n);   }   else if(n%4 == 0)   {     evenx2(out,n);   }   else   {     evenx1(out,n);   }   return out; }

     BOOST_FORCEINLINE result_type operator()(A0& a0, A1& s, A2& u,A3 &vt, A4 jobu, A5 jobvt) const     {        result_type that;        details::workspace<typename A0::value_type> w;        nt2_la_int  m  = nt2::height(a0);        nt2_la_int  n  = nt2::width(a0);        nt2_la_int  ld = a0.leading_size();        nt2_la_int ldu = u.leading_size() > 1 ? u.leading_size() : 1 ;        nt2_la_int ldvt= vt.leading_size() > 1 ? vt.leading_size() : 1 ;        NT2_F77NAME(sgesvd) (&jobu,&jobvt,&m, &n, 0, &ld, 0, 0, &ldu                            , 0, &ldvt, w.main()                            , details::query(), &that                            );        w.resize_main(5*nt2::max(m, n));        nt2::gesvd(a0,s,u,vt,jobu,jobvt,w);        return that;     }

     BOOST_FORCEINLINE result_type operator()(A0& a0, A1& s) const     {        result_type info;        details::workspace<typename A0::value_type> w;        nt2_la_int  m  = nt2::height(a0);        nt2_la_int  n  = nt2::width(a0);        nt2_la_int  ld = a0.leading_size();        nt2_la_int ldu = 1 ;        nt2_la_int ldvt= 1 ;        char job =  'N';        NT2_F77NAME(sgesvd) (&job, &job, &m, &n, 0, &ld, 0, 0, &ldu                            , 0, &ldvt, w.main()                            , details::query(), &info                            );        w.prepare_main();        w.resize_main(5*nt2::max(m, n));        nt2::gesvd_w(a0,s,w);        return info;     }