C++ DVectorSlice类代码示例

void MatrixSymDiagStd::V_InvMtV(
  DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1
  , const SpVectorSlice& sv_rhs2) const
{
  const DVectorSlice diag = this->diag();
  size_type n = diag.size();

  // y = inv(op(A)) * x
  //
  // A is symmetric and diagonal A = diag(diag) so:
  //
  // y(j) = x(j) / diag(j), for j = 1...n
  //
  // x is sparse so take account of this.
  
  DenseLinAlgPack::Vp_MtV_assert_sizes( vs_lhs->size()
    , n, n, trans_rhs1, sv_rhs2.size() );

  for(   SpVectorSlice::const_iterator x_itr = sv_rhs2.begin()
     ; x_itr != sv_rhs2.end()
     ; ++x_itr )
  {
    (*vs_lhs)(x_itr->indice() + sv_rhs2.offset())
      = x_itr->value() / diag(x_itr->indice() + sv_rhs2.offset());
      // Note: The indice x(i) invocations are ranged check
      // if this is compiled into the code.
  }
}

bool DenseLinAlgPack::assert_print_nan_inf( const DVectorSlice& v
  , const std::string & name, bool throw_excpt, std::ostream* out )
{

  bool has_nan_or_inf = false;
  bool printed_header = false;

  for( DVectorSlice::const_iterator v_itr = v.begin(); v_itr != v.end(); ++v_itr ) {
    if( RTOp_is_nan_inf(*v_itr) ) {
      if(out) {
        if(!printed_header) {
          *out
            << "The vector \"" << name
            << "\" has the following NaN or Inf entries\n";
          printed_header = true;
        }
        *out
          << name << "(" << v_itr - v.begin() + 1 << ") = "
          << *v_itr << std::endl;
      }
      has_nan_or_inf = true;
    }
  }
  if( has_nan_or_inf && throw_excpt ) {
    if(out)
      out->flush();
    std::ostringstream omsg;
    omsg
      << "assert_print_nan_inf(...) : Error, the vector named "
      << name << " has at least one element which is NaN or Inf";
    throw NaNInfException( omsg.str() );
  }

  return !has_nan_or_inf;
}

bool DenseLinAlgPack::comp_less(const DVectorSlice& vs, value_type alpha)
{
  DVectorSlice::const_iterator vs_itr = vs.begin();
  const value_type denom = my_max( ::fabs(alpha), 1.0 );
  for(; vs_itr != vs.end(); ++vs_itr)
    if( *vs_itr > alpha ) return false;
  return true;
}

void DenseLinAlgPack::syr(value_type alpha, const DVectorSlice& vs_rhs, DMatrixSliceSym* sym_lhs)
{
  assert_gms_square(sym_lhs->gms());
  MtV_assert_sizes( sym_lhs->gms().rows(), sym_lhs->gms().cols()
    , BLAS_Cpp::no_trans, vs_rhs.dim() );
  BLAS_Cpp::syr( sym_lhs->uplo(), vs_rhs.dim(), alpha, vs_rhs.raw_ptr()
    , vs_rhs.stride(), sym_lhs->gms().col_ptr(1), sym_lhs->gms().max_rows() );
}

bool DenseLinAlgPack::comp(const DVectorSlice& vs1, const DVectorSlice& vs2) {
  DVectorSlice::const_iterator
    vs1_itr = vs1.begin(),
    vs2_itr = vs2.begin();
  for(; vs1_itr != vs1.end() && vs2_itr != vs2.end(); ++vs1_itr, ++vs2_itr)
    if( !_comp(*vs1_itr,*vs2_itr) ) return false;
  return true;
}

