AbstractLinAlgPack_VectorMutable.cpp

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 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include "AbstractLinAlgPack_VectorMutable.hpp"
00030 #include "AbstractLinAlgPack_VectorMutableSubView.hpp"
00031 #include "AbstractLinAlgPack_VectorSpace.hpp"
00032 #include "RTOp_TOp_assign_scalar.h"
00033 #include "RTOp_TOp_assign_vectors.h"
00034 #include "RTOp_TOp_axpy.h"
00035 #include "RTOp_TOp_set_sub_vector.h"
00036 #include "RTOpPack_RTOpC.hpp"
00037 #include "Teuchos_TestForException.hpp"
00038 
00039 namespace {
00040 
00041 // vector scalar assignment operator
00042 static RTOpPack::RTOpC          assign_scalar_op;
00043 // vector assignment operator
00044 static RTOpPack::RTOpC          assign_vec_op;
00045 // set element operator
00046 static RTOpPack::RTOpC          set_ele_op;
00047 // set sub-vector operator
00048 static RTOpPack::RTOpC          set_sub_vector_op;
00049 // axpy operator
00050 static RTOpPack::RTOpC          axpy_op;
00051 
00052 // Simple class for an object that will initialize the operator objects
00053 class init_rtop_server_t {
00054 public:
00055   init_rtop_server_t() {
00056     // Vector scalar assignment operator
00057     TEST_FOR_EXCEPT(0!=RTOp_TOp_assign_scalar_construct(0.0,&assign_scalar_op.op()));
00058     // Vector assignment operator
00059     TEST_FOR_EXCEPT(0!=RTOp_TOp_assign_vectors_construct(&assign_vec_op.op()));
00060     // Set sub-vector operator
00061     RTOp_SparseSubVector spc_sub_vec;
00062     RTOp_sparse_sub_vector_null(&spc_sub_vec);
00063     TEST_FOR_EXCEPT(0!=RTOp_TOp_set_sub_vector_construct(&spc_sub_vec,&set_sub_vector_op.op()));
00064     // axpy operator
00065     TEST_FOR_EXCEPT(0!=RTOp_TOp_axpy_construct(0.0,&axpy_op.op()));
00066   }
00067 }; 
00068 
00069 // When the program starts, this object will be created and the RTOp_Server object will
00070 // be initialized before main() gets underway!
00071 init_rtop_server_t  init_rtop_server;
00072 
00073 } // end namespace
00074 
00075 namespace AbstractLinAlgPack {
00076 
00077 // VectorMutable
00078 
00079 VectorMutable& VectorMutable::operator=(value_type alpha)
00080 {
00081   TEST_FOR_EXCEPT(0!=RTOp_TOp_assign_scalar_set_alpha(alpha,&assign_scalar_op.op()));
00082   VectorMutable* targ_vecs[1] = { this };
00083   AbstractLinAlgPack::apply_op(assign_scalar_op,0,NULL,1,targ_vecs,NULL);
00084   return *this;
00085 }
00086 
00087 VectorMutable& VectorMutable::operator=(const Vector& vec)
00088 {
00089   if( dynamic_cast<const void*>(&vec) == dynamic_cast<const void*>(this) )
00090     return *this; // Assignment to self!
00091   const Vector*   vecs[1]      = { &vec };
00092   VectorMutable*  targ_vecs[1] = { this };
00093   AbstractLinAlgPack::apply_op(assign_vec_op,1,vecs,1,targ_vecs,NULL);
00094   return *this;
00095 }
00096 
00097 VectorMutable& VectorMutable::operator=(const VectorMutable& vec)
00098 {
00099   return this->operator=(static_cast<const Vector&>(vec));
00100 }
00101 
00102 void VectorMutable::set_ele( index_type i, value_type alpha )
00103 {
00104   TEST_FOR_EXCEPT(0!=RTOp_TOp_assign_scalar_set_alpha(alpha,&assign_scalar_op.op()));
00105   VectorMutable* targ_vecs[1] = { this };
00106   AbstractLinAlgPack::apply_op(
00107     assign_scalar_op,0,NULL,1,targ_vecs,NULL
00108     ,i,1,0 // first_ele, sub_dim, global_offset
00109     );
00110 }
00111 
00112 VectorMutable::vec_mut_ptr_t
00113 VectorMutable::sub_view( const Range1D& rng_in )
00114 {
00115   namespace rcp = MemMngPack;
00116   const index_type dim = this->dim();
00117   const Range1D    rng = rng_in.full_range() ? Range1D(1,dim) : rng_in;
00118 #ifdef TEUCHOS_DEBUG
00119   TEST_FOR_EXCEPTION(
00120     rng.ubound() > dim, std::out_of_range
00121     ,"VectorMutable::sub_view(rng): Error, rng = ["<<rng.lbound()<<","<<rng.ubound()<<"] "
00122     "is not in the range [1,this->dim()] = [1,"<<dim<<"]" );
00123 #endif  
00124   if( rng.lbound() == 1 && rng.ubound() == dim )
00125     return Teuchos::rcp( this, false );
00126   // We are returning a view that could change this vector so we had better
00127   // wipe out the cache
00128   //this->has_changed();  // I don't think this line is needed!
00129   return Teuchos::rcp(
00130     new VectorMutableSubView(
00131       Teuchos::rcp( this, false )
00132       ,rng ) );
00133 }
00134 
00135 void VectorMutable::zero()
00136 {
00137   this->operator=(0.0);
00138 }
00139 
00140 void VectorMutable::axpy( value_type alpha, const Vector& x )
00141 {
00142   TEST_FOR_EXCEPT(0!=RTOp_TOp_axpy_set_alpha(alpha,&axpy_op.op()));
00143   const Vector*  vecs[1]      = { &x   };
00144   VectorMutable* targ_vecs[1] = { this };
00145   AbstractLinAlgPack::apply_op(axpy_op,1,vecs,1,targ_vecs,NULL);
00146 }
00147 
00148 void VectorMutable::get_sub_vector( const Range1D& rng, RTOpPack::MutableSubVector* sub_vec_inout )
00149 {
00150   //
00151   // Here we get a copy of the data for the sub-vector that the
00152   // client will modify.  We must later commit these changes to the
00153   // actual vector when the client calls commitDetachedView(...).
00154   // Note, this implementation is very dependent on the behavior of
00155   // the default implementation of constant version of
00156   // Vector<Scalar>::get_sub_vector(...) and the implementation of
00157   // Vector<Scalar>::set_sub_vector(...)!
00158   //
00159   RTOpPack::SubVector sub_vec;
00160   Vector::get_sub_vector( rng, &sub_vec );
00161   sub_vec_inout->initialize(
00162     sub_vec.globalOffset(),sub_vec.subDim(),const_cast<value_type*>(sub_vec.values()),sub_vec.stride());
00163 }
00164 
00165 void VectorMutable::commit_sub_vector( RTOpPack::MutableSubVector* sub_vec_inout )
00166 {
00167   RTOpPack::SparseSubVector spc_sub_vec(
00168     sub_vec_inout->globalOffset(), sub_vec_inout->subDim()
00169     ,sub_vec_inout->values(), sub_vec_inout->stride()
00170     );
00171   VectorMutable::set_sub_vector( spc_sub_vec );            // Commit the changes!
00172   RTOpPack::SubVector sub_vec(*sub_vec_inout);
00173   Vector::free_sub_vector( &sub_vec );                     // Free the memory!
00174   sub_vec_inout->set_uninitialized();                      // Make null as promised!
00175 }
00176 
00177 void VectorMutable::set_sub_vector( const RTOpPack::SparseSubVector& sub_vec )
00178 {
00179   RTOp_SparseSubVector spc_sub_vec;
00180   if(sub_vec.indices()) {
00181     RTOp_sparse_sub_vector(
00182       sub_vec.globalOffset(), sub_vec.subDim(), sub_vec.subNz()
00183       ,sub_vec.values(), sub_vec.valuesStride(), sub_vec.indices(), sub_vec.indicesStride()
00184       ,sub_vec.localOffset(), sub_vec.isSorted()
00185       ,&spc_sub_vec
00186       );
00187   }
00188   else {
00189     RTOp_SubVector _sub_vec;
00190     RTOp_sub_vector(
00191       sub_vec.globalOffset(), sub_vec.subDim(), sub_vec.values(), sub_vec.valuesStride()
00192       ,&_sub_vec
00193       );
00194     RTOp_sparse_sub_vector_from_dense( &_sub_vec, &spc_sub_vec );
00195   }
00196   RTOpPack::RTOpC  set_sub_vector_op;
00197   TEST_FOR_EXCEPT(0!=RTOp_TOp_set_sub_vector_construct(&spc_sub_vec,&set_sub_vector_op.op()));
00198   VectorMutable* targ_vecs[1] = { this };
00199   AbstractLinAlgPack::apply_op(
00200     set_sub_vector_op,0,NULL,1,targ_vecs,NULL
00201     ,sub_vec.globalOffset()+1,sub_vec.subDim(),sub_vec.globalOffset() // first_ele, sub_dim, global_offset
00202     );
00203 }
00204 
00205 void VectorMutable::Vp_StMtV(
00206   value_type                       alpha
00207   ,const GenPermMatrixSlice        &P
00208   ,BLAS_Cpp::Transp                P_trans
00209   ,const Vector                    &x
00210   ,value_type                      beta
00211   )
00212 {
00213   TEST_FOR_EXCEPTION(
00214     true, std::logic_error
00215     ,"VectorMutable::Vp_StMtV(...): Error, this function has not yet been "
00216     "given a default implementation and has not been overridden for the "
00217     "subclass \'" << typeid(*this).name() << "\'!"
00218     );
00219   // ToDo: Implement using reduction or transformation operators that will
00220   // work with any type of vector.
00221 }
00222 
00223 // Overridden from Vector
00224 
00225 Vector::vec_ptr_t
00226 VectorMutable::sub_view( const Range1D& rng ) const
00227 {
00228   namespace rcp = MemMngPack;
00229   return const_cast<VectorMutable*>(this)->sub_view(rng);
00230 }
00231 
00232 } // end namespace AbstractLinAlgPack

Generated on Thu Sep 18 12:33:52 2008 for AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects by doxygen 1.3.9.1