C++ VectorView类代码示例

    template <> void LapHessenberg(
        MatrixView<double> A, VectorView<double> Ubeta)
    {
        TMVAssert(A.iscm());
        TMVAssert(A.colsize() == A.rowsize());
        TMVAssert(Ubeta.size() == A.rowsize()-1);
        TMVAssert(A.ct()==NonConj);

        int n = A.rowsize();
        int ilo = 1;
        int ihi = n;
        int lda = A.stepj();
        int Lap_info=0;
#ifndef LAPNOWORK
        int lwork = n*LAP_BLOCKSIZE;
        double* work = LAP_DWork(lwork);
#endif
        LAPNAME(dgehrd) (
            LAPCM LAPV(n),LAPV(ilo),LAPV(ihi),
            LAPP(A.ptr()),LAPV(lda),LAPP(Ubeta.ptr())
            LAPWK(work) LAPVWK(lwork) LAPINFO);
#ifdef LAPNOWORK
        LAP_Results(Lap_info,"dgehrd");
#else
        LAP_Results(Lap_info,int(work[0]),m,n,lwork,"dgehrd");
#endif
    }

    static void NonBlockHessenberg(
        MatrixView<T> A, VectorView<T> Ubeta)
    {
#ifdef XDEBUG
        cout<<"Start NonBlock Hessenberg Reduction: A = "<<A<<endl;
        Matrix<T> A0(A);
#endif
        // Decompose A into U H Ut
        // H is a Hessenberg Matrix
        // U is a Unitary Matrix
        // On output, H is stored in the upper-Hessenberg part of A
        // U is stored in compact form in the rest of A along with 
        // the vector Ubeta.
        const ptrdiff_t N = A.rowsize();

        TMVAssert(A.colsize() == A.rowsize());
        TMVAssert(N > 0);
        TMVAssert(Ubeta.size() == N-1);
        TMVAssert(A.iscm() || A.isrm());
        TMVAssert(!Ubeta.isconj());
        TMVAssert(Ubeta.step()==1);

        // We use Householder reflections to reduce A to the Hessenberg form:
        T* Uj = Ubeta.ptr();
        T det = 0; // Ignore Householder det calculations
        for(ptrdiff_t j=0;j<N-1;++j,++Uj) {
#ifdef TMVFLDEBUG
            TMVAssert(Uj >= Ubeta._first);
            TMVAssert(Uj < Ubeta._last);
#endif
            *Uj = Householder_Reflect(A.subMatrix(j+1,N,j,N),det);
            if (*Uj != T(0))
                Householder_LMult(A.col(j+2,N),*Uj,A.subMatrix(0,N,j+1,N).adjoint());
        }

#ifdef XDEBUG
        Matrix<T> U(N,N,T(0));
        U.subMatrix(1,N,1,N) = A.subMatrix(1,N,0,N-1);
        U.upperTri().setZero();
        Vector<T> Ubeta2(N);
        Ubeta2.subVector(1,N) = Ubeta;
        Ubeta2(0) = T(0);
        GetQFromQR(U.view(),Ubeta2);
        Matrix<T> H = A;
        if (N>2) LowerTriMatrixViewOf(H).offDiag(2).setZero();
        Matrix<T> AA = U*H*U.adjoint();
        if (Norm(A0-AA) > 0.001*Norm(A0)) {
            cerr<<"NonBlock Hessenberg: A = "<<Type(A)<<"  "<<A0<<endl;
            cerr<<"A = "<<A<<endl;
            cerr<<"Ubeta = "<<Ubeta<<endl;
            cerr<<"U = "<<U<<endl;
            cerr<<"H = "<<H<<endl;
            cerr<<"UHUt = "<<AA<<endl;
            abort();
        }
#endif
    }

 void StateModel::simulate_initial_state(VectorView eta)const{
   if(eta.size() != state_dimension()){
     std::ostringstream err;
     err << "output vector 'eta' has length " << eta.size()
         << " in StateModel::simulate_initial_state.  Expected length "
         << state_dimension();
     report_error(err.str());
   }
   eta = rmvn(initial_state_mean(), initial_state_variance());
 }

 inline ::Eigen::Map<const ::Eigen::VectorXd, ::Eigen::Unaligned,
                     ::Eigen::InnerStride<::Eigen::Dynamic>>
 EigenMap(const VectorView &view) {
   return ::Eigen::Map<const ::Eigen::VectorXd,
                       ::Eigen::Unaligned,
                       ::Eigen::InnerStride<::Eigen::Dynamic>>(
       view.data(),
       view.size(),
       ::Eigen::InnerStride<::Eigen::Dynamic>(view.stride()));
 }

