RTOpPack_SUNDIALS_Ops.hpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // RTOp: Interfaces and Support Software for Vector Reduction Transformation
00005 //       Operations
00006 //                Copyright (2006) Sandia Corporation
00007 // 
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //  
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //  
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00026 // 
00027 // ***********************************************************************
00028 // @HEADER
00029 
00030 #ifndef RTOPPACK_SUNDIALS_OPS_HPP
00031 #define RTOPPACK_SUNDIALS_OPS_HPP
00032 
00033 #include "RTOpPack_RTOpTHelpers.hpp"
00034 #include "RTOpPack_RTOpT.hpp"
00035 #include "RTOpPack_Types.hpp"
00036 
00037 
00038 namespace RTOpPack 
00039 {
00040   /* ------------------------------------------------------------------------------
00041    *
00042    * RTOp implementation of operations required by the SUNDIALS N_Vector interface
00043    * but not present in the Thyra StdVectorOps set. A few of these methods are
00044    * minor variations on the standard ops and could be composed from them, but 
00045    * composing them in that way requires vector copies not required by the operations
00046    * themselves. Thus for the sake of performance, we've implemented specialized RTOps
00047    * for all N_Vector functions not supported exactly as standard ops. 
00048    *
00049    * The following methods can be implemented directly with standard operations:
00050    * 
00051    * N_VLinearSum   -->  Thyra::linear_combination()
00052    * N_VConst       -->  Thyra::put_scalar()
00053    * N_VAbs         -->  Thyra::abs()
00054    * N_VInv         -->  Thyra::reciprocal() 
00055    * N_VDotProd     -->  Thyra::dot() 
00056    * N_VL1Norm      -->  Thyra::norm_1() 
00057    * N_VMaxNorm     -->  Thyra::norm_inf()
00058    * N_VMin         -->  Thyra::min() 
00059    *
00060    * The following, while similar to functions in the standard ops set, must 
00061    * be implemented here for optimal efficiency. 
00062    * 
00063    * N_VProd            - elementwise product
00064    * N_VDiv             - elementwise divide
00065    * N_VScale           - scale by constant
00066    * N_VAddConst        - add a constant to all elements
00067    *
00068    *
00069    * Weighted norms:
00070    *
00071    * The Thyra standard ops do include a weighted norm, but N_Vector uses a different
00072    * definition of the weighting factor. Thyra defines the weighted L2Norm as 
00073    * Sqrt[Sum[w_i x_i^2]],  but SUNDIALS wants Sqrt[Sum[w_i^2 x_i^2]]. 
00074    * 
00075    * N_VWL2Norm         - weighted 2-norm
00076    * N_VWrmsNorm        - weighted norm, normalized by vector length
00077    * N_VWrmsNormMask    - weighted rms norm, ignoring certain masked elements
00078    *
00079    * Notice that because RTOps work with chunks of data, 
00080    * the normalization by vector length *cannot* be performed inside
00081    * the RTOp and must be deferred until the Thyra wrapper functions. 
00082    *
00083    * 
00084    * Finally, there are several methods unlike those in the standard ops, which
00085    * are of course implemented anew:
00086    *
00087    * N_VConstrMask
00088    * N_VInvTest
00089    * N_VCompare
00090    * N_VMinQuotient
00091    * 
00092    * We do not implement the data access methods
00093    *
00094    * N_VGetArrayPointer
00095    * N_VSetArrayPointer
00096    *
00097    * which means that the SUNDIALS dense and banded solvers cannot be used. 
00098    *
00099    * \author Kevin Long
00100    * \date 10 December 2005
00101    *
00102    * ------------------------------------------------------------------------------*/
00103 
00104 
00105 
00109   template<class Scalar>
00110   class SUNDIALS_VProd : public ROpScalarTransformationBase<Scalar> 
00111   {
00112   public:
00114     SUNDIALS_VProd()
00115       : RTOpT<Scalar>("SUNDIALS_VProd"), ROpScalarTransformationBase<Scalar>() 
00116     {}
00120     void apply_op(const int num_vecs, 
00121                   const ConstSubVectorView<Scalar> sub_vecs[],
00122                   const int  num_targ_vecs,  
00123                   const SubVectorView<Scalar> targ_sub_vecs[],
00124                   ReductTarget *reduct_obj) const
00125     {
00126       RTOP_APPLY_OP_2_1(num_vecs,sub_vecs,num_targ_vecs,targ_sub_vecs);
00127       for( index_type i = 0; i < subDim; 
00128            ++i,  v0_val += v0_s,  v1_val += v1_s, z0_val += z0_s ) 
00129         {
00130           *z0_val = (*v0_val) * (*v1_val);
00131         }
00132     }
00134   }; 
00135 
00139   template<class Scalar>
00140   class SUNDIALS_VDiv : public ROpScalarTransformationBase<Scalar> 
00141   {
00142   public:
00144     SUNDIALS_VDiv()
00145       : RTOpT<Scalar>("SUNDIALS_VDiv"), ROpScalarTransformationBase<Scalar>() 
00146     {}
00150     void apply_op(const int num_vecs, 
00151                   const ConstSubVectorView<Scalar> sub_vecs[],
00152                   const int  num_targ_vecs,  
00153                   const SubVectorView<Scalar> targ_sub_vecs[],
00154                   ReductTarget *reduct_obj) const
00155     {
00156       RTOP_APPLY_OP_2_1(num_vecs,sub_vecs,num_targ_vecs,targ_sub_vecs);
00157       for( index_type i = 0; i < subDim; 
00158            ++i,  v0_val += v0_s,  v1_val += v1_s, z0_val += z0_s ) 
00159         {
00160           *z0_val = (*v0_val) / (*v1_val);
00161         }
00162     }
00164   }; 
00165 
00166 
00170   template<class Scalar>
00171   class SUNDIALS_VScale : public ROpScalarTransformationBase<Scalar> 
00172   {
00173   public:
00175     void alpha( const Scalar& alpha ) { this->scalarData(alpha); }
00177     Scalar alpha() const { return this->scalarData(); }
00179     SUNDIALS_VScale(const Scalar& alpha)
00180       : RTOpT<Scalar>("SUNDIALS_VScale"), ROpScalarTransformationBase<Scalar>(alpha) 
00181     {}
00185     void apply_op(const int num_vecs, 
00186                   const ConstSubVectorView<Scalar> sub_vecs[],
00187                   const int  num_targ_vecs,  
00188                   const SubVectorView<Scalar> targ_sub_vecs[],
00189                   ReductTarget *reduct_obj) const
00190     {
00191       RTOP_APPLY_OP_1_1(num_vecs,sub_vecs,num_targ_vecs,targ_sub_vecs);
00192       Scalar a = alpha();
00193       for( index_type i = 0; i < subDim; 
00194            ++i,  v0_val += v0_s,  z0_val += z0_s ) 
00195         {
00196           *z0_val = a * (*v0_val);
00197         }
00198     }
00200   }; 
00201 
00202 
00203 
00207   template<class Scalar>
00208   class SUNDIALS_VAddConst : public ROpScalarTransformationBase<Scalar> 
00209   {
00210   public:
00212     void alpha( const Scalar& alpha ) { this->scalarData(alpha); }
00214     Scalar alpha() const { return this->scalarData(); }
00216     SUNDIALS_VAddConst(const Scalar& alpha)
00217       : RTOpT<Scalar>("SUNDIALS_VAddConst"), ROpScalarTransformationBase<Scalar>(alpha) 
00218     {}
00222     void apply_op(const int num_vecs, 
00223                   const ConstSubVectorView<Scalar> sub_vecs[],
00224                   const int  num_targ_vecs,  
00225                   const SubVectorView<Scalar> targ_sub_vecs[],
00226                   ReductTarget *reduct_obj) const
00227     {
00228       RTOP_APPLY_OP_1_1(num_vecs,sub_vecs,num_targ_vecs,targ_sub_vecs);
00229       Scalar a = alpha();
00230       for( index_type i = 0; i < subDim; 
00231            ++i,  v0_val += v0_s,  z0_val += z0_s ) 
00232         {
00233           *z0_val = (*v0_val) + a;
00234         }
00235     }
00237   }; 
00238 
00239 
00240 
00241 
00254   template <class Scalar>
00255   class SUNDIALS_VWL2Norm : public ROpScalarReductionBase<Scalar>
00256   {
00257   public:
00259     SUNDIALS_VWL2Norm()
00260       : RTOpT<Scalar>(""),
00261         ROpScalarReductionBase<Scalar>(0.0)
00262     {;}
00263 
00264 
00265     
00266     Scalar operator()(const ReductTarget& reduct_obj ) const 
00267     { return this->getRawVal(reduct_obj); }
00270       
00272     void apply_op(
00273                   const int num_vecs,    
00274                   const ConstSubVectorView<Scalar> sub_vecs[],
00275                   const int num_targ_vecs, 
00276                   const SubVectorView<Scalar> targ_sub_vecs[],
00277                   ReductTarget *reduct_obj) const
00278     {
00279       ReductTargetScalar<Scalar> &_reduct_obj 
00280         = Teuchos::dyn_cast<ReductTargetScalar<Scalar> >(*reduct_obj); 
00281       RTOP_APPLY_OP_2_0(num_vecs,sub_vecs,num_targ_vecs,targ_sub_vecs);
00282       Scalar sum = this->getRawVal(*reduct_obj);
00283       for( index_type i = 0; i < subDim; ++i,  
00284              v0_val += v0_s,  
00285              v1_val += v1_s) 
00286         {
00287           double x = (*v0_val) * (*v1_val);
00288           sum += x*x;
00289         }
00290       _reduct_obj.set(sum);
00291       
00292     }
00294   }; 
00295 
00296 
00309   template <class Scalar>
00310   class SUNDIALS_VWrmsMaskNorm : public ROpScalarReductionBase<Scalar>
00311   {
00312   public:
00313     
00315     SUNDIALS_VWrmsMaskNorm()
00316       : RTOpT<Scalar>(""),
00317         ROpScalarReductionBase<Scalar>(0.0)
00318     {;}
00319 
00320 
00321     
00322     Scalar operator()(const ReductTarget& reduct_obj ) const 
00323     { return this->getRawVal(reduct_obj); }
00326       
00328     void apply_op(
00329                   const int num_vecs,    
00330                   const ConstSubVectorView<Scalar> sub_vecs[],
00331                   const int num_targ_vecs, 
00332                   const SubVectorView<Scalar> targ_sub_vecs[],
00333                   ReductTarget *reduct_obj) const
00334     {
00335       ReductTargetScalar<Scalar> &_reduct_obj 
00336         = Teuchos::dyn_cast<ReductTargetScalar<Scalar> >(*reduct_obj); 
00337       RTOP_APPLY_OP_3_0(num_vecs,sub_vecs,num_targ_vecs,targ_sub_vecs);
00338       Scalar sum = this->getRawVal(*reduct_obj);
00339       const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
00340       for( index_type i = 0; i < subDim; ++i,  
00341              v0_val += v0_s,  
00342              v1_val += v1_s,  
00343              v2_val += v2_s) 
00344         {
00345           if (*v2_val > zero)
00346             {
00347               double x = (*v0_val) * (*v1_val);
00348               sum += x*x;
00349             }
00350         }
00351       _reduct_obj.set(sum);
00352       
00353     }
00355   }; 
00356 
00357 
00382   template<class Scalar>
00383   class SUNDIALS_VConstrMask : public RTOpBoolReduceAndTransform<Scalar>
00384   {
00385   public:
00386     
00388     SUNDIALS_VConstrMask()
00389       : RTOpT<Scalar>(""), 
00390         RTOpBoolReduceAndTransform<Scalar>() 
00391     {}
00392     
00396     void apply_op(
00397                   const int num_vecs,    
00398                   const ConstSubVectorView<Scalar> sub_vecs[],
00399                   const int num_targ_vecs, 
00400                   const SubVectorView<Scalar> targ_sub_vecs[],
00401                   ReductTarget *reduct_obj) const
00402     {
00403       ReductTargetScalar<index_type>& _reduct_obj 
00404         = Teuchos::dyn_cast<ReductTargetScalar<index_type> >(*reduct_obj); 
00405       RTOP_APPLY_OP_2_1(num_vecs,sub_vecs,num_targ_vecs,targ_sub_vecs);
00406       const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
00407       bool feasible = _reduct_obj.get();
00408       for( index_type i = 0; i < subDim; ++i,  
00409              v0_val += v0_s,  
00410              v1_val += v1_s,  
00411              z0_val += z0_s ) 
00412         {
00413           Scalar x_i = *v0_val;
00414           Scalar c_i = *v1_val;
00415           bool ok = true;
00416           if (c_i > 1.5 && c_i < 2.5)
00417             {
00418               ok = x_i > zero;
00419             }
00420           else if (c_i > 0.5 && c_i < 1.5)
00421             {
00422               ok = x_i >= zero;
00423             }
00424           else if (c_i < -0.5 && c_i > -1.5)
00425             {
00426               ok = x_i <= zero;
00427             }
00428           else if (c_i < -1.5)
00429             {
00430               ok = x_i < zero;
00431             }
00432           else
00433             {
00434               TEST_FOR_EXCEPTION(true, runtime_error,
00435                                  "illegal constraint flag = " << c_i
00436                                  << ". Allowed values are {-2.0, -1.0, 1.0, 2.0}");
00437             }
00438 
00439           /* Set the result element to 1 if the point is infeasible. */
00440           if (ok)
00441             {
00442               *z0_val = 0.0;
00443             }
00444           else
00445             {
00446               feasible = false;
00447               *z0_val = 1.0;
00448             }
00449         }
00450       _reduct_obj.set(feasible);
00451     }
00453   }; 
00454 
00455 
00456 
00471   template <class Scalar>
00472   class SUNDIALS_VMinQuotient : public ROpScalarReductionBase<Scalar>
00473   {
00474   public:
00475     
00477     SUNDIALS_VMinQuotient()
00478       : RTOpT<Scalar>(""),
00479         ROpScalarReductionBase<Scalar>(Teuchos::ScalarTraits<Scalar>::rmax())
00480     {;}
00481 
00482 
00483     
00484     Scalar operator()(const ReductTarget& reduct_obj ) const 
00485     { return this->getRawVal(reduct_obj); }
00489     void reduce_reduct_objs(const ReductTarget& in_reduct_obj, 
00490                             ReductTarget* inout_reduct_obj) const
00491     {
00492       const Scalar in_min_ele    = this->getRawVal(in_reduct_obj);
00493       const Scalar inout_min_ele = this->getRawVal(*inout_reduct_obj);
00494       this->setRawVal( in_min_ele < inout_min_ele ? in_min_ele : inout_min_ele, inout_reduct_obj );
00495     }
00496       
00498     void apply_op(
00499                   const int num_vecs,    
00500                   const ConstSubVectorView<Scalar> sub_vecs[],
00501                   const int num_targ_vecs, 
00502                   const SubVectorView<Scalar> targ_sub_vecs[],
00503                   ReductTarget *reduct_obj) const
00504     {
00505       ReductTargetScalar<Scalar> &_reduct_obj 
00506         = Teuchos::dyn_cast<ReductTargetScalar<Scalar> >(*reduct_obj); 
00507       RTOP_APPLY_OP_2_0(num_vecs,sub_vecs,num_targ_vecs,targ_sub_vecs);
00508       Scalar min_ele = this->getRawVal(*reduct_obj);
00509       const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
00510       for( index_type i = 0; i < subDim; ++i,  
00511              v0_val += v0_s,  
00512              v1_val += v1_s) 
00513         {
00514           if (*v1_val != zero)
00515             {
00516               Scalar x = (*v0_val)/(*v1_val);
00517               if (x < min_ele) min_ele = x;
00518             }
00519         }
00520       _reduct_obj.set(min_ele);
00521       
00522     }
00524   }; 
00525 
00526 
00527 
00528 
00538   template<class Scalar>
00539   class SUNDIALS_VInvTest : public RTOpBoolReduceAndTransform<Scalar>
00540   {
00541   public:
00542     
00544     SUNDIALS_VInvTest()
00545       : RTOpT<Scalar>(""), 
00546         RTOpBoolReduceAndTransform<Scalar>() 
00547     {}
00548     
00552     void apply_op(
00553                   const int num_vecs,    
00554                   const ConstSubVectorView<Scalar> sub_vecs[],
00555                   const int num_targ_vecs, 
00556                   const SubVectorView<Scalar> targ_sub_vecs[],
00557                   ReductTarget *reduct_obj) const
00558     {
00559       ReductTargetScalar<index_type>& _reduct_obj 
00560         = Teuchos::dyn_cast<ReductTargetScalar<index_type> >(*reduct_obj); 
00561       RTOP_APPLY_OP_1_1(num_vecs,sub_vecs,num_targ_vecs,targ_sub_vecs);
00562       const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00563       const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
00564       bool pass =  _reduct_obj.get();
00565       for( index_type i = 0; i < subDim; ++i,  
00566              v0_val += v0_s,  
00567              z0_val += z0_s ) 
00568         {
00569           if (*v0_val != zero) 
00570             {
00571               *z0_val = one/(*v0_val);
00572             }
00573           else 
00574             {
00575               *z0_val = 0.0;
00576               pass = false;
00577             }
00578         }
00579       _reduct_obj.set(pass);
00580     }
00582   }; 
00583 
00584 
00585 
00594   template<class Scalar>
00595   class SUNDIALS_VCompare : public ROpScalarTransformationBase<Scalar> {
00596   public:
00598     void alpha( const Scalar& alpha ) { this->scalarData(alpha); }
00600     Scalar alpha() const { return this->scalarData(); }
00602     SUNDIALS_VCompare( const Scalar &alpha)
00603       : RTOpT<Scalar>("SUNDIALS_Compare"), 
00604         ROpScalarTransformationBase<Scalar>(alpha) 
00605     {}
00609     void apply_op(
00610                   const int num_vecs,    
00611                   const ConstSubVectorView<Scalar> sub_vecs[],
00612                   const int num_targ_vecs, 
00613                   const SubVectorView<Scalar> targ_sub_vecs[],
00614                   ReductTarget *reduct_obj) const
00615     {
00616       RTOP_APPLY_OP_1_1(num_vecs,sub_vecs,num_targ_vecs,targ_sub_vecs);
00617       const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00618       const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
00619       for( index_type i = 0; i < subDim; ++i,  v0_val += v0_s,  z0_val += z0_s ) 
00620         {
00621           if (fabs(*v0_val) >= alpha()) *z0_val = one;
00622           else *z0_val = zero;
00623         }
00624     }
00626   };
00627 
00628 }
00629 
00630 
00631 #endif 

Generated on Thu Sep 18 12:30:43 2008 for RTOp Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1