MOOCHO (Single Doxygen Collection) Version of the Day
AbstractLinAlgPack_VectorMutable.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include "AbstractLinAlgPack_VectorMutable.hpp"
00043 #include "AbstractLinAlgPack_VectorMutableSubView.hpp"
00044 #include "AbstractLinAlgPack_VectorSpace.hpp"
00045 #include "RTOp_TOp_assign_scalar.h"
00046 #include "RTOp_TOp_assign_vectors.h"
00047 #include "RTOp_TOp_axpy.h"
00048 #include "RTOp_TOp_set_sub_vector.h"
00049 #include "RTOpPack_RTOpC.hpp"
00050 #include "Teuchos_Assert.hpp"
00051 
00052 namespace {
00053 
00054 // vector scalar assignment operator
00055 static RTOpPack::RTOpC          assign_scalar_op;
00056 // vector assignment operator
00057 static RTOpPack::RTOpC          assign_vec_op;
00058 // set element operator
00059 static RTOpPack::RTOpC          set_ele_op;
00060 // set sub-vector operator
00061 static RTOpPack::RTOpC          set_sub_vector_op;
00062 // axpy operator
00063 static RTOpPack::RTOpC          axpy_op;
00064 
00065 // Simple class for an object that will initialize the operator objects
00066 class init_rtop_server_t {
00067 public:
00068   init_rtop_server_t() {
00069     // Vector scalar assignment operator
00070     TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_assign_scalar_construct(0.0,&assign_scalar_op.op()));
00071     // Vector assignment operator
00072     TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_assign_vectors_construct(&assign_vec_op.op()));
00073     // Set sub-vector operator
00074     RTOp_SparseSubVector spc_sub_vec;
00075     RTOp_sparse_sub_vector_null(&spc_sub_vec);
00076     TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_set_sub_vector_construct(&spc_sub_vec,&set_sub_vector_op.op()));
00077     // axpy operator
00078     TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_axpy_construct(0.0,&axpy_op.op()));
00079   }
00080 }; 
00081 
00082 // When the program starts, this object will be created and the RTOp_Server object will
00083 // be initialized before main() gets underway!
00084 init_rtop_server_t  init_rtop_server;
00085 
00086 } // end namespace
00087 
00088 namespace AbstractLinAlgPack {
00089 
00090 // VectorMutable
00091 
00092 VectorMutable& VectorMutable::operator=(value_type alpha)
00093 {
00094   TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_assign_scalar_set_alpha(alpha,&assign_scalar_op.op()));
00095   VectorMutable* targ_vecs[1] = { this };
00096   AbstractLinAlgPack::apply_op(assign_scalar_op,0,NULL,1,targ_vecs,NULL);
00097   return *this;
00098 }
00099 
00100 VectorMutable& VectorMutable::operator=(const Vector& vec)
00101 {
00102   if( dynamic_cast<const void*>(&vec) == dynamic_cast<const void*>(this) )
00103     return *this; // Assignment to self!
00104   const Vector*   vecs[1]      = { &vec };
00105   VectorMutable*  targ_vecs[1] = { this };
00106   AbstractLinAlgPack::apply_op(assign_vec_op,1,vecs,1,targ_vecs,NULL);
00107   return *this;
00108 }
00109 
00110 VectorMutable& VectorMutable::operator=(const VectorMutable& vec)
00111 {
00112   return this->operator=(static_cast<const Vector&>(vec));
00113 }
00114 
00115 void VectorMutable::set_ele( index_type i, value_type alpha )
00116 {
00117   TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_assign_scalar_set_alpha(alpha,&assign_scalar_op.op()));
00118   VectorMutable* targ_vecs[1] = { this };
00119   AbstractLinAlgPack::apply_op(
00120     assign_scalar_op,0,NULL,1,targ_vecs,NULL
00121     ,i,1,0 // first_ele, sub_dim, global_offset
00122     );
00123 }
00124 
00125 VectorMutable::vec_mut_ptr_t
00126 VectorMutable::sub_view( const Range1D& rng_in )
00127 {
00128   namespace rcp = MemMngPack;
00129   const index_type dim = this->dim();
00130   const Range1D    rng = rng_in.full_range() ? Range1D(1,dim) : rng_in;
00131 #ifdef TEUCHOS_DEBUG
00132   TEUCHOS_TEST_FOR_EXCEPTION(
00133     rng.ubound() > dim, std::out_of_range
00134     ,"VectorMutable::sub_view(rng): Error, rng = ["<<rng.lbound()<<","<<rng.ubound()<<"] "
00135     "is not in the range [1,this->dim()] = [1,"<<dim<<"]" );
00136 #endif  
00137   if( rng.lbound() == 1 && rng.ubound() == dim )
00138     return Teuchos::rcp( this, false );
00139   // We are returning a view that could change this vector so we had better
00140   // wipe out the cache
00141   //this->has_changed();  // I don't think this line is needed!
00142   return Teuchos::rcp(
00143     new VectorMutableSubView(
00144       Teuchos::rcp( this, false )
00145       ,rng ) );
00146 }
00147 
00148 void VectorMutable::zero()
00149 {
00150   this->operator=(0.0);
00151 }
00152 
00153 void VectorMutable::axpy( value_type alpha, const Vector& x )
00154 {
00155   TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_axpy_set_alpha(alpha,&axpy_op.op()));
00156   const Vector*  vecs[1]      = { &x   };
00157   VectorMutable* targ_vecs[1] = { this };
00158   AbstractLinAlgPack::apply_op(axpy_op,1,vecs,1,targ_vecs,NULL);
00159 }
00160 
00161 void VectorMutable::get_sub_vector( const Range1D& rng, RTOpPack::MutableSubVector* sub_vec_inout )
00162 {
00163   //
00164   // Here we get a copy of the data for the sub-vector that the
00165   // client will modify.  We must later commit these changes to the
00166   // actual vector when the client calls commitDetachedView(...).
00167   // Note, this implementation is very dependent on the behavior of
00168   // the default implementation of constant version of
00169   // Vector<Scalar>::get_sub_vector(...) and the implementation of
00170   // Vector<Scalar>::set_sub_vector(...)!
00171   //
00172   RTOpPack::SubVector sub_vec;
00173   Vector::get_sub_vector( rng, &sub_vec );
00174   sub_vec_inout->initialize(
00175     sub_vec.globalOffset(),sub_vec.subDim(),const_cast<value_type*>(sub_vec.values()),sub_vec.stride());
00176 }
00177 
00178 void VectorMutable::commit_sub_vector( RTOpPack::MutableSubVector* sub_vec_inout )
00179 {
00180   RTOpPack::SparseSubVector spc_sub_vec(
00181     sub_vec_inout->globalOffset(), sub_vec_inout->subDim()
00182     ,sub_vec_inout->values(), sub_vec_inout->stride()
00183     );
00184   VectorMutable::set_sub_vector( spc_sub_vec );            // Commit the changes!
00185   RTOpPack::SubVector sub_vec(*sub_vec_inout);
00186   Vector::free_sub_vector( &sub_vec );                     // Free the memory!
00187   sub_vec_inout->set_uninitialized();                      // Make null as promised!
00188 }
00189 
00190 void VectorMutable::set_sub_vector( const RTOpPack::SparseSubVector& sub_vec )
00191 {
00192   RTOp_SparseSubVector spc_sub_vec;
00193   if (!is_null(sub_vec.indices())) {
00194     RTOp_sparse_sub_vector(
00195       sub_vec.globalOffset(), sub_vec.subDim(), sub_vec.subNz(),
00196       sub_vec.values().get(), sub_vec.valuesStride(),
00197       sub_vec.indices().get(), sub_vec.indicesStride(),
00198       sub_vec.localOffset(), sub_vec.isSorted(),
00199       &spc_sub_vec
00200       );
00201   }
00202   else {
00203     RTOp_SubVector _sub_vec;
00204     RTOp_sub_vector(
00205       sub_vec.globalOffset(), sub_vec.subDim(),
00206       sub_vec.values().get(), sub_vec.valuesStride(),
00207       &_sub_vec
00208       );
00209     RTOp_sparse_sub_vector_from_dense( &_sub_vec, &spc_sub_vec );
00210   }
00211   RTOpPack::RTOpC  set_sub_vector_op;
00212   TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_set_sub_vector_construct(&spc_sub_vec,&set_sub_vector_op.op()));
00213   VectorMutable* targ_vecs[1] = { this };
00214   AbstractLinAlgPack::apply_op(
00215     set_sub_vector_op,0,NULL,1,targ_vecs,NULL
00216     ,sub_vec.globalOffset()+1,sub_vec.subDim(),sub_vec.globalOffset() // first_ele, sub_dim, global_offset
00217     );
00218 }
00219 
00220 void VectorMutable::Vp_StMtV(
00221   value_type                       alpha
00222   ,const GenPermMatrixSlice        &P
00223   ,BLAS_Cpp::Transp                P_trans
00224   ,const Vector                    &x
00225   ,value_type                      beta
00226   )
00227 {
00228   TEUCHOS_TEST_FOR_EXCEPTION(
00229     true, std::logic_error
00230     ,"VectorMutable::Vp_StMtV(...): Error, this function has not yet been "
00231     "given a default implementation and has not been overridden for the "
00232     "subclass \'" << typeName(*this) << "\'!"
00233     );
00234   // ToDo: Implement using reduction or transformation operators that will
00235   // work with any type of vector.
00236 }
00237 
00238 // Overridden from Vector
00239 
00240 Vector::vec_ptr_t
00241 VectorMutable::sub_view( const Range1D& rng ) const
00242 {
00243   namespace rcp = MemMngPack;
00244   return const_cast<VectorMutable*>(this)->sub_view(rng);
00245 }
00246 
00247 } // end namespace AbstractLinAlgPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines