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

Generated on Tue Jul 13 09:30:50 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7