void DenseLinAlgPack::Vp_StMtV(DVectorSlice* vs_lhs, value_type alpha, const DMatrixSlice& gms_rhs1
  , BLAS_Cpp::Transp trans_rhs1, const DVectorSlice& vs_rhs2, value_type beta)
{
  Vp_MtV_assert_sizes(vs_lhs->dim(), gms_rhs1.rows()	, gms_rhs1.cols(), trans_rhs1
    , vs_rhs2.dim());
  BLAS_Cpp::gemv(trans_rhs1,gms_rhs1.rows(),gms_rhs1.cols(),alpha,gms_rhs1.col_ptr(1)
    ,gms_rhs1.max_rows(), vs_rhs2.raw_ptr(),vs_rhs2.stride(),beta,vs_lhs->raw_ptr()
    ,vs_lhs->stride());
}

void MatrixSymDiagStd::Vp_StMtV(
  DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
  , const DVectorSlice& vs_rhs2, value_type beta) const
{
  const DVectorSlice diag = this->diag();
  size_type n = diag.size();

  //
  // y = b*y + a * op(A) * x
  //
  DenseLinAlgPack::Vp_MtV_assert_sizes(
    vs_lhs->size(), n, n, trans_rhs1, vs_rhs2.size() );
  //
  // A is symmetric and diagonal A = diag(diag) so:
  //
  // y(j) += a * diag(j) * x(j), for j = 1...n
  //
  if( vs_rhs2.stride() == 1 && vs_lhs->stride() == 1 ) {
    // Optimized implementation
    const value_type
      *d_itr      = diag.raw_ptr(),
      *x_itr      = vs_rhs2.raw_ptr();
    value_type
      *y_itr      = vs_lhs->raw_ptr(),
      *y_end      = y_itr + vs_lhs->size();

    if( beta == 0.0 ) {
      while( y_itr != y_end )
        *y_itr++ = alpha * (*d_itr++) * (*x_itr++);
    }
    else if( beta == 1.0 ) {
      while( y_itr != y_end )
        *y_itr++ += alpha * (*d_itr++) * (*x_itr++);
    }
    else {
      for( ; y_itr != y_end; ++y_itr )
        *y_itr = beta * (*y_itr) + alpha * (*d_itr++) * (*x_itr++);
    }
  }
  else {
    // Generic implementation
    DVectorSlice::const_iterator
      d_itr = diag.begin(),
      x_itr = vs_rhs2.begin();
    DVectorSlice::iterator
      y_itr = vs_lhs->begin(),
      y_end = vs_lhs->end();
    for( ; y_itr != y_end; ++y_itr, ++d_itr, ++x_itr ) {
#ifdef LINALGPACK_CHECK_RANGE
      TEST_FOR_EXCEPT( !(  d_itr < diag.end()  ) );
      TEST_FOR_EXCEPT( !(  x_itr < vs_rhs2.end()  ) );
      TEST_FOR_EXCEPT( !(  y_itr < vs_lhs->end()  ) );
#endif
      *y_itr = beta * (*y_itr) + alpha * (*d_itr) * (*x_itr);
    }
  }
}

void MatrixSymPosDefLBFGS::Vp_DtV( DVectorSlice* y, const DVectorSlice& x ) const
{
  DenseLinAlgPack::Vp_MtV_assert_sizes(
    y->dim(), m_bar_, m_bar_, BLAS_Cpp::no_trans, x.dim() );

  DVectorSlice::const_iterator
    d_itr	= STY_.diag(0).begin(),
    x_itr	= x.begin();
  DVectorSlice::iterator
    y_itr	= y->begin();

  while( y_itr != y->end() )
    *y_itr++ += (*d_itr++) * (*x_itr++);		
}

void AbstractLinAlgPack::add_elements( SpVector* sv_lhs, value_type alpha, const DVectorSlice& vs_rhs
                   , size_type offset, bool add_zeros )
{
  typedef SpVector::element_type ele_t;
  const bool assume_sorted = !sv_lhs->nz() || ( sv_lhs->nz() && sv_lhs->is_sorted() );
  DVectorSlice::const_iterator
    itr = vs_rhs.begin();
  if(add_zeros) {
    for( size_type i = 1; i <= vs_rhs.dim(); ++i )
      sv_lhs->add_element( ele_t( i + offset, alpha * (*itr++) ) );
  }
  else {
    for( size_type i = 1; i <= vs_rhs.dim(); ++i, ++itr )
      if( *itr != 0.0 )
        sv_lhs->add_element( ele_t( i + offset, alpha * (*itr) ) );
  }
  sv_lhs->assume_sorted(assume_sorted);
}

