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

Generated on Thu Sep 18 12:35:13 2008 for MOOCHO (Single Doxygen Collection) by doxygen 1.3.9.1