inline
void invnorm2( const VectorView & x ,
               const ValueView  & r ,
               const ValueView  & r_inv )
{
  Kokkos::parallel_reduce( x.dimension_0() , InvNorm2< VectorView , ValueView >( x , r , r_inv ) );
}

inline
void dot_neg( const VectorView & x ,
              const VectorView & y ,
              const ValueView  & r ,
              const ValueView  & r_neg )
{
  Kokkos::parallel_reduce( x.dimension_0() , DotM< VectorView , ValueView >( x , y , r , r_neg ) );
}

void VectorView::set_to_product(const MatrixView& m, const VectorView& v,
                                const bool transpose)
{
  CBLAS_TRANSPOSE tr;
  if (transpose){
    tr = CblasTrans;
    assert(m.cols() == length());
    assert(m.rows() == v.length());
  } else {
    tr = CblasNoTrans;
    assert(m.cols() == v.length());
    assert(m.rows() == length());
  }

  cblas_dgemv(CblasColMajor, tr, m.rows(), m.cols(), 1.0, m.data(),
              m.stride(), v.data(), 1, 0.0, data_, 1);
}

 double ZGS::increment_log_prior_gradient(const ConstVectorView &parameters,
                                          VectorView gradient) const {
   if (parameters.size() != 1 || gradient.size() != 1) {
     report_error(
         "Wrong size arguments passed to "
         "ZeroMeanGaussianConjSampler::increment_log_prior_gradient.");
   }
   return log_prior(parameters[0], &gradient[0], nullptr);
 }

    static inline void Hessenberg(
        MatrixView<T> A, VectorView<T> Ubeta)
    {
        TMVAssert(A.colsize() == A.rowsize());
        TMVAssert(Ubeta.size() == A.rowsize()-1);
        TMVAssert(A.isrm() || A.iscm());
        TMVAssert(A.ct()==NonConj);
        TMVAssert(Ubeta.step() == 1);

        if (A.rowsize() > 0) {
#ifdef LAP
            if (A.iscm()) 
                LapHessenberg(A,Ubeta);
            else 
#endif
                NonLapHessenberg(A,Ubeta);
        }
    }

double compute_allele_freq(const VectorView snp_genotypes)
{
  double macount = 0;
  size_t ngenos = 0;
  for (size_t i = 0; i < snp_genotypes.length(); ++i){
    if (snp_genotypes(i) >= 0){
      macount += snp_genotypes(i);
      ++ngenos;
    }
  }
  return macount / (2.0 * ngenos);
}

void LMAT::Tmult(VectorView lhs, const ConstVectorView &rhs)const {
    if(lhs.size()!=3) {
        report_error("lhs is the wrong size in LMAT::Tmult");
    }
    if(rhs.size()!=3) {
        report_error("rhs is the wrong size in LMAT::Tmult");
    }
    lhs[0] = rhs[0];
    double phi = phi_->value();
    lhs[1] = rhs[0] + phi * rhs[1];
    lhs[2] = (1-phi) * rhs[1] + rhs[2];
}

