MOOCHO (Single Doxygen Collection) Version of the Day
AbstractLinAlgPack_Vector.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 
00043 #include "AbstractLinAlgPack_VectorMutable.hpp"
00044 #include "AbstractLinAlgPack_VectorSubView.hpp"
00045 #include "RTOp_ROp_dot_prod.h"
00046 #include "RTOp_ROp_sum.h"
00047 #include "RTOp_ROp_norms.h"
00048 #include "RTOp_ROp_num_nonzeros.h"
00049 #include "RTOp_ROp_get_sub_vector.h"
00050 #include "RTOpPack_RTOpC.hpp"
00051 #include "RTOpPack_print_sub_vector.hpp"
00052 #include "Teuchos_dyn_cast.hpp"
00053 #include "Teuchos_Assert.hpp"
00054 #include "Teuchos_FancyOStream.hpp"
00055 
00056 #include <limits>
00057 #include <ostream>
00058 
00059 #include <assert.h>
00060 
00061 
00062 // Uncomment to ignore cache of reduction data
00063 //#define ALAP_VECTOR_IGNORE_CACHE_DATA
00064 
00065 
00066 namespace {
00067 
00068 
00069 // get element operator
00070 static RTOpPack::RTOpC                               sum_op;
00071 static Teuchos::RCP<RTOpPack::ReductTarget>  sum_targ;
00072 // number nonzros
00073 static RTOpPack::RTOpC                               num_nonzeros_op;
00074 static Teuchos::RCP<RTOpPack::ReductTarget>  num_nonzeros_targ;
00075 // Norm 1
00076 static RTOpPack::RTOpC                               norm_1_op;
00077 static Teuchos::RCP<RTOpPack::ReductTarget>  norm_1_targ;
00078 // Norm 2
00079 static RTOpPack::RTOpC                               norm_2_op;
00080 static Teuchos::RCP<RTOpPack::ReductTarget>  norm_2_targ;
00081 // Norm inf
00082 static RTOpPack::RTOpC                               norm_inf_op;
00083 static Teuchos::RCP<RTOpPack::ReductTarget>  norm_inf_targ;
00084 // get sub-vector operator
00085 static RTOpPack::RTOpC                               get_sub_vector_op;
00086 
00087 
00088 // Simple class for an object that will initialize the RTOp_Server and operators.
00089 class init_rtop_server_t {
00090 public:
00091   init_rtop_server_t() {
00092     // Operator and target for getting a vector element
00093     TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_ROp_sum_construct(&sum_op.op()));
00094     sum_targ = sum_op.reduct_obj_create();
00095     // Operator and target for norm 1
00096     TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_ROp_num_nonzeros_construct(&num_nonzeros_op.op()));
00097     num_nonzeros_targ = num_nonzeros_op.reduct_obj_create();
00098     // Operator and target for norm 1
00099     TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_ROp_norm_1_construct(&norm_1_op.op()));
00100     norm_1_targ = norm_1_op.reduct_obj_create();
00101     // Operator and target for norm 1
00102     TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_ROp_norm_2_construct(&norm_2_op.op()));
00103     norm_2_targ = norm_2_op.reduct_obj_create();
00104     // Operator and target for norm 1
00105     TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_ROp_norm_inf_construct(&norm_inf_op.op()));
00106     norm_inf_targ = norm_inf_op.reduct_obj_create();
00107     // Get sub-vector operator
00108     TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_ROp_get_sub_vector_construct(1,1,&get_sub_vector_op.op()));
00109   }
00110 }; 
00111 
00112 
00113 // When the program starts, this object will be created
00114 init_rtop_server_t  init_rtop_server;
00115 
00116 
00117 } // end namespace
00118 
00119 
00120 namespace AbstractLinAlgPack {
00121 
00122 
00123 Vector::Vector()
00124   :num_nonzeros_(-1), norm_1_(-1.0), norm_2_(-1.0), norm_inf_(-1.0) // uninitalized
00125 {}
00126 
00127 
00128 index_type Vector::dim() const
00129 {
00130   return this->space().dim();
00131 }
00132 
00133 
00134 index_type Vector::nz() const
00135 {
00136 #ifdef ALAP_VECTOR_IGNORE_CACHE_DATA
00137   {
00138 #else
00139   if( num_nonzeros_ < 0 ) {
00140 #endif
00141     num_nonzeros_op.reduct_obj_reinit(num_nonzeros_targ.ptr());
00142     const Vector *vecs[1] = { this };
00143     AbstractLinAlgPack::apply_op(num_nonzeros_op,1,vecs,0,NULL,&*num_nonzeros_targ);
00144     num_nonzeros_ = RTOp_ROp_num_nonzeros_val(num_nonzeros_op(*num_nonzeros_targ));
00145   }
00146   return num_nonzeros_;
00147 }
00148 
00149 
00150 std::ostream& Vector::output(
00151   std::ostream& out_arg, bool print_dim, bool newline,
00152   index_type global_offset
00153   ) const
00154 {
00155   Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::getFancyOStream(Teuchos::rcp(&out_arg,false));
00156   Teuchos::OSTab tab(out);
00157   RTOpPack::SubVector sub_vec;
00158   this->get_sub_vector( Range1D(), &sub_vec );
00159   RTOpPack::SubVector sub_vec_print( sub_vec.globalOffset() + global_offset, sub_vec.subDim(), sub_vec.values(), sub_vec.stride() );
00160   RTOpPack::output(*out,sub_vec_print,print_dim,newline);
00161   this->free_sub_vector( &sub_vec );
00162   return out_arg;
00163 }
00164 
00165 
00166 VectorMutable::vec_mut_ptr_t Vector::clone() const
00167 {
00168   vec_mut_ptr_t
00169     vec = this->space().create_member();
00170   *vec = *this;
00171   return vec;
00172 }
00173 
00174 
00175 value_type Vector::get_ele(index_type i) const
00176 {
00177   sum_op.reduct_obj_reinit(sum_targ.ptr());
00178   const Vector *vecs[1] = { this };
00179   AbstractLinAlgPack::apply_op(
00180     sum_op,1,vecs,0,NULL,&*sum_targ
00181     ,i,1,0 // first_ele, sub_dim, global_offset
00182     );
00183    return RTOp_ROp_sum_val(sum_op(*sum_targ));
00184 }
00185 
00186 
00187 value_type Vector::norm_1() const
00188 {
00189 #ifdef ALAP_VECTOR_IGNORE_CACHE_DATA
00190   {
00191 #else
00192   if( norm_1_ < 0.0 ) {
00193 #endif
00194     norm_1_op.reduct_obj_reinit(norm_1_targ.ptr());
00195     const Vector *vecs[1] = { this };
00196     AbstractLinAlgPack::apply_op(norm_1_op,1,vecs,0,NULL,&*norm_1_targ);
00197     norm_1_ = RTOp_ROp_norm_1_val(norm_1_op(*norm_1_targ));
00198   }
00199   return norm_1_;
00200 }
00201 
00202 
00203 value_type Vector::norm_2() const
00204 {
00205 #ifdef ALAP_VECTOR_IGNORE_CACHE_DATA
00206   {
00207 #else
00208   if( norm_2_ < 0.0 ) {
00209 #endif
00210     norm_2_op.reduct_obj_reinit(norm_2_targ.ptr());
00211     const Vector *vecs[1] = { this };
00212     AbstractLinAlgPack::apply_op(norm_2_op,1,vecs,0,NULL,&*norm_2_targ);
00213     norm_2_ = RTOp_ROp_norm_2_val(norm_2_op(*norm_2_targ));
00214   }
00215   return norm_2_;
00216 }
00217 
00218 
00219 value_type Vector::norm_inf() const
00220 {
00221 #ifdef ALAP_VECTOR_IGNORE_CACHE_DATA
00222   {
00223 #else
00224   if( norm_inf_ < 0.0 ) {
00225 #endif
00226     norm_inf_op.reduct_obj_reinit(norm_inf_targ.ptr());
00227     const Vector *vecs[1] = { this };
00228     AbstractLinAlgPack::apply_op(norm_inf_op,1,vecs,0,NULL,&*norm_inf_targ);
00229     norm_inf_ = RTOp_ROp_norm_inf_val(norm_inf_op(*norm_inf_targ));
00230   }
00231   return norm_inf_;
00232 }
00233 
00234 
00235 value_type Vector::inner_product(  const Vector& v ) const
00236 {
00237   return this->space().inner_prod()->inner_prod(*this,v);
00238 }
00239 
00240 Vector::vec_ptr_t
00241 Vector::sub_view( const Range1D& rng_in ) const
00242 {
00243   namespace rcp = MemMngPack;
00244   const index_type dim = this->dim();
00245   const Range1D    rng = rng_in.full_range() ? Range1D(1,dim) : rng_in;
00246 #ifdef TEUCHOS_DEBUG
00247   TEUCHOS_TEST_FOR_EXCEPTION(
00248     rng.ubound() > dim, std::out_of_range
00249     ,"Vector::sub_view(rng): Error, rng = ["<<rng.lbound()<<","<<rng.ubound()<<"] "
00250     "is not in the range [1,this->dim()] = [1,"<<dim<<"]" );
00251 #endif  
00252   if( rng.lbound() == 1 && rng.ubound() == dim )
00253     return vec_ptr_t( this, false );
00254   return Teuchos::rcp(
00255     new VectorSubView(
00256       vec_ptr_t( this, false )
00257       ,rng ) );
00258 }
00259 
00260 
00261 void Vector::get_sub_vector( const Range1D& rng_in, RTOpPack::SubVector* sub_vec_inout ) const
00262 {
00263   const Range1D rng = rng_in.full_range() ? Range1D(1,this->space().dim()) : rng_in;
00264 #ifdef TEUCHOS_DEBUG
00265   TEUCHOS_TEST_FOR_EXCEPTION(
00266     this->space().dim() < rng.ubound(), std::out_of_range
00267     ,"Vector::get_sub_vector(rng,...): Error, rng = ["<<rng.lbound()<<","<<rng.ubound()
00268     <<"] is not in range = [1,"<<this->space().dim()<<"]" );
00269 #endif
00270   // Free sub_vec if needed (note this is dependent on the implemenation of this operator class!)
00271   if( sub_vec_inout->values() ) {
00272     std::free( (void*)sub_vec_inout->values()  );
00273   }
00274   // Initialize the operator
00275   RTOpPack::RTOpC get_sub_vector_op;
00276   TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_ROp_get_sub_vector_construct(rng.lbound(),rng.ubound(),&get_sub_vector_op.op()));
00277   // Create the reduction object (another sub_vec)
00278   Teuchos::RCP<RTOpPack::ReductTarget> reduct_obj = get_sub_vector_op.reduct_obj_create(); // This is really of type RTOpPack::ConstSubVectorView<Scalar>!
00279   // Perform the reduction (get the sub-vector requested)
00280   const size_t  num_vecs = 1;
00281   const Vector* sub_vecs[num_vecs] = { this };
00282   AbstractLinAlgPack::apply_op(
00283     get_sub_vector_op,num_vecs,sub_vecs,0,NULL,&*reduct_obj
00284     ,rng.lbound(),rng.size(),rng.lbound()-1 // first_ele, sub_dim, global_offset
00285     );
00286   // Set the sub-vector.  Note reduct_obj will go out of scope so the sub_vec parameter will
00287   // own the memory allocated within get_sub_vector_op.create_reduct_obj_raw(...).  This is okay
00288   // since the client is required to call release_sub_vector(...) so release memory!
00289   RTOp_ReductTarget reduct_obj_raw = get_sub_vector_op(*reduct_obj);
00290   RTOp_SubVector sub_vec = RTOp_ROp_get_sub_vector_val(reduct_obj_raw);
00291   sub_vec_inout->initialize(sub_vec.global_offset,sub_vec.sub_dim,sub_vec.values,sub_vec.values_stride);\
00292   reduct_obj.release();  // Do not allow delete to be called!
00293   std::free(reduct_obj_raw); // Now *sub_vec owns the values[] and indices[] arrays!
00294 }
00295 
00296 
00297 void Vector::free_sub_vector( RTOpPack::SubVector* sub_vec ) const
00298 {
00299   // Free sub_vec if needed (note this is dependent on the implemenation of this operator class!)
00300   if( sub_vec->values() )
00301     std::free( (void*)sub_vec->values() );
00302   sub_vec->set_uninitialized();
00303 }
00304 
00305 
00306 void Vector::has_changed() const
00307 {
00308   num_nonzeros_= -1;  // uninitalized;
00309   norm_1_ = norm_2_ = norm_inf_ = -1.0;
00310 }
00311 
00312 
00313 // protected
00314 
00315 
00316 void Vector::finalize_apply_op(
00317   const size_t num_targ_vecs, VectorMutable** targ_vecs
00318   ) const
00319 {
00320   for( size_t k = 0; k < num_targ_vecs; ++k )
00321     targ_vecs[k]->has_changed();
00322 }
00323 
00324 
00325 } // end namespace AbstractLinAlgPack
00326 
00327 
00328 // nonmember functions
00329 
00330 
00331 void AbstractLinAlgPack::apply_op(
00332   const RTOpPack::RTOp       &op
00333   ,const size_t              num_vecs
00334   ,const Vector*             vecs[]
00335   ,const size_t              num_targ_vecs
00336   ,VectorMutable*            targ_vecs[]
00337   ,RTOpPack::ReductTarget    *reduct_obj
00338   ,const index_type          first_ele
00339   ,const index_type          sub_dim
00340   ,const index_type          global_offset
00341   )
00342 {
00343   if(num_vecs)
00344     vecs[0]->apply_op(op,num_vecs,vecs,num_targ_vecs,targ_vecs,reduct_obj,first_ele,sub_dim,global_offset);
00345   else if (num_targ_vecs)
00346     targ_vecs[0]->apply_op(op,num_vecs,vecs,num_targ_vecs,targ_vecs,reduct_obj,first_ele,sub_dim,global_offset);
00347 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines