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;
}

 void DRSM::increment_expected_gradient(
     VectorView gradient, int t, const ConstVectorView &state_error_mean,
     const ConstSubMatrix &state_error_variance) {
   if (gradient.size() != xdim_ || state_error_mean.size() != xdim_ ||
       state_error_variance.nrow() != xdim_ ||
       state_error_variance.ncol() != xdim_) {
     report_error(
         "Wrong size arguments passed to "
         "DynamicRegressionStateModel::increment_expected_gradient.");
   }
   for (int i = 0; i < xdim_; ++i) {
     double mean = state_error_mean[i];
     double var = state_error_variance(i, i);
     double sigsq = DynamicRegressionStateModel::sigsq(i);
     ;
     double tmp = (var + mean * mean) / (sigsq * sigsq) - 1.0 / sigsq;
     gradient[i] += .5 * tmp;
   }
 }

 void swap(VectorView<T>& in_oLHS, VectorView<T>& in_oRHS) { in_oLHS.swap(in_oRHS); }

 void DynamicRegressionStateModel::simulate_state_error(VectorView eta, int t)const{
   check_size(eta.size());
   for (int i = 0; i < eta.size(); ++i) {
     eta[i] = rnorm(0, coefficient_transition_model_[i]->sigma());
   }
 }

    static void BlockHessenberg(
        MatrixView<T> A, VectorView<T> Ubeta)
    {
        // Much like the block version of Bidiagonalize, we try to maintain
        // the operation of several successive Householder matrices in
        // a block form, where the net Block Householder is I - YZYt.
        //
        // But as with the bidiagonlization algorithm (and unlike a simple
        // block QR decomposition), we update the matrix from both the left 
        // and the right, so we also need to keep track of the product
        // ZYtm in addition.
        //
        // The block update at the end of the block loop is
        // m' = (I-YZYt) m (I-YZtYt)
        //
        // The Y matrix is stored in the first K columns of m,
        // and the Hessenberg portion of these columns is updated as we go.
        // For the right-hand-side update, m -= mYZtYt, the m on the right
        // needs to be the full original matrix m, including the original
        // versions of these K columns.  Therefore, we can't wait until 
        // the end for this calculation.  
        //
        // Instead, we keep track of mYZt as we progress, so the final update
        // is:
        //
        // m' = (I-YZYt) (m - mYZt Y)
        //
        // We also need to do this same calculation for each column as we
        // progress through the block.
        //
        const ptrdiff_t N = A.rowsize();

#ifdef XDEBUG
        Matrix<T> A0(A);
#endif

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

        ptrdiff_t ncolmax = MIN(HESS_BLOCKSIZE,N-1);
        Matrix<T,RowMajor> mYZt_full(N,ncolmax);
        UpperTriMatrix<T,NonUnitDiag|ColMajor> Z_full(ncolmax);

        T det(0); // Ignore Householder Determinant calculations
        T* Uj = Ubeta.ptr();
        for(ptrdiff_t j1=0;j1<N-1;) {
            ptrdiff_t j2 = MIN(N-1,j1+HESS_BLOCKSIZE);
            ptrdiff_t ncols = j2-j1;
            MatrixView<T> mYZt = mYZt_full.subMatrix(0,N-j1,0,ncols);
            UpperTriMatrixView<T> Z = Z_full.subTriMatrix(0,ncols);

            for(ptrdiff_t j=j1,jj=0;j<j2;++j,++jj,++Uj) { // jj = j-j1

                // Update current column of A
                //
                // m' = (I - YZYt) (m - mYZt Yt)
                // A(0:N,j)' = A(0:N,j) - mYZt(0:N,0:j) Y(j,0:j)t
                A.col(j,j1+1,N) -= mYZt.Cols(0,j) * A.row(j,0,j).Conjugate();
                //
                // A(0:N,j)'' = A(0:N,j) - Y Z Yt A(0:N,j)'
                // 
                // Let Y = (L)     where L is unit-diagonal, lower-triangular,
                //         (M)     and M is rectangular
                //
                LowerTriMatrixView<T> L = 
                    LowerTriMatrixViewOf(A.subMatrix(j1+1,j+1,j1,j),UnitDiag);
                MatrixView<T> M = A.subMatrix(j+1,N,j1,j);
                // Use the last column of Z as temporary storage for Yt A(0:N,j)'
                VectorView<T> YtAj = Z.col(jj,0,jj);
                YtAj = L.adjoint() * A.col(j,j1+1,j+1);
                YtAj += M.adjoint() * A.col(j,j+1,N);
                YtAj = Z.subTriMatrix(0,jj) * YtAj;
                A.col(j,j1+1,j+1) -= L * YtAj;
                A.col(j,j+1,N) -= M * YtAj;

                // Do the Householder reflection 
                VectorView<T> u = A.col(j,j+1,N);
                T bu = Householder_Reflect(u,det);
#ifdef TMVFLDEBUG
                TMVAssert(Uj >= Ubeta._first);
                TMVAssert(Uj < Ubeta._last);
#endif
                *Uj = bu;

                // Save the top of the u vector, which isn't actually part of u
                T& Atemp = *u.cptr();
                TMVAssert(IMAG(Atemp) == RealType(T)(0));
                RealType(T) Aorig = REAL(Atemp);
                Atemp = RealType(T)(1);

                // Update Z
                VectorView<T> Zj = Z.col(jj,0,jj);
                Zj = -bu * M.adjoint() * u;
                Zj = Z * Zj;
                Z(jj,jj) = -bu;

                // Update mYtZt:
//.........这里部分代码省略.........

 // static methods
 static double dotproduct(const VectorView &v1,
                              const VectorView &v2)
 {
   assert(v1.length() == v2.length());
   return cblas_ddot(v1.length(), v1.data_, 1, v2.data_, 1);
 }