Thyra_MultiVectorStdOps.hpp

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

Generated on Tue Oct 20 12:47:27 2009 for Thyra Package Browser (Single Doxygen Collection) by doxygen 1.4.7