ROL_Vector.hpp

Go to the documentation of this file.

// @HEADER
// ************************************************************************
//
//               Rapid Optimization Library (ROL) Package
//                 Copyright (2014) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact lead developers:
//              Drew Kouri   (dpkouri@sandia.gov) and
//              Denis Ridzal (dridzal@sandia.gov)
//
// ************************************************************************
// @HEADER

#ifndef ROL_VECTOR_H
#define ROL_VECTOR_H

#define ROL_UNUSED(x) (void) x

#include <ostream>
#include <vector>
#include <algorithm>

#include "ROL_Elementwise_Function.hpp"

#include "ROL_Ptr.hpp"
#include "ROL_Stream.hpp"

namespace ROL {

template <class Real>
class Vector
#ifdef ENABLE_PYROL
 : public std::enable_shared_from_this<Vector<Real>> 
#endif 
{
public:

  virtual ~Vector() {}


  virtual void plus( const Vector &x ) = 0;


  virtual void scale( const Real alpha ) = 0;


  virtual Real dot( const Vector &x ) const = 0;


  virtual Real norm() const = 0;


  virtual ROL::Ptr<Vector> clone() const = 0;


  virtual void axpy( const Real alpha, const Vector &x ) {
    ROL::Ptr<Vector> ax = x.clone();
    ax->set(x);
    ax->scale(alpha);
    this->plus(*ax);
  }

  virtual void zero() {
    this->scale( (Real)0 );
  }


  virtual ROL::Ptr<Vector> basis( const int i ) const {
    ROL_UNUSED(i);
    return ROL::nullPtr;
  }


  virtual int dimension() const {return 0;}


  virtual void set( const Vector &x ) {
    this->zero();
    this->plus(x);
  }


  virtual const Vector & dual() const {
    return *this;
  }

  virtual Real apply(const Vector<Real> &x) const {
    return this->dot(x.dual());
  }

  virtual void applyUnary( const Elementwise::UnaryFunction<Real> &f ) {
    ROL_UNUSED(f);
    ROL_TEST_FOR_EXCEPTION( true, std::logic_error,
      "The method applyUnary was called, but not implemented" << std::endl);
  }

  virtual void applyBinary( const Elementwise::BinaryFunction<Real> &f, const Vector &x ) {
    ROL_UNUSED(f);
    ROL_UNUSED(x);
    ROL_TEST_FOR_EXCEPTION( true, std::logic_error,
      "The method applyBinary was called, but not implemented" << std::endl);
  }

  virtual Real reduce( const Elementwise::ReductionOp<Real> &r ) const {
    ROL_UNUSED(r);
    ROL_TEST_FOR_EXCEPTION( true, std::logic_error,
      "The method reduce was called, but not implemented" << std::endl); 
  }

  virtual void print( std::ostream &outStream ) const {
    outStream << "The method print was called, but not implemented" << std::endl;
  }

  virtual void setScalar( const Real C ) {
    this->applyUnary(Elementwise::Fill<Real>(C));
  }

  virtual void randomize( const Real l = 0.0, const Real u = 1.0 ) {
    Elementwise::UniformlyRandom<Real> ur(l,u);
    this->applyUnary(ur);
  }

  virtual std::vector<Real> checkVector( const Vector<Real> &x,
                                         const Vector<Real> &y,
                                         const bool printToStream = true,
                                         std::ostream & outStream = std::cout ) const {
    Real zero =  0.0;
    Real one  =  1.0;
    Real a    =  1.234;
    Real b    = -0.4321;
    int width =  94;
    std::vector<Real> vCheck;

    ROL::nullstream bhs; // outputs nothing

    ROL::Ptr<std::ostream> pStream;
    if (printToStream) {
      pStream = ROL::makePtrFromRef(outStream);
    } else {
      pStream = ROL::makePtrFromRef(bhs);
    }

    // Save the format state of the original pStream.
    ROL::nullstream oldFormatState, headerFormatState;
    oldFormatState.copyfmt(*pStream);

    ROL::Ptr<Vector> v    = this->clone();
    ROL::Ptr<Vector> vtmp = this->clone();
    ROL::Ptr<Vector> xtmp = x.clone();
    ROL::Ptr<Vector> ytmp = y.clone();

    *pStream << "\n" << std::setw(width) << std::left << std::setfill('*') << "********** Begin verification of linear algebra. " << "\n\n";
    headerFormatState.copyfmt(*pStream);

    // Commutativity of addition.
    v->set(*this); xtmp->set(x); ytmp->set(y);
    v->plus(x); xtmp->plus(*this); v->axpy(-one, *xtmp); vCheck.push_back(v->norm());
    *pStream << std::scientific << std::setprecision(12) << std::setfill('>');
    *pStream << std::setw(width) << std::left << "Commutativity of addition. Consistency error: " << " " << vCheck.back() << "\n";

    // Associativity of addition.
    v->set(*this); xtmp->set(x); ytmp->set(y);
    ytmp->plus(x); v->plus(*ytmp); xtmp->plus(*this); xtmp->plus(y); v->axpy(-one, *xtmp); vCheck.push_back(v->norm());
    *pStream << std::setw(width) << std::left << "Associativity of addition. Consistency error: " << " " << vCheck.back() << "\n";

    // Identity element of addition.
    v->set(*this); xtmp->set(x); ytmp->set(y);
    v->zero(); v->plus(x); v->axpy(-one, x); vCheck.push_back(v->norm());
    *pStream << std::setw(width) << std::left << "Identity element of addition. Consistency error: " << " " << vCheck.back() << "\n";

    // Inverse elements of addition.
    v->set(*this); xtmp->set(x); ytmp->set(y);
    v->scale(-one); v->plus(*this); vCheck.push_back(v->norm());
    *pStream << std::setw(width) << std::left << "Inverse elements of addition. Consistency error: " << " " << vCheck.back() << "\n";

    // Identity element of scalar multiplication.
    v->set(*this); xtmp->set(x); ytmp->set(y);
    v->scale(one); v->axpy(-one, *this); vCheck.push_back(v->norm());
    *pStream << std::setw(width) << std::left << "Identity element of scalar multiplication. Consistency error: " << " " << vCheck.back() << "\n";

    // Consistency of scalar multiplication with field multiplication.
    v->set(*this); vtmp->set(*this);
    v->scale(b); v->scale(a); vtmp->scale(a*b); v->axpy(-one, *vtmp); vCheck.push_back(v->norm());
    *pStream << std::setw(width) << std::left << "Consistency of scalar multiplication with field multiplication. Consistency error: " << " " << vCheck.back() << "\n";

    // Distributivity of scalar multiplication with respect to field addition.
    v->set(*this); vtmp->set(*this);
    v->scale(a+b); vtmp->scale(a); vtmp->axpy(b, *this); v->axpy(-one, *vtmp); vCheck.push_back(v->norm());
    *pStream << std::setw(width) << std::left << "Distributivity of scalar multiplication with respect to field addition. Consistency error: " << " " << vCheck.back() << "\n";

    // Distributivity of scalar multiplication with respect to vector addition.
    v->set(*this); xtmp->set(x); ytmp->set(y);
    v->plus(x); v->scale(a); xtmp->scale(a); xtmp->axpy(a, *this); v->axpy(-one, *xtmp); vCheck.push_back(v->norm());
    *pStream << std::setw(width) << std::left << "Distributivity of scalar multiplication with respect to vector addition. Consistency error: " << " " << vCheck.back() << "\n";

    // Commutativity of dot (inner) product over the field of reals.
    vCheck.push_back(std::abs(this->dot(x) - x.dot(*this)));
    *pStream << std::setw(width) << std::left << "Commutativity of dot (inner) product over the field of reals. Consistency error: " << " " << vCheck.back() << "\n";

    // Additivity of dot (inner) product.
    xtmp->set(x);
    xtmp->plus(y); vCheck.push_back(std::abs(this->dot(*xtmp) - this->dot(x) - this->dot(y))/std::max({static_cast<Real>(std::abs(this->dot(*xtmp))), static_cast<Real>(std::abs(this->dot(x))), static_cast<Real>(std::abs(this->dot(y))), one}));
    *pStream << std::setw(width) << std::left << "Additivity of dot (inner) product. Consistency error: " << " " << vCheck.back() << "\n";

    // Consistency of scalar multiplication and norm.
    v->set(*this);
    Real vnorm = v->norm();
    if (vnorm == zero) {
      v->scale(a);
      vCheck.push_back(std::abs(v->norm() - zero));
    } else {
      v->scale(one/vnorm);
      vCheck.push_back(std::abs(v->norm() - one));
    }
    *pStream << std::setw(width) << std::left << "Consistency of scalar multiplication and norm. Consistency error: " << " " << vCheck.back() << "\n";

    // Reflexivity.
    v->set(*this);
    xtmp = ROL::constPtrCast<Vector>(ROL::makePtrFromRef(this->dual()));
    ytmp = ROL::constPtrCast<Vector>(ROL::makePtrFromRef(xtmp->dual()));
    v->axpy(-one, *ytmp); vCheck.push_back(v->norm());
    *pStream << std::setw(width) << std::left << "Reflexivity. Consistency error: " << " " << vCheck.back() << "\n";

    // Consistency of apply and dual.
    v->set(*this);
    xtmp = x.dual().clone(); xtmp->set(x.dual());
    Real vx  = v->apply(*xtmp);
    Real vxd = v->dot(xtmp->dual());
    Real vdx = xtmp->dot(v->dual());
    if (vx == zero) {
      vCheck.push_back(std::max(std::abs(vx-vxd),std::abs(vx-vdx)));
    }
    else {
      vCheck.push_back(std::max(std::abs(vx-vxd),std::abs(vx-vdx))/std::abs(vx));
    }
    *pStream << std::setw(width) << std::left << "Consistency of apply and dual:" << " " << vCheck.back() << "\n\n";

    //*pStream << "************   End verification of linear algebra.\n\n";

    // Restore format state of pStream used for the header info.
    pStream->copyfmt(headerFormatState);
    *pStream << std::setw(width) << std::left << "********** End verification of linear algebra. " << "\n\n";

    // Restore format state of the original pStream.
    pStream->copyfmt(oldFormatState);

    return vCheck;
  }

}; // class Vector

} // namespace ROL

#endif