void DenseLinAlgPack::V_InvMtV(DVectorSlice* vs_lhs, const DMatrixSliceTri& tri_rhs1, BLAS_Cpp::Transp trans_rhs1
  , const DVectorSlice& vs_rhs2)
{
  assert_gms_square(tri_rhs1.gms());
  MtV_assert_sizes(tri_rhs1.gms().rows(), tri_rhs1.gms().cols(), trans_rhs1, vs_rhs2.dim());
  Vp_V_assert_sizes( vs_lhs->dim(), tri_rhs1.gms().rows() );
  (*vs_lhs) = vs_rhs2;
  BLAS_Cpp::trsv(tri_rhs1.uplo(),trans_rhs1,tri_rhs1.diag(),tri_rhs1.gms().rows()
    ,tri_rhs1.gms().col_ptr(1),tri_rhs1.gms().max_rows(), vs_lhs->raw_ptr(),vs_lhs->stride());
}

void DenseLinAlgPack::ger(
  value_type alpha, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2
  , DMatrixSlice* gms_lhs )
{
  Vp_MtV_assert_sizes( vs_rhs2.dim(),  gms_lhs->rows(), gms_lhs->cols()
    , BLAS_Cpp::no_trans, vs_rhs1.dim() );
  BLAS_Cpp::ger(
    gms_lhs->rows(), gms_lhs->cols(), alpha
    ,vs_rhs1.raw_ptr(), vs_rhs1.stride()
    ,vs_rhs2.raw_ptr(), vs_rhs2.stride()
    ,gms_lhs->col_ptr(1), gms_lhs->max_rows() );
}

void MatrixHessianRelaxed::Vp_StMtV(
    DVectorSlice* y, value_type a, BLAS_Cpp::Transp M_trans
  , const DVectorSlice& x, value_type b ) const
{
  using BLAS_Cpp::no_trans;
  using BLAS_Cpp::trans;
  using AbstractLinAlgPack::Vp_StMtV;
  //
  // y = b*y + a * M * x
  // 
  //   = b*y + a * [ H  0    ] * [ x1 ]
  //               [ 0  bigM ]   [ x2 ]
  //               
  // =>              
  //               
  // y1 = b*y1 + a*H*x1
  // 
  // y2 = b*y2 + bigM * x2
  //
  LinAlgOpPack::Vp_MtV_assert_sizes(y->size(),rows(),cols(),M_trans,x.size());

  DVectorSlice
    y1 = (*y)(1,n_);
  value_type
    &y2 = (*y)(n_+1);
  const DVectorSlice
    x1 = x(1,n_);
  const value_type
    x2 = x(n_+1);

  // y1 = b*y1 + a*H*x1
  Vp_StMtV( &y1, a, *H_, no_trans, x1, b );

  // y2 = b*y2 + bigM * x2
  if( b == 0.0 )
    y2 = bigM_ * x2;
  else
    y2 = b*y2 + bigM_ * x2;
  
}

QPSolverStats::ESolutionType
QPSolverRelaxedLOQO::imp_solve_qp(
      std::ostream* out, EOutputLevel olevel, ERunTests test_what
    , const DVectorSlice& g, const MatrixOp& G
    , value_type etaL
    , const SpVectorSlice& dL, const SpVectorSlice& dU
    , const MatrixOp* E, BLAS_Cpp::Transp trans_E, const DVectorSlice* b
      , const SpVectorSlice* eL, const SpVectorSlice* eU
    , const MatrixOp* F, BLAS_Cpp::Transp trans_F, const DVectorSlice* f
    , value_type* obj_d
    , value_type* eta, DVectorSlice* d
    , SpVector* nu
    , SpVector* mu, DVectorSlice* Ed
    , DVectorSlice* lambda, DVectorSlice* Fd
  )
{
  using Teuchos::Workspace;
  Teuchos::WorkspaceStore* wss = wsp::default_workspace_store.get();

  const value_type inf_bnd  = std::numeric_limits<value_type>::max();
//	const value_type real_big = 1e+20;
  const value_type real_big = HUGE_VAL;

  const size_type
    nd   = g.size(),
    m_in = E ? b->size() : 0,
    m_eq = F ? f->size() : 0;

  //
  // Create a LOQO QP definition struct
  //

  LOQO *loqo_lp = openlp();
  TEUCHOS_TEST_FOR_EXCEPT( !(  loqo_lp  ) );

  //
  // Setup loqo_r and loqo_b and count the number of actual
  // constraints.
  //

  // LOQO's b vector storage
  MALLOC( loqo_lp->b, m_in+m_eq, double ); // May not use all of this storage
  DVectorSlice loqo_b( loqo_lp->b, m_in+m_eq );
  // LOQO's r vector storage
  MALLOC( loqo_lp->r, m_in+m_eq, double ); // May not use all of this storage
  DVectorSlice loqo_r( loqo_lp->r, m_in+m_eq );
  // Gives status of b.
  //                  /  j : if eL(j) > -inf_bnd
  // loqo_b_stat(k) = |
  //                  \ -j : if eL(j) <= -inf_bnd && eU(j) < +inf_bnd
  //
  // , for k = 1...num_inequal
  //
  Workspace<int>               loqo_b_stat_ws(wss,m_in); // May not use all of this
  DenseLinAlgPack::VectorSliceTmpl<int>  loqo_b_stat(&loqo_b_stat_ws[0],loqo_b_stat_ws.size());
  std::fill( loqo_b_stat.begin(), loqo_b_stat.end(), 0 ); // Initialize to zero

  // Fill up loqo_b, loqo_r and loqo_b_stat
  size_type num_inequal = 0; // The actual number of bouned general inequalities
  if(E) {
    // Read iterators
    AbstractLinAlgPack::sparse_bounds_itr
      eLU_itr( eL->begin(), eL->end(), eL->offset()
           , eU->begin(), eU->end(), eU->offset(), inf_bnd );
    // written iterators
    DVectorSlice::iterator
      b_itr		= loqo_b.begin(),
      r_itr		= loqo_r.begin();
    DenseLinAlgPack::VectorSliceTmpl<int>::iterator
      b_stat_itr  = loqo_b_stat.begin();
    // loop
    for( int k = 1; !eLU_itr.at_end(); ++k, ++eLU_itr, ++b_itr, ++r_itr, ++b_stat_itr, ++num_inequal )
    {
      const size_type j = eLU_itr.indice();
      if(eLU_itr.lbound() > -inf_bnd) {
        *b_itr = eLU_itr.lbound();
        *r_itr = eLU_itr.ubound() >= inf_bnd ? real_big : eLU_itr.ubound() - eLU_itr.lbound();
        *b_stat_itr = j; // We need to make A(k,:) = [ +op(E)(j,:), -b(j) ]
      }
      else {
        TEUCHOS_TEST_FOR_EXCEPT( !( eLU_itr.ubound() < +inf_bnd ) );
        *b_itr = -eLU_itr.ubound();
        *r_itr = eLU_itr.lbound() <= -inf_bnd ? real_big : - eLU_itr.lbound() + eLU_itr.ubound();
        *b_stat_itr = -j; // We need to make A(k,:) = [ -op(E)(j,:), +b(j) ]
      }
    }
  }
  if(F) {
    LinAlgOpPack::V_StV( &loqo_b(num_inequal+1,num_inequal+m_eq), -1.0, *f );
    loqo_r(num_inequal+1,num_inequal+m_eq) = 0.0;
  }

  //
  // Setup the QP dimensions
  //

  loqo_lp->n = nd+1;
  loqo_lp->m = num_inequal + m_eq;

  //
//.........这里部分代码省略.........

void QPInitFixedFreeStd::initialize(
    const DVectorSlice                   &g
    ,const MatrixSymOp              &G
    ,const MatrixOp                 *A
    ,size_type                          n_R
    ,const size_type                    i_x_free[]
    ,const size_type                    i_x_fixed[]
    ,const EBounds                      bnd_fixed[]
    ,const DVectorSlice                  &b_X
    ,const MatrixSymOpNonsing    &Ko
    ,const DVectorSlice                  &fo
    ,Constraints                        *constraints
    ,std::ostream                       *out
    ,bool                               test_setup
    ,value_type                         warning_tol
    ,value_type                         error_tol
    ,bool                               print_all_warnings
)
{
    namespace GPMSTP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;

    if(!constraints)
        throw std::invalid_argument( "QPInitFixedFreeStd::initialize(...) : Error, "
                                     "constraints == NULL is not allowed." );

    // Validate the sizes of the input arguments
    const size_type
    n = constraints->n(),
    n_X = n - n_R;
    if( n_R > n )
        throw std::invalid_argument( "QPInitFixedFreeStd::initialize(...) : Error, "
                                     "n_R > constraints->n() is not allowed." );
    if(g.dim() !=n)
        throw std::invalid_argument( "QPInitFixedFreeStd::initialize(...) : Error, "
                                     "g.dim() != constraints->n()." );
    if(G.rows() != n || G.cols() !=  n)
        throw std::invalid_argument( "QPInitFixedFreeStd::initialize(...) : Error, "
                                     "G.rows() != constraints->n() or G.cols() !=  constraints->n()." );
    size_type
    m = 0;
    if(A) {
        m = A->cols();
        if( A->rows() != n )
            throw std::invalid_argument( "QPInitFixedFreeStd::initialize(...) : Error, "
                                         "A->rows() != constraints->n()." );
    }
    if(b_X.dim() != n_X)
        throw std::invalid_argument( "QPInitFixedFreeStd::initialize(...) : Error, "
                                     "b_X.dim() != constraints->n() - n_R." );
    if(Ko.rows() != n_R+m || Ko.cols() !=  n_R+m)
        throw std::invalid_argument( "QPInitFixedFreeStd::initialize(...) : Error, "
                                     "Ko.rows() != n_R+A->cols() or Ko.cols() !=  n_R+A->cols()." );
    if(fo.dim() != n_R+m)
        throw std::invalid_argument( "QPInitFixedFreeStd::initialize(...) : Error, "
                                     "fo.dim() != n_R+A->cols()." );

    // Setup x_init, l_x_X_map, i_x_X_map

    const int NOT_SET_YET = -9999; // Can't be FREE, LOWER, UPPER, or EQUALITY
    if(test_setup)
        x_init_.assign( n, (EBounds)NOT_SET_YET );
    else
        x_init_.resize(n);
    l_x_X_map_.assign(n,0);
    i_x_X_map_.assign(n_X,0);

    // Set free portion of x_init
    if( i_x_free == NULL ) {
        for( size_type i = 0; i < n_R; ++i )
            x_init_[i] = FREE;
    }
    else {
        if(test_setup) {
            for( const size_type *i_x_R = i_x_free; i_x_R != i_x_free + n_R; ++i_x_R ) {
                if( *i_x_R < 1 || *i_x_R > n ) {
                    std::ostringstream omsg;
                    omsg
                            << "QPInitFixedFreeStd::initialize(...) : Error, "
                            << "i_x_free[" << i_x_R-i_x_free << "] = "
                            << (*i_x_R) << " is out of bounds";
                    throw std::invalid_argument( omsg.str() );
                }
                if( x_init_(*i_x_R) != NOT_SET_YET ) {
                    std::ostringstream omsg;
                    omsg
                            << "QPInitFixedFreeStd::initialize(...) : Error, "
                            << "Duplicate entries for i_x_free[i] = "
                            << (*i_x_R);
                    throw std::invalid_argument( omsg.str() );
                }
                x_init_(*i_x_R) = FREE;
            }
        }
        else {
            for( const size_type *i_x_R = i_x_free; i_x_R != i_x_free + n_R; ++i_x_R ) {
                x_init_(*i_x_R) = FREE;
            }
        }
    }

//.........这里部分代码省略.........

void QPSchurInitKKTSystemHessianRelaxed::initialize_kkt_system(
    const DVectorSlice&    g
    ,const MatrixOp&  G
    ,value_type           etaL
    ,const SpVectorSlice& dL
    ,const SpVectorSlice& dU
    ,const MatrixOp*  F
    ,BLAS_Cpp::Transp     trans_F
    ,const DVectorSlice*   f
    ,const DVectorSlice&   d
    ,const SpVectorSlice& nu
    ,size_type*           n_R
    ,i_x_free_t*          i_x_free
    ,i_x_fixed_t*         i_x_fixed
    ,bnd_fixed_t*         bnd_fixed
    ,j_f_decomp_t*        j_f_decomp
    ,DVector*              b_X
    ,Ko_ptr_t*            Ko
    ,DVector*              fo
) const
{
    using BLAS_Cpp::trans;

    // Validate type of and convert G
    const MatrixSymHessianRelaxNonSing
    *G_relax_ptr = dynamic_cast<const MatrixSymHessianRelaxNonSing*>(&G);

    if( G_relax_ptr == NULL ) {
        init_kkt_full_.initialize_kkt_system(
            g,G,etaL,dL,dU,F,trans_F,f,d,nu,n_R,i_x_free,i_x_fixed,bnd_fixed
            ,j_f_decomp,b_X,Ko,fo);
        return;
    }

    const MatrixSymHessianRelaxNonSing
    &G_relax = *G_relax_ptr;

    // get some stuff
    const MatrixSymWithOpFactorized
    &G_orig = G_relax.G(),
     &M      = G_relax.M();
    const size_type
    nd = g.size(),
    no = G_orig.rows(),
    nr = M.rows();
    TEST_FOR_EXCEPT( !(  no + nr == nd  ) );

    // Setup output arguments

    // n_R = nd_R
    *n_R = no;
    // i_x_free.size() == 0 and i_x_free is implicitly identity
    i_x_free->resize(no);
{   for(size_type l = 1; l <= no; ++l ) {
            (*i_x_free)[l-1] = l;
        }
    }
    // i_x_fixed[]
    i_x_fixed->resize(nr+1);
    if(nr) {
        // i_x_fixed[l-1] = no + l, l = 1...nr
        for( size_type l = 1; l <= nr; ++l )
            (*i_x_fixed)[l-1] = no+l;
    }
    (*i_x_fixed)[nr] = nd+1; // extra relaxation is always initially active
    // bnd_fixed[]
    bnd_fixed->resize(nr+1);
    if(nr) {
        // bnd_fixed[l-1] = LOWER, l = 1...nr
        std::fill_n( bnd_fixed->begin(), nr, LOWER );
    }
    (*bnd_fixed)[nr] = LOWER; // relaxation is always initially active
    // j_f_decomp[]
    j_f_decomp->resize(0);
    // b_X
    b_X->resize(nr+1);
    if(nr) {
        // b_X[l-1] = dL(no+l), l = 1...nr
        LinAlgOpPack::assign( &(*b_X)(1,nr), dL(no+1,no+nr) );
    }
    (*b_X)[nr] = etaL; // relaxation is always initially active
    // Ko = G.G
    *Ko = G_relax.G_ptr(); // now B_RR is a shared object
    // fo = - *g(1:no)
    LinAlgOpPack::V_StV( fo, -1.0, g(1,no) );

}