// TODO(stevescott):  test
void ASSR::simulate_initial_state(VectorView state0)const {
    // First, simulate the initial state of the client state vector.
    VectorView client_state(state0, 0, state0.size()-2);
    StateSpaceModelBase::simulate_initial_state(client_state);

    // Next simulate the initial value of the first latent weekly
    // observation.
    double mu = StateSpaceModelBase::observation_matrix(0).dot(client_state);
    state0[state_dimension() - 2] = rnorm(mu, regression_->sigma());

    // Finally, the initial state of the cumulator variable is zero.
    state0[state_dimension() - 1] = 0;
}

    void MultMV(
        const T alpha, const GenDiagMatrix<Ta>& A, const GenVector<Tx>& x,
        VectorView<T> y)
        // y (+)= alpha * A * x 
        // yi (+)= alpha * Ai * xi
    {
        TMVAssert(A.size() == x.size());
        TMVAssert(A.size() == y.size());
#ifdef XDEBUG
        //cout<<"MultMV: \n";
        //cout<<"alpha = "<<alpha<<endl;
        //cout<<"A = "<<TMV_Text(A)<<"  "<<A<<endl;
        //cout<<"x = "<<TMV_Text(x)<<"  "<<x<<endl;
        //cout<<"y = "<<TMV_Text(y)<<"  "<<y<<endl;
        Vector<T> y0 = y;
        Vector<Tx> x0 = x;
        Matrix<Ta> A0 = A;
        Vector<T> y2 = alpha*A0*x0;
        if (add) y2 += y0;
#endif

        ElemMultVV<add>(alpha,A.diag(),x,y); 

#ifdef XDEBUG
        if (!(Norm(y-y2) <=
              0.001*(TMV_ABS(alpha)*Norm(A0)*Norm(x0)+
                     (add?Norm(y0):TMV_RealType(T)(0))))) {
            cerr<<"MultMV: alpha = "<<alpha<<endl;
            cerr<<"add = "<<add<<endl;
            cerr<<"A = "<<TMV_Text(A)<<" step "<<A.diag().step()<<"  "<<A0<<endl;
            cerr<<"x = "<<TMV_Text(x)<<" step "<<x.step()<<"  "<<x0<<endl;
            cerr<<"y = "<<TMV_Text(y)<<" step "<<y.step()<<"  "<<y0<<endl;
            cerr<<"-> y = "<<y<<endl;
            cerr<<"y2 = "<<y2<<endl;
            abort();
        }
#endif
    }

    static inline void NonLapHessenberg(
        MatrixView<T> A, VectorView<T> Ubeta)
    {
        TMVAssert(A.rowsize() == A.colsize());
        TMVAssert(A.rowsize() > 0);
        TMVAssert(Ubeta.size() == A.rowsize()-1);

#if 0
        if (A.rowsize() > HESS_BLOCKSIZE)
            BlockHessenberg(A,Ubeta,Vbeta,D,E,det);
        else
#endif
            NonBlockHessenberg(A,Ubeta);
    }

/****************************************************************************
*
*   StreamInfo1::Configure( )
*
*   Set configuration of stream information.
*
****************************************************************************/
bool StreamInfo1::Configure
(
    VectorView& sieves,                     // vector of sieve sizes
    std::vector<PMineralInfo1>& minerals,   // collection of minerals
    double liquidSG                         // density of liquid phase
)
{
    // ensure reasonable number of sieves
    if( sieves.size() <= 2 )
        goto initFailed;

    // ensure reasonable number of minerals
    if( minerals.size() < 1 )
        goto initFailed;

    // set implementation to supplied size distribution
    nSize_ = sieves.size();
    sizes_.resize( nSize_ );
    sizes_ = sieves;

    // set implementation to supplied mineral collection
    nType_ = static_cast<long>( minerals.size() );
    minerals_.clear( );
    minerals_.assign( minerals.begin(), minerals.end() );

    // should probably check for NULL
    //  pointer entries in minerals_

    // succeeded
    return true;

initFailed:

    // Initialization failed - object should not be used
    return false;
}