Thyra Package Browser (Single Doxygen Collection): Thyra

00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) 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 Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef THYRA_MULTI_VECTOR_STD_OPS_HPP
00030 #define THYRA_MULTI_VECTOR_STD_OPS_HPP
00031 
00032 #include "Thyra_MultiVectorStdOpsDecl.hpp"
00033 #include "Thyra_VectorStdOps.hpp"
00034 #include "Thyra_VectorSpaceBase.hpp"
00035 #include "Thyra_VectorStdOps.hpp"
00036 #include "Thyra_MultiVectorBase.hpp"
00037 #include "RTOpPack_ROpDotProd.hpp"
00038 #include "RTOpPack_ROpNorm1.hpp"
00039 #include "RTOpPack_ROpNormInf.hpp"
00040 #include "RTOpPack_TOpAssignScalar.hpp"
00041 #include "RTOpPack_TOpAssignVectors.hpp"
00042 #include "RTOpPack_TOpAXPY.hpp"
00043 #include "RTOpPack_TOpLinearCombination.hpp"
00044 #include "RTOpPack_TOpScaleVector.hpp"
00045 #include "Teuchos_TestForException.hpp"
00046 
00047 template<class Scalar>
00048 void Thyra::norms( const MultiVectorBase<Scalar>& V, typename Teuchos::ScalarTraits<Scalar>::magnitudeType norms[] )
00049 {
00050   typedef Teuchos::ScalarTraits<Scalar> ST;
00051   const int m = V.domain()->dim();
00052   std::vector<Scalar> prods(m);
00053   V.range()->scalarProds(V,V,&prods[0]);
00054   for( int j = 0; j < m; ++j )
00055     norms[j] = ST::magnitude(ST::squareroot(prods[j]));
00056 }
00057 
00058 template<class Scalar, class NormOp>
00059 void Thyra::reductions( const MultiVectorBase<Scalar>& V, const NormOp &op, typename Teuchos::ScalarTraits<Scalar>::magnitudeType norms[] )
00060 {
00061   int kc;
00062   const int m = V.domain()->dim();
00063   std::vector<Teuchos::RefCountPtr<RTOpPack::ReductTarget> >  rcp_op_targs(m);
00064   std::vector<RTOpPack::ReductTarget*>                        op_targs(m);
00065   for( kc = 0; kc < m; ++kc ) {
00066     rcp_op_targs[kc] = op.reduct_obj_create();
00067     op_targs[kc] = &*rcp_op_targs[kc];
00068   }
00069   const MultiVectorBase<Scalar>* multi_vecs[] = { &V,};
00070   applyOp<Scalar>(op,1,multi_vecs,0,static_cast<MultiVectorBase<Scalar>**>(NULL),&op_targs[0]);
00071   for( kc = 0; kc < m; ++kc ) {
00072     norms[kc] = op(*op_targs[kc]);
00073   }
00074 }
00075 
00076 template<class Scalar>
00077 void Thyra::dots( const MultiVectorBase<Scalar>& V1, const MultiVectorBase<Scalar>& V2, Scalar dots[] )
00078 {
00079   int kc;
00080   const int m = V1.domain()->dim();
00081   RTOpPack::ROpDotProd<Scalar> dot_op;
00082   std::vector<Teuchos::RefCountPtr<RTOpPack::ReductTarget> >  rcp_dot_targs(m);
00083   std::vector<RTOpPack::ReductTarget*>                        dot_targs(m);
00084   for( kc = 0; kc < m; ++kc ) {
00085     rcp_dot_targs[kc] = dot_op.reduct_obj_create();
00086     dot_targs[kc] = &*rcp_dot_targs[kc];
00087   }
00088   const MultiVectorBase<Scalar>* multi_vecs[] = { &V1, &V2 };
00089   applyOp<Scalar>(dot_op,2,multi_vecs,0,static_cast<MultiVectorBase<Scalar>**>(NULL),&dot_targs[0]); // Cast required by sun compiler
00090   for( kc = 0; kc < m; ++kc ) {
00091     dots[kc] = dot_op(*dot_targs[kc]);
00092   }
00093 }
00094 
00095 template<class Scalar>
00096 void Thyra::sums( const MultiVectorBase<Scalar>& V, Scalar sums[] )
00097 {
00098   int kc;
00099   const int m = V.domain()->dim();
00100   RTOpPack::ROpSum<Scalar> sum_op;
00101   std::vector<Teuchos::RefCountPtr<RTOpPack::ReductTarget> >  rcp_op_targs(m);
00102   std::vector<RTOpPack::ReductTarget*>                        op_targs(m);
00103   for( kc = 0; kc < m; ++kc ) {
00104     rcp_op_targs[kc] = sum_op.reduct_obj_create();
00105     op_targs[kc] = &*rcp_op_targs[kc];
00106   }
00107   const MultiVectorBase<Scalar>* multi_vecs[] = { &V,};
00108   applyOp<Scalar>(sum_op,1,multi_vecs,0,static_cast<MultiVectorBase<Scalar>**>(NULL),&op_targs[0]);
00109   for( kc = 0; kc < m; ++kc ) {
00110     sums[kc] = sum_op(*op_targs[kc]);
00111   }
00112 }
00113 
00114 template<class Scalar>
00115 typename Teuchos::ScalarTraits<Scalar>::magnitudeType
00116 Thyra::norm_1( const MultiVectorBase<Scalar>& V )
00117 {
00118   // Primary column-wise reduction (sum of absolute values)
00119   RTOpPack::ROpNorm1<Scalar> sum_abs_op;
00120   // Secondary reduction (max over all columns = induced norm_1 matrix  norm)
00121   RTOpPack::ROpNormInf<Scalar> max_op;
00122   // Reduction object (must be same for both sum_abs and max_targ objects)
00123   Teuchos::RefCountPtr<RTOpPack::ReductTarget>
00124     max_targ = max_op.reduct_obj_create();
00125   // Perform the reductions
00126   const MultiVectorBase<Scalar>* multi_vecs[] = { &V };
00127   applyOp<Scalar>(sum_abs_op,max_op,1,multi_vecs,0,static_cast<MultiVectorBase<Scalar>**>(NULL),&*max_targ); // Sun complier requires this cast
00128   // Return the final value
00129   return max_op(*max_targ);
00130 }
00131 
00132 template<class Scalar>
00133 void Thyra::scale( Scalar alpha, MultiVectorBase<Scalar>* V )
00134 {
00135   typedef Teuchos::ScalarTraits<Scalar> ST;
00136 #ifdef TEUCHOS_DEBUG
00137   TEST_FOR_EXCEPTION(V==NULL,std::logic_error,"assign(...), Error!");
00138 #endif
00139   if(alpha==ST::zero()) {
00140     assign( V, ST::zero() );
00141     return;
00142   }
00143   if(alpha==ST::one()) {
00144     return;
00145   }
00146   RTOpPack::TOpScaleVector<Scalar> scale_vector_op(alpha);
00147   MultiVectorBase<Scalar>* targ_multi_vecs[] = { V };
00148   applyOp<Scalar>(
00149     scale_vector_op,0,(const MultiVectorBase<Scalar>**)NULL // The SUN compiler requires these casts!
00150     ,1,targ_multi_vecs,(RTOpPack::ReductTarget**)NULL
00151     );
00152 }
00153 
00154 template<class Scalar>
00155 void Thyra::scaleUpdate( const VectorBase<Scalar>& a, const MultiVectorBase<Scalar>& U, MultiVectorBase<Scalar>* V )
00156 {
00157 #ifdef TEUCHOS_DEBUG
00158   TEST_FOR_EXCEPTION(V==NULL,std::logic_error,"update(...), Error!");
00159   bool is_compatible = U.range()->isCompatible(*a.space());
00160   TEST_FOR_EXCEPTION(
00161     !is_compatible,Exceptions::IncompatibleVectorSpaces
00162     ,"update(...), Error, U.range()->isCompatible(*a.space())==false");
00163   is_compatible = U.range()->isCompatible(*V->range());
00164   TEST_FOR_EXCEPTION(
00165     !is_compatible,Exceptions::IncompatibleVectorSpaces
00166     ,"update(...), Error, U.range()->isCompatible((V->range())==false ");
00167   is_compatible = U.domain()->isCompatible(*V->domain());
00168   TEST_FOR_EXCEPTION(
00169     !is_compatible,Exceptions::IncompatibleVectorSpaces
00170     ,"update(...), Error, U.domain().isCompatible(V->domain())==false ");
00171 #endif
00172   const int m = U.domain()->dim();
00173   for( int j = 0; j < m; ++j ) {
00174     ele_wise_prod( Scalar(1.0), a, *U.col(j), &*V->col(j) ); 
00175   }
00176 }
00177 
00178 template<class Scalar>
00179 void Thyra::assign( MultiVectorBase<Scalar>* V, Scalar alpha )
00180 {
00181 #ifdef TEUCHOS_DEBUG
00182   TEST_FOR_EXCEPTION(V==NULL,std::logic_error,"assign(...), Error!");
00183 #endif
00184   RTOpPack::TOpAssignScalar<Scalar> assign_scalar_op(alpha);
00185   MultiVectorBase<Scalar>* targ_multi_vecs[] = { V };
00186   applyOp<Scalar>(
00187     assign_scalar_op,0,(const MultiVectorBase<Scalar>**)NULL // The SUN compiler requires these casts!
00188     ,1,targ_multi_vecs,(RTOpPack::ReductTarget**)NULL
00189     );
00190 }
00191 
00192 template<class Scalar>
00193 void Thyra::assign( MultiVectorBase<Scalar>* V, const MultiVectorBase<Scalar>& U )
00194 {
00195 #ifdef TEUCHOS_DEBUG
00196   TEST_FOR_EXCEPTION(V==NULL,std::logic_error,"assign(...), Error!");
00197 #endif
00198   RTOpPack::TOpAssignVectors<Scalar> assign_vectors_op;
00199   const MultiVectorBase<Scalar>* multi_vecs[]      = { &U };
00200   MultiVectorBase<Scalar>*       targ_multi_vecs[] = { V   };
00201   applyOp<Scalar>(
00202     assign_vectors_op,1,multi_vecs,1,targ_multi_vecs
00203     ,(RTOpPack::ReductTarget**)NULL // The SUN compiler requires this cast!
00204     );
00205 }
00206 
00207 template<class Scalar>
00208 void Thyra::update( Scalar alpha, const MultiVectorBase<Scalar>& U, MultiVectorBase<Scalar>* V )
00209 {
00210 #ifdef TEUCHOS_DEBUG
00211   TEST_FOR_EXCEPTION(V==NULL,std::logic_error,"update(...), Error!");
00212 #endif
00213   RTOpPack::TOpAXPY<Scalar> axpy_op(alpha);
00214   const MultiVectorBase<Scalar>* multi_vecs[]       = { &U };
00215   MultiVectorBase<Scalar>*       targ_multi_vecs[]  = { V  };
00216   applyOp<Scalar>(axpy_op,1,multi_vecs,1,targ_multi_vecs,(RTOpPack::ReductTarget**)NULL); // Sun compiler requires this cast
00217 }
00218 
00219 template<class Scalar>
00220 void Thyra::update( Scalar alpha[], Scalar beta, const MultiVectorBase<Scalar>& U, MultiVectorBase<Scalar>* V )
00221 {
00222 #ifdef TEUCHOS_DEBUG
00223   TEST_FOR_EXCEPTION(V==NULL,std::logic_error,"update(...), Error!");
00224   bool is_compatible = U.range()->isCompatible(*V->range());
00225   TEST_FOR_EXCEPTION(
00226     !is_compatible,Exceptions::IncompatibleVectorSpaces
00227     ,"update(...), Error, U.range()->isCompatible((V->range())==false ");
00228   is_compatible = U.domain()->isCompatible(*V->domain());
00229   TEST_FOR_EXCEPTION(
00230     !is_compatible,Exceptions::IncompatibleVectorSpaces
00231     ,"update(...), Error, U.domain().isCompatible(V->domain())==false ");
00232 #endif
00233   const int m = U.domain()->dim();
00234   for( int j = 0; j < m; ++j )
00235     Vp_StV( V->col(j).get(), alpha[j]*beta, *U.col(j) );
00236 }
00237 
00238 
00239 template<class Scalar>
00240 void Thyra::update( const MultiVectorBase<Scalar>& U, Scalar alpha[], Scalar beta, MultiVectorBase<Scalar>* V )
00241 {
00242 #ifdef TEUCHOS_DEBUG
00243     TEST_FOR_EXCEPTION(V==NULL,std::logic_error,"update(...), Error!");
00244   bool is_compatible = U.range()->isCompatible(*V->range());
00245     TEST_FOR_EXCEPTION(
00246     !is_compatible,Exceptions::IncompatibleVectorSpaces
00247     ,"update(...), Error, U.range()->isCompatible((V->range())==false ");
00248   is_compatible = U.domain()->isCompatible(*V->domain());
00249     TEST_FOR_EXCEPTION(
00250     !is_compatible,Exceptions::IncompatibleVectorSpaces
00251     ,"update(...), Error, U.domain().isCompatible(V->domain())==false ");
00252 #endif
00253   const int m = U.domain()->dim();
00254   for( int j = 0; j < m; ++j ) {
00255     Vt_S( V->col(j).get(), alpha[j]*beta );
00256     Vp_StV( V->col(j).get(), 1.0, *U.col(j) );
00257   }
00258 }
00259 
00260 template<class Scalar>
00261 void Thyra::linear_combination(
00262   const int                         m
00263   ,const Scalar                     alpha[]
00264   ,const MultiVectorBase<Scalar>*   X[]
00265   ,const Scalar                     &beta
00266   ,MultiVectorBase<Scalar>          *Y
00267   )
00268 {
00269 #ifdef TEUCHOS_DEBUG
00270   TEST_FOR_EXCEPTION(Y==NULL,std::logic_error,"linear_combination(...), Error!");
00271 #endif
00272   if( beta == Teuchos::ScalarTraits<Scalar>::one() && m == 1 ) {
00273     update( alpha[0], *X[0], Y );
00274     return;
00275   }
00276   else if( m == 0 ) {
00277     scale( beta, Y );
00278     return;
00279   }
00280   RTOpPack::TOpLinearCombination<Scalar> lin_comb_op(m,alpha,beta);
00281   MultiVectorBase<Scalar>* targ_multi_vecs[] = { Y };
00282   Thyra::applyOp<Scalar>(lin_comb_op,m,X,1,targ_multi_vecs,(RTOpPack::ReductTarget**)NULL);  // Cast returned by sun compiler
00283 }
00284 
00285 template<class Scalar>
00286 void Thyra::randomize( Scalar l, Scalar u, MultiVectorBase<Scalar>* V )
00287 {
00288 #ifdef TEUCHOS_DEBUG
00289     TEST_FOR_EXCEPTION(V==NULL,std::logic_error,"randomize(...), Error!");
00290 #endif
00291   const int m = V->domain()->dim();
00292   for( int j = 0; j < m; ++j ) {
00293     randomize( l, u, V->col(j).get() ); // Todo: call applyOp(...) directly!
00294   }
00295 }
00296 
00297 template<class Scalar>
00298 void Thyra::V_VpV( MultiVectorBase<Scalar>* Z, const MultiVectorBase<Scalar>& X, const MultiVectorBase<Scalar>& Y )
00299 {
00300   typedef Teuchos::ScalarTraits<Scalar> ST;
00301   linear_combination(
00302     2,Teuchos::arrayArg<Scalar>(ST::one(),ST::one())()
00303     ,Teuchos::arrayArg<const MultiVectorBase<Scalar>*>(&X,&Y)()
00304     ,ST::zero(),Z
00305     );
00306 }
00307 
00308 template<class Scalar>
00309 void Thyra::V_VmV( MultiVectorBase<Scalar>* Z, const MultiVectorBase<Scalar>& X, const MultiVectorBase<Scalar>& Y )
00310 {
00311   typedef Teuchos::ScalarTraits<Scalar> ST;
00312   linear_combination(
00313     2,Teuchos::arrayArg<Scalar>(ST::one(),Scalar(-ST::one()))()
00314     ,Teuchos::arrayArg<const MultiVectorBase<Scalar>*>(&X,&Y)()
00315     ,ST::zero(),Z
00316     );
00317 }
00318 
00319 #endif // THYRA_MULTI_VECTOR_STD_OPS_HPP
Thyra_MultiVectorStdOps.